Bassi

E Studer

A Beccantini

S Gounand

F Dabbene

J P Magnaud

H Paillère

I Limaiem

F Damian

H Golfier

C Bassi

J C Garnier

CAST3M/ARCTURUS: A coupled heat transfer CFD code for thermal-hydraulic analyzes of gas cooled reactors

CAST3M/ARCTURUS: A coupled heat transfer CFD

code for thermal-hydraulic analyzes of gas cooled reactors

Introduction

The safety of gas cooled reactors must be ensured by active or passive systems which fulfill the task of keeping loads on component (fuel) and structures (vessel, containment) within acceptable limits under accidental conditions at all moments. To achieve this objective, thermal-hydraulics computer codes are necessary tools to design, enhance the performance and ensure a high safety level of the different reactors. ARC-TURUS is a multi-dimensional (2D and 3D) heat transfer and thermal-hydraulics simulation tool developed by CEA for appli-cations to gas cooled reactors. Models, elementary test cases, benchmark exercises and applications are developed and consolidated in the CAST3M platform environment, which is a general purpose structural and fluid mechanics code. The first version of ARCTURUS was released 2 years ago [START_REF] Studer | Gas cooled reactor thermal-hydraulics using CAST3Mand CRONOS2 codes[END_REF] and was applied to the High Temperature Engineering Reactor (HTTR) Vessel Cooling benchmark and MHTGR steady state calculations. Some new developments, verifications and applications are reported in the present article.

The first section is dedicated to GFR containment analyses of a depressurization accident. This includes the development of a lumped parameter thermal-hydraulics analyses using a depressurization accident scenario coming from CATHARE code simulations. The second section focuses on GCR invessel thermal-hydraulics including the description of the turbulent contribution 1, 2, 3 location in the simplified depressurization model ARCTURUS "equivalent porous media" and 1D fluid/3D thermal models, elementary verification of these models and application to the HTR-10 benchmark exercise. In the third section, we describe some coupling analyses between thermal-hydraulics and neutronics for HTR reactors in order to investigate the importance of the various parameters in feedback calculations. Conclusions follow.

GFR Containment thermal-hydraulics

Objectives

Current concepts of Gas Cooled Fast Reactors use helium gas as coolant under about 70 bars of total pressure. Reactor vessel depressurization accident is one of the accidents that has to be considered to design the reactor core and the safety systems. The advantages of helium gas drastically disappear when the pressure decreases. So, as far as passive residual heat removal is concerned, helium is not the most interesting gas unless a sufficient pressure is maintained within the reactor vessel (including the guard containment) or a heavier gas is added. A leak-tight containment is first chosen to prevent any off-site radioactive release in post-accidental situations. An important safety issue is to properly justify the containment (pre)design, taking into account the thermal and mechanical loadings which could occur in postulated depressurization accidents. In previous exploratory studies [START_REF] Poette | Advanced GCFR preliminary design 300 MWe -Project status and trends for a higher unit power selection[END_REF][START_REF] Dumaz | The thermal-hydraulic studies in support to the GFR pre-conceptual design[END_REF], the containment also had the safety function described above which is to ensure, in depressurization accidents, a back-up pressure high enough to cool the core with pure gas natural circulation (guard containment option for passive decay heat removal). Therefore, the guard containment consists in a pre-stressed concrete structure -single-cavity cylinder, as close as possible to the primary circuit. It is (pre)designed to sustain a static pressure of about 2 MPa, whereas, in normal operating conditions, it is pressurized with dry air at about 1.4 MPa (therefore allowing for potential pressure increase in accidental conditions). The objective of the present study is two-fold:

• to prepare and test a preliminary model of the thermalhydraulics phenomena in the guard containment looking at heat transfer to the walls and structures and pressure evolution versus time, • to consolidate the design presented in [START_REF] Poette | Advanced GCFR preliminary design 300 MWe -Project status and trends for a higher unit power selection[END_REF] by performing some scoping studies.

The studies are conducted with the ARCTURUS code using boundary conditions (mass and energy flow rates into the containment) coming from CATHARE system code calculations. In the first part, the design of the 600 MW GFR containment including free volumes, walls area, flow paths between different volumes is described. These data are a necessary part of the containment code input deck and the lumped parameter approach adopted for the ARCTURUS code is also briefly described in the same section. Then, the second part is devoted to an analysis of a CATHARE depressurization scenario (1 in. break Loss-Of-Coolant Accident (LOCA) cumulated with station black-out).

Containment geometry and lumped parameter model description

A general presentation of the 600 MW GFR reactor is presented in Fig. 1. The containment is mainly divided into four main compartments: the reactor pit with the reactor vessel inside, the compartment above the reactor pit with the control rods mechanisms, the Power Conversion System (PCS) compartment including mainly the gas turbine, and a surrounding compartment with the Residual Heat Removal System (RHRS composed of three heat exchangers). These fluid volumes are enclosed by concrete walls and the different volumes are connected through atmospheric junctions (doors, stairs, . . .). All this system is included in a reinforced concrete containment of about 14000 m 3 . In Table 1, the main characteristics of this containment are reported including free gaseous and heat slabs volumes. The ARCTURUS Lumped Parameter (LP) approach of the reactor containment fluid mechanics is based on mass and energy balance in sub-compartment and an inertia type equation is solved to connect together the different sub-compartments [START_REF] Studer | On the use of the MISTRA coupled effect test facility for the validation of containment thermal-hydraulics codes[END_REF].

V k dρ l k dt + j ρ l kj ρ kj Q kj = 0 (1) ρ k - l ρ l k = 0 (2) L kj S kj dQ kj dt + (p j -p k ) + K kj |Q kj |Q kj 2ρ kj S 2 kj -ρ k g(z k -z kj ) +ρ j g(z j -z kj ) = 0 (3) dρ k V k e k dt + j [h kj Q kj ] = 0 (4)
where V k represents the volume of compartment k, ρ l kj the density of species l in the junction between compartment k and j, ρ kj the total mass density in this junction, Q kj the total mass flow rate through the junction and ρ l k the density of species l in the compartment k. In the momentum equation, K kj represents a friction coefficient in the junction between compartment k and j, of area S kj and length L kj , and z is the vertical coordinate.

A free volume enclosed in concrete walls may be divided in several sub-compartments in order to describe the flow pattern inside the containment. Usually, the sub-grid of the lumped parameter approach is constructed by engineering judgment and often verified by some CFD simulations. In the present study, this last item has not been performed because we are at the design stage and many hypotheses have already been taken into account due to the lack of information regarding the design drawings. Nevertheless, the experience gained by simulating Pressurized Water Reactor (PWR) containment thermal-hydraulics [START_REF] Studer | On the use of the MISTRA coupled effect test facility for the validation of containment thermal-hydraulics codes[END_REF] can be applied to GFR containment. A schematic view of the proposed nodalization of the GFR 600 MW containment is given in Figs. 2 and3. The reactor cavity and the PCS containment are divided in two radial sectors and five axial levels. The three axial levels between +8 m and +24 m for the cavity or +27.5 m for the PCS compartment are composed of an inner zone close to the vessels and an outer zone close to the concrete walls. The RHRS compartment is also composed of five axial levels but four radial sectors (only 3 have a RHRS heat exchanger). The control rods compartment is discretized with a single gaseous volume. Finally, 45 compartments are used to describe the gaseous volume of the containment. Four different kinds of heat slabs are considered in the present model:

• the steel plugs;

• the concrete inner walls;

• the component walls (reactor vessel, PCS, RHRS exchangers and pipes) used as thermal boundary conditions; • the outer containment walls composed of three layers: a steel liner of 6 mm, an air gap of 1 mm and finally the external concrete wall. Detailed distributions of the different walls between the 45 compartments are not reported in this article but this gives finally 118 walls. Connections between the different compartments are modeled by atmospheric junctions. In the proposed nodalization two situations may be identified. The first refers to a virtual nodalization of an open volume such as the reactor cavity or Fig. 3. GFR 600 MW-Proposed nodalization: upper view. the power conversion system compartment. For these, all the possible connections have been taken into account and the junction surface is taken as the real open surface, the length as the distance between the center of each compartment and the friction coefficient as 0.5 in vertical directions and 1.0 in horizontal directions. These numbers have been determined during verification processes (NUPEC M2-2 test and also MISTRA tests [START_REF] Studer | On the use of the MISTRA coupled effect test facility for the validation of containment thermal-hydraulics codes[END_REF]). The second situation refers to real openings between compartments: doors, gratings in floors, holes for pipes, . . .. Regarding this latter point, the proposed design is not sufficient to characterize such geometrical details. So, it has been assumed that a limited number of connections exists between the different compartments and compared to French PWR containment a more confined geometry is adopted. Finally, 90 atmospheric junctions are used for the flow pattern nodalization.

ARCTURUS/CATHARE depressurization scenario

The fast depressurization of the core vessel occurring after a LOCA in a GCR study is a complex problem. Indeed, complex shock structures traveling and reflecting on the external containment may occur, threatening the integrity of the reactor systems. Detailed CFD simulations are needed for such analysis and this has been performed using compressible flow models of the ARCTURUS code for the safety evaluation of high pressure experimental loops. Here, we describe and use a simplified model to determine whether the average increase of pressure and temperature can threaten the integrity of the containment. The selected scenario corresponds to a small break LOCA (1 in. equivalent diameter) cumulated with station blackout. This depressurization accident scenario has been computed with the system code CATHARE [START_REF] Dumaz | The thermal-hydraulic studies in support to the GFR pre-conceptual design[END_REF] using constant conditions for the containment (pure helium at 10 bar total pressure). The main hypothesis regarding the accident scenario are the following:

• 2.5 cm small break LOCA in the cross-duct cumulated with station black-out. • Safety core/turbine by-pass valve is in operation.

• Natural convection loop with a helium/CO2 heat exchanger (HEX) which comes into operation when the turbine velocity goes to zero. • Inertia of the turbomachinary is taken into account and so, 10 bar of back-up pressure are estimated to be sufficient for residual heat removal by natural convection (preliminary study).

The main time schedule of the depressurization accident is:

• t 0 : reactor shut-down;

• t 0 + 10 s: start of reduction of PCS rotation velocity;

• t 0 + 54 s: stop of the PCS rotation;

• t 0 + 64 s: He/CO 2 HEX in operation;

• t 0 + 6180 s: equilibrium P vessel = 10 bar. and the mass and enthalpy at the break are plotted in Figs. 4 and5. The injected energy during 10,000 s corresponds to about 9.9 GJ. The main containment thermal-hydraulic results are reported in Table 2 and Fig. 6. Due to slow transient (compared to large break LOCA), the calculated values are far from the adiabatic ones. Part of the energy is absorbed by the walls. Indeed, the maximum overpressure is computed at the end of the sonic phase and the maximum gas temperature in the break zone is reached earlier mainly due to the heat losses to the walls. Regarding gas mixing, gradients are maintained until the injected mass flow rate decreases (at about 6000 s). Thermal loads to the steel liner are very small for such scenario. Nevertheless, control volumes used in LP codes lead to homogeneisation and local effect may be more challenging (CFD simulation needed). There is also another important point in the model used in the present simulation: the heat transfer correlation used by the computer code (key parameter [START_REF] Macnab | The CONTEMPT-G computer program and its application to HTGR containments[END_REF]). We use a turbulent natural convection correlation along vertical plate (Nu = a(Gr.Pr) 1/3 ). Such correlation, free of length scale (due to Gr 1/3 ) is suitable for a LP approach. Nevertheless, in a fast depressurization accident, forced convection (or at least mixed convection) occurs near the injection zone. Thus, the heat transfer is increased and available correlation introduces a Reynolds number (an additional length scale and velocity scale). Selection of such "user defined" parameters is not straightforward in a LP approach and this can be performed by several ways (height of the compartment, weighted velocities in the atmospheric junctions, . . .). Indeed, a thorough validation process is needed to implement such extension and this has not been performed in the present study.

The next step will the feed-back of these containment thermal-hydraulic results on the CATHARE simulation (transient behavior of the back-up pressure, the gas temperature and the mixture composition). This last example concludes the containment thermal-hydraulic part. Pressure, gas temperature and concentration transient behaviors have been estimated for the guard containment. These results will be used first to check the possible modification of the in-vessel transient (onset of natural convection, possible air ingress, . . .) and secondly to evaluate the structural behavior of the containment.

GCR in-vessel thermal-hydraulics

Objectives

This section is dedicated to in-vessel thermal hydraulics studies with ARCTURUS/CAST3M code for both HTR and GFR reactors. Different thermal hydraulic models are available in the code including 3D thermal/1D hydraulic and 3D equivalent porous media formulation. After a short description of these two formulations, a cross-comparison is detailed on a simplified test case, a single heated channel under forced convection (GFR fuel assembly). Then, some recent results obtained with the equivalent porous media formulation at the reactor scale are reported (HTR-10 thermal hydraulics benchmark).

ARCTURUS models description

Equivalent porous media formulation

The asymptotic low Mach number flow model using equivalent porous media formulation developed for the ARCTURUS code solves the following set of equations: let us first introduce the porous media velocity

u D = β u (5)
where β represents the porosity (ratio of fluid volume to total volume) and u the velocity. The dynamic pressure p (P(x, t) = P(t) + p (x, t)) is recovered from the velocity divergence constrain:

∇ • u D = - 1 ρ f β ∂ρ f ∂t + u D • ∇ρ f (6)
where ρ f represents the fluid density. The velocity field u D is obtained from the momentum equation:

ρ f β ∂ u D ∂t + ρ f β 2 u D • ∇ u D = F + ρ f g -∇p + μ + μ T β u D (7)
where F represents the frictional pressure loss, g the gravity and μ T the turbulent viscosity. The fluid temperature T f comes from the fluid energy balance:

ρ f Cp f β ∂T f ∂t + ρ f Cp f u D • ∇T f = β∇ • (λ + λ T ) ∇T f + W fs (8)
where λ T represents the turbulent heat conductivity and W fs the heat flux inside the porous media determined by correlations for the heat exchange coefficient. Then, the solid temperature T s is expressed from the solid energy balance:

(

1 -β)ρ s Cp s ∂T s ∂t = λ eq T s -W fs + Φ rad + Φ neutro (9)
where λ eq represents the equivalent thermal conductivity of the homogenized porous media, Φ rad the power exchanged by thermal radiation and Φ neutro the nuclear power. Finally, the equation of state is used with constant thermodynamic pressure (steadystate) to recover the fluid density ρ f :

ρ f T f = constant (10)
The spatial discretization uses finite elements (Q2P1: quadratic for velocity and linear for pressure) and the numerical method uses semi implicit incremental projection method [START_REF] Studer | On the use of the MISTRA coupled effect test facility for the validation of containment thermal-hydraulics codes[END_REF]. The linear system is solved by using classical iterative methods. Thermal radiation is based on classical grey surfaces emissivity in enclosures methods available in the CAST3M code.

3D thermal/1D hydraulic formulation

A simplified steady state model has been implemented in the CAST3M/ARCTURUS code in order to study the thermal behavior of a vertical fuel assembly under normal conditions or in case of accidental blocked channel. In this model, the thermal behavior is supposed to be 3D and the flow is assumed 1D. The thermal coupling between the both is performed via correlations. Thermal radiation and pressure losses are neglected in this geometry configuration. With the above assumptions, the system of equation reduces to:

• energy conservation in the solid part

∇ • (-λ s ∇T s ) = W s (11)
where λ s is the homogenized fuel heat conductivity, T s the fuel temperature and W s is the power generated in the reactor core.

• mass conservation in the fluid region

S channel ∂ρ f w ∂z = 0 ( 12
)
where z is the vertical coordinate, ρ f the coolant density, w the coolant vertical velocity and S channel is the fluid channel cross section. • energy conservation in the fluid region

S channel ∂ρ f Cp f wT f ∂z = S channel ∂ ∂z λ f ∂T f ∂z + W f ( 13
)
where Cp f is the coolant heat capacity, T f the coolant temperature, λ f the coolant heat conductivity and W f the volumetric heat source described below.

Regarding boundary conditions, the coolant temperature (T f ) and the coolant mass flow rate (ρ f w) at the inlet are imposed, zero flux conditions are supposed on the outer solid surface and for the surface between coolant and solid

-λ s ∂T s ∂n = h(T s -T f ) ( 14
)
where n is the normal vector and h a heat exchange coefficient described below. The coolant properties including ρ f are supposed to be temperature-dependent [START_REF] Melese | Thermal and Flow design of Helium Cooled Reactors[END_REF].

Coupling between fuel and coolant channels

From the Colburn correlation,

Nu = 0.023Re 0.8 Pr 1/3 (15)
Nu = hd H λ f (16) Re = ρ f wd H μ f (17) Pr = μ f Cp f λ f (18)
these numbers are evaluated at a mean temperature:

T m = T s + T f 2 (19)
one can deduce the heat exchange coefficient h. Then, the volumetric heat source W f can easily be expressed from the surface integration of ( 14) at a given elevation z 0 .

Elementary test case: GFR fuel assembly under normal conditions

Comparison of the two previous models on a simplified test case is performed in this section. The test case correspond to GFR prismatic fuel element with the main characteristics defined in Table 3. Additional hypotheses are assumed:

• a cosine shape of the power density is selected with the same axial peaking factor (1.5), • no irradiation effect has been taken into account on the thermal conductivity of the fuel element, • helium is supposed to behave as a perfect gas.

The maximum temperatures are compared in Table 4. A difference of about 90 • C is computed between the two models. If one compares the homogenized temperature of the 1D-3D model (surface homogenization at each axial level), the difference decreases to about 40 • C (Fig. 7). The main reason for that difference is the non azimuthal symmetry of the 1D-3D model (Fig. 8). The maximum fuel temperature is located at the same location for the both models. The computed pressure drop through the reactor core is about 0.15 bar (1.95 m in height) in the equivalent porous media model and this is comparable to the value estimated by classical correlations. By using the friction coefficient correlation (Blasius type)

f = 0.316Re -0.25 (20)
the pressure drops can be estimated by a simplified momentum balance (inlet and outlet pressure drops are neglected):

p out -p in = [ ρf -ρ f,m ]gL -[ρ f,out u 2 out -ρ f,in u 2 in ] - 1 2 ρ f,m f m L d H u 2 m ( 21
)
where the subscript 'in' refers to the inlet conditions, 'out' to the outlet conditions, ρ f,m corresponds to the mean value of the density in the reactor core and ρf to the external conditions (it is assumed regarding the plenum volumes that the external condition is composed of 80% of the hot helium and 20% of the cold helium). This gives the following results for the different terms of the pressure difference:

• gravity term: 6 Pa,

• inertia term: 3976 Pa,

• friction term: 11,669 Pa.

As expected, the friction term represents the main part of the pressure difference under forced convection conditions. So, the total pressure difference is about 0.156 bar. Finally, we can conclude that our equivalent porous media gives accurate results and has been checked under forced convection conditions. Verification of this model is the subject of the next section.

HTR-10 benchmark exercise

HTR-10 is a pebble bed High Temperature gas cooled reactor designed, constructed and operated by the Institute of Nuclear Energy Technology (INET) in China. Full power operation with a core outlet temperature at 700 • C was achieved in January 2003. This reactor operates with a primary helium pressure of about 3 MPa and a core inlet temperature of 250 • C. The benchmark proposed within the IAEA/CRP5 framework [START_REF] Dong | Benchmark problem of the HTR-10 steady state temperature distribution for the full power initial core[END_REF] is devoted to steady state temperature distribution inside the reactor vessel for Full Power Initial Core (FPIC). It has been proposed both for code-to-code and for code-toexperiment benchmarking. Thermocouples have been installed in the reactor cavity to provide the experimental results. CEA contributes to this benchmark exercise with the ARCTURUS CFD code (equivalent porous media formulation). A vertical cross-section of the HTR-10 Reactor Pressure Vessel (RPV) is described in depth and plotted in Fig. 9. The active core is cylinder-shaped and cone-shaped at the bottom to unload fuel elements. It is surrounded by graphite reflectors categorized as top, side and bottom reflectors. The active core and the reflectors are enclosed in a core vessel radially and upper and lower core support plates axially. Stagnant helium between this core vessel and the reactor pressure vessel acts as insulation and most of the heat is transfered by thermal radiation. After entering the RPV, the cold helium flows downwards through the annular gap between the core vessel and the RPV. Then it changes its flow direction to flow upwards. A small part of this cold helium goes through the cone-shaped discharging tube to cool the fuel elements and after merges with the hot helium. The major part of helium goes around the core support structures and enters into the cold helium channels in the side reflector. At the top of this reflector, cold helium is collected in the upper plenum and then for a second time, it changes its flow direction. Part of this helium goes through the control rods channels and is collected in a small plenum above. The main stream of helium flows through the active core and is collected in the lower hot plenum (Fig. 10). Correlations (friction pressure drops, heat transfer in the pebble bed, . . .), delivered nuclear power distribution and material properties (solids and helium coolant) have been provided in the specification document [START_REF] Dong | Benchmark problem of the HTR-10 steady state temperature distribution for the full power initial core[END_REF]. Constitutive relations are very important but are not discussed in the present article because there have been provided in the benchmark specification document. Illustrations of the ARCTURUS results compared to the provided experimental results [START_REF] Dong | Additionnal Informations to Benchmark problem of the HTR-10 steady state temperature distribution for the full power initial core[END_REF]) are given below. Interpretation of the differences are still under discussion and so, only a few comments are reported. Three radial solid temperature profiles have been selected for comparison at Z = 0 cm, Z = 80 cm and Z = 170 cm in Fig. 11 (the reference point (0,0) is located in the upper center of the pebble bed). First, the computed results do not start at 0.0 m radius because of symmetry axis translation needed by quadratic finite element discretization. Second, the Fig. 12. HTR-10-CAST3M/ARCTURUS mesh description. computed profile at Z = 170 cm is not reported in the active core due to grid deformation (see Fig. 12). The comparison with the experimental results shows a rather poor comparison especially in the middle core axial level. Nevertheless, it seems that the slopes are rather well predicted (correctly computed conduction fluxes). One important reason for that deviation could be the thermal effect of the by-pass through the control rods. In this region, a by-pass of 2% of the total mass flow rate creates a laminar flow through the control rods channels and in the computation, the heat extracted by this convective transfer may be perhaps too low compared to the reactor case. In the hot helium plenum (Table 5), the calculated temperatures are lower than the measured ones (about 20 • C) and the radial difference is adequate. One has to remember that the computed values correspond to homogenized values of the solid temperature and no reconstruction is proposed in the present study to recover the solid temperature. A 2D view of the computed solid temperature and velocity fields in the core region is given in Fig. 13. This benchmark exercise has demonstrated the capabilities of CAST3M/ARCTURUS code to simulate steady state temperature fields including complex convective, conductive and radiative heat transfers. Some dis- crepancies between the calculated and the measured values still exist especially where convective heat transfer is important (bypass regions). Some adjustment may be done after discussions of the present results with the experimental team.

Coupling between thermal-hydraulics and neutronics in HTR

One of the principal objectives of CEA's R&D on HTR is the development of a global simulation system. This model will be used for the design of different reactor configurations, the evaluation of their performances and finally the optimization of the different design parameters. Coupling thermal-hydraulic with the neutronic modeling allows to determine the impact of the temperature feed-back on the neutronic parameters such as cross-sections, power distribution but also the temperature coefficient (Doppler, moderator), critical position of control rods [START_REF] Yadigaroglu | Trends and needs in experimentation and numerical simulation for LWR safety[END_REF][START_REF] Rademer | Coupling of neutronics and thermal-Hydraulics codes for the simulation of transients of pebble bed HTR reactors[END_REF][START_REF] Reitsma | The PBMR steadystate and coupled kinetics core thermal-hydraulics benchmark test problems[END_REF]. A preliminary version of neutronic and thermal-hydraulic coupling system, based on a generic thermal-hydraulic model [START_REF] Studer | Gas cooled reactor thermal-hydraulics using CAST3Mand CRONOS2 codes[END_REF], and v0 neutronic calculation scheme [START_REF] Limaiem | VHTR core modelling: Coupling between neutronics and thermal-hydraulics[END_REF] was developed in the past. In the next paragraph we present the new coupling system. This new simulation system is based on the one hand on an improved thermal-hydraulic model (described in the previous paragraphs) and on the other hand on the v2 neutronic calculation scheme. In this section, first, the new coupling system is detailed. Then, some results of the sensitivity study are discussed. And finally an example of coupling is detailed.

Coupling system architecture

The architecture of the coupling system is built on the combination of two computer codes: CRONOS2 for the core neutronic calculation and CAST3M/ARCTURUS for the core thermal-hydraulics calculation (Fig. 14). Linux PIPE system files supervise the data exchange between the different calculations. The coupling system is very complex. In fact, it depends on many numerical and physical parameters [START_REF] Limaiem | VHTR core modelling: Coupling between neutronics and thermal-hydraulics[END_REF]. In order to evaluate the importance of each parameter and to reduce the system dependency on those parameters, we carry out a sensitivity study to evaluate the weight of each parameter and the accuracy of system.

Coupling system parameter sensitivity study

The sensitivity study involves many numerical and physical parameters. An exhaustive review of the system variables has been established. Only the main parameters are preserved for Results dealing only with radial and axial thermal-hydraulic meshing are detailed below. The geometry used for this study relates to an annular prismatic blocks core type. For each sensitivity study a reference configuration is defined. It allows comparison between different test cases. We choose this method since at present time there are no available thermal and neutronic experimental results for HTR reactor. A reference case is characterized by a fine meshing and optimal set of input parameters for both neutronic and thermal-hydraulic calculations.

Radial meshing study for the thermal-hydraulics modeling

We define a reference thermal-hydraulics radial meshing described in Fig. 15. The same radial meshing as in the neutronic modeling is used [START_REF] Limaiem | VHTR core modelling: Coupling between neutronics and thermal-hydraulics[END_REF]. Then, we compare these results to a coarse grid configuration described in Fig. 16. In the reference configuration, there are 36 meshes in the standard fuel and reflector assemblies (61 in the control assemblies). The number of meshes in the test configuration is reduced to 6 in the standard fuel and reflector assemblies (32 in the control assemblies). On the one hand, for the reference case, we exchange 54 power values and 36 temperatures per assembly. On the other hand, we exchange only 6 power values and 6 temperature values in the test configuration. The number of temperature values exchanged between the neutronic and the thermal-hydraulics calculation is directly related to the number of depletion regions defined in the neutronic modeling. The temperature and the power values exchanged between codes are calculated by homogenization procedures. The comparison between reference and test case is made on both detailed and integral neutronic and thermal-hydraulics results of the coupling calculation such as:

• neutronic eigenvalue, • 1D/2D/3D power and temperature distribution, • power axial offset, • energy distribution in the core.

For all coupling calculations made in the sensitivity study, all neutronic modeling parameters are fixed so we could evaluate properly the influence of the thermal-hydraulic meshing. The comparison of the 2D power and thermal distributions between the reference (fine mesh) and the test configuration (N1 -coarse mesh) shows variations under 1% of the power for 80% of the meshes in Fig. 17. For the 2D thermal distribution variations are under 2 K on 70% of the meshes in Fig. 17. For all the investigated test configurations, variations never exceed 3% of the power and 30 K. These differences are located in regions characterized by an important power and temperature gradient (between fuel and reflector assemblies). If we only take into account the coupling calculation CPU time of each configuration, we could choose test configuration N1 of thermal-hydraulic meshing for the 3D calculation. Indeed, the reference configuration takes twice as much CPU time calculation of test configuration N1 (thermal-hydraulic part). Moreover, the dis-crepancies observed on the neutronic power are under 1% for 80% of the calculation meshes. This value of 1% is acceptable since we define this value as convergence criteria for the power distribution during the coupling calculation. In other words, variations less than 1% have no interest because the precision in the power distribution is 1%. The temperature variation concerning the test configuration N1 is also acceptable since we are under 2 K on 70% of the meshes.

Axial meshing study for the thermal-hydraulics modeling

We use the same method as the one defined for the radial mesh. We select an axial meshing reference for thermal-hydraulics modeling in Fig. 18. Then we made comparison to coarser grid configurations in Fig. 19 (we present only one test case). All thermal-hydraulics configurations are coupled with a reference neutronic calculation modeling. The results of the comparison between the reference axial meshing and the meshing test configuration N2 shows differences in temperature that reaches 30 K in the fuel assemblies zone of the core in Fig. 20 (at each axial level, a 2D examination of the power and temperature field is performed). In the other part of the core temperature difference varies from -15 K to +15 K. These temperature variations cause differences in the power distribution of 5% in certain regions of the core. We demonstrate that the use of four axial calculation points per assembly in the thermal-hydraulics modeling reduces the temperature discrepancies to 5 K. In this study, the differences on power distribution are under 2%.

The radial and axial meshing study performed in the thermal-hydraulic modeling of the coupling system allows us to evaluate the impact of the meshing on the neutronic and thermal results of the coupling system. Thus, we define an optimal thermal-hydraulics model for the coupling system described as the following:

• 6 radial meshes for standard assemblies, • 32 radial meshes for control assemblies, • 5 axial meshes for standard and control assemblies, Fig. 17. Radial meshing study: N1 test configuration results compared to the reference ones. The differences are expressed in terms of percentages of the grid meshes with a certain level of power and temperature differences.

• Exchange of 6 radial power and temperature values for standard assemblies (32 for the control assemblies), • Exchange of 1 or 5 power and temperature values for standard and control assemblies.

This optimal model is used today to calculate different steady state situation of the GTM-HR core reactor.

Some neutronic and thermal-hydraulics results of the coupling system

These results concern the coupling configuration defined in the last paragraph in Fig. 20. The convergence of the coupling iterative calculation is obtained when the variations of the 3D power and temperature distributions are less than 1% for the power and 1 K for the temperature (between two successive iter-Fig. 18. Reference geometry for the axial meshing study (6 radial meshes and 5 axial meshes). ation of the coupling calculation). Only some results are reported below. The temperature inside the core varies from 480 (helium inlet temperature at the top of the core) to 1210 • C in the bottom of the core (Fig. 21). The hottest point of the core is located at 2.5 m above the bottom of the core. The thermal feedback impacts the axial power distribution. It induces a shift of the power peaking factot from the middle of the core (in a constant temperature hypothesis) to the top of the core (Fig. 21). The average helium temperature at the outlet of the core is 915 • C (480 • C at the core inlet). In steady-state conditions, the helium power in the bottom of the core is 100 MW. Ninety-four percent of the energy is taken from the fuel assemblies, 4% from control rod cavity (large cavities where control rods are inserted) and 2% from the reflector by-pass. Helium velocity varies from 23 to 48 m/s. Due to thermal dilatation, the helium velocity is more important in the hottest region of the core.

Conclusions

This paper has described some flows and heat transfer models of the CAST3M/ARCTURUS code and applications to simulate the thermal-hydraulic behavior of gas cooled reactors (containment and in-vessel). Different models are available depending on the geometrical length scale. A thorough validation process is being conducted using basic separate effect tests, in-house experiments and international benchmarks. Surface radiation and interaction between convective and radiative heat transfer are important phenomena that have to be accurately modeled to quantify the different temperature levels or extracted power suitable to design, enhance the performance and ensure a high safety level of the different reactors. Containment thermal-hydraulics following a depressurization accident is also important to quantify pressure and thermal loads on components (guard containment) and to provide transient back-up pressure for the onset of natural convection inside the vessel (GFR concepts). This last point is a difficult task that will be investigated in the future with the CAST3M/ARCTURUS code. Coupling with neutronics is a key aspect of gas cooled reactors and CAST3M/CRONOS2 allows the study of these strongly coupled phenomena. Sensitivity studies have demonstrated the weight of the different parameters. Verification of this coupling will be the next step when some experimental results (HTTR) are available.

Fig. 1 .

 1 Fig. 1. GFR 600 MW concept: general view of the containment.

Fig. 2 .

 2 Fig. 2. GFR 600 MW-Proposed nodalization: front view.

Fig. 4 .

 4 Fig. 4. GFR 600 MW-Injected mass flow rate: CATHARE result of 1 in. break scenario.

Fig. 6 .

 6 Fig. 6. GFR 600 MW-ARCTURUS containment results of 1 in. break scenario: Break zone overpressure, maximum temperature and helium volumetric concentration.

Fig. 7 .

 7 Fig. 7. Comparison 1D-3D model and Equivalent Porous media modelvertical profiles (Top: helium temperature, Middle: homogenized fuel temperature, Bottom: fluid velocity.

Fig. 8 .

 8 Fig. 8. Verification of proposed model-3D temperature field: Solid Temperature between 480 and 1162 • C (1D-3D model).

Fig. 9 .

 9 Fig. 9. HTR-10-Vessel view.

Fig. 10 .

 10 Fig. 10. HTR-10-Specified mass flow-rate distribution.

Fig. 11 .

 11 Fig. 11. HTR-10-CAST3M/ARCTURUS results: radial temperature profiles (in the solid region).

Fig. 13 .

 13 Fig. 13. HTR-10-CAST3M/ARCTURUS: solid temperature and Velocity field (zoom).

Fig. 14 .

 14 Fig. 14. Schematic description of the thermal-hydraulic/neutronic coupling architecture.

Fig. 15 .

 15 Fig. 15. Reference meshing for the radial meshing sensitivity study.

Fig. 16 .

 16 Fig. 16. Different test configurations investigated in the sensitivity study.

Fig. 19 .

 19 Fig. 19. Geometry of the examined test configuration.

Fig. 20 .

 20 Fig. 20. Axial distribution of temperature differences between coarse and fine axial thermal-hydraulic grids.

Table 2

 2 GFR 600 MW -Summary of the thermal-hydraulics results for 1 in. break LOCA cumulated with station black-out (maximum overpressure and gaseous temperature, temperature and helium concentration levels in the injection zone after the equilibrium of pressure between containment and reactor vessel and steel liner maximum temperature increase)

	Case	P max (bar)	T max (K)	T rel.t>teq (K)	[He] rel.t>teq (vol%)	T max -WRH39 (K)
	1 in. break	1.7	393	≤ 375	≤ 13	4 (5000 s)

Table 3

 3 Main characteristics of the GFR 600 MW "CerCer performant"

	Thermal power (MW)	600
	Power density (MW/m 3 )	103
	Diameter of the helium channels (cm)	1.65
	Porosity (%)	55
	Height of the reactor core (m)	1.95
	Diameter of the reactor core (m)	1.95
	Injection helium mass flow rate (kg/s)	312.5
	Inlet helium temperature ( • C)	480
	where the non-dimensional numbers (Nu Nusselt number, Re
	Reynolds number and Pr Prandtl number) are:

Table 4

 4 Verification of the proposed model

	Characteristic	1D-3D model	2D axis model
	Fuel maximum temperature ( • C)	1161.9	1072.9
	Fluid maximum temperature ( • C)	848.8	849.6
	Nusselt number (min/max)	135-170	140-175
	Heat exchange coefficient (min/max)	3086-3276	3072-3249
	(W/m 2 /K)

Table 5

 5 HTR-10 benchmark: hot helium plenum temperature near the outlet (T sol : solid temperature, T flu : fluid temperature, dT/dr: radial solid temperature gradient and dT/dz axial solid temperature gradient)

	Nb.	R (cm)	Z (cm)	Calculated T sol ( • C)	Calculated T flu ( • C)	dT/dr ( • C/m)	dT/dz ( • C/m)	Experimental T ( • C)
	14	40	234	783	785	-159	-8	800.2
	15	60	234	740	742	-220	-55	763.1


H Lagrave

G Ranc

C Gallé

S Durand

Durability design of heated concrete structures Methodology and application to long-term interim storage

The operation of civil engineering structures subjected to thermal and mechanical loading has led the CEA to examine temperature-dependent variations in the concrete properties and the processes affecting the durability of these structures.

A new approach has been undertaken to specify the thermal, hydric and mechanical history of these structures. This technical approach is based on three areas of research: material characterization, modelling to identify weaknesses in the structure and validation by experimental tests on heavily instrumented structures subjected to representative loads.

The procedure adopted for long-term interim storage facilities [1, 2] can also be applied to other domains.

INTRODUCTION TO THERMAL HYGRO MECHANICAL (THM) PROGRAM

The main goals are studying the evolution with temperature of microstructural properties of concrete, understanding damage processes in heated concrete and during heat cycles and studying effect of temperature gradients and pore pressures on concrete strength.

The mechanical properties of concrete are determined by its mechanical, hydric, thermal and chemical history (degradation, transformation of hydrates, formation of new species) and loading time. The chemical transformations are accompanied by increasing porosity, loss of mechanical properties and increasing transfer properties with the temperature [START_REF] Bazant | Concrete at High Temperatures, Material properties and mathematical models[END_REF] . In the 80-200°C range which concerns us particularly, there is no clear threshold for the mechanical properties of concrete with a water-to-cement weight ratio below 0.5, with Portland cement.

Phenomenological models exist within the CEA to estimate the effect of drying on the mechanical behaviour and thermo-hydro-mechanical behaviour of the concrete.

These tools have been extensively validated by tests at centimetre scale but full-scale validation tests must still be carried out certain cases: this is the objective of the GALATEE program.

MATERIAL PARAMETERS CHARACTERIZATION

Experimental characterization focused on concrete thermo-hydric-mechanicals residual properties after thermal loading between 60 and 250°C (and up to 450°C in some cases). In it considered that for some concepts, high level waste canisters will be placed in concrete structures. In a normal evolution scenario, cooling systems will keep the temperature of these structures around 60°C. But in case of their temporary failure, high temperatures could be reached -up to 250°C - [START_REF] Ranc | Comportement à haute température du béton de structure. Impact de la fissuration sur les transferts hydriques[END_REF] . In such conditions, thermally related degradations, like drying, cracking or spalling, could drastically reduce concrete mechanical, confinement, and durability performances. It is essential therefore to characterize and to study the evolution of thermo-hydric-mechanical (T-H-M) properties of materials which could be used under such severe temperature conditions. All these results were part of a large data bank dedicated to concrete thermo-hydric-mechanical properties in severe thermal conditions. Concrete behaviour at high temperature depends on thermo-mechanical processes related to temperature gradients generating thermal expansion gradients [START_REF] Ulm | The channel fire. II: analysis of concrete damage[END_REF][START_REF] Noumowe | Effets de hautes températures (20 -600°C) sur le béton[END_REF] and on thermo-hydric processes related to water -liquid and vapour -movements [START_REF] Bazant | Concrete at High Temperatures, Material properties and mathematical models[END_REF][START_REF] Kalifa | Spalling and pore pressure in HPC at high temperatures[END_REF] . The kinetics and amplitude of such processes are controlled by micro-and macroscopic parameters. At material scale, those parameters are: thermal conductivity (), heat capacity, water and gas (k g ) permeability, porosity (), mechanicalcompressive (f c ) and tensile strength, elastic modulus (E) -and thermo-mechanical propertiesthermal expansion. Such properties were characterized for various high-strength concretes in order to constitute a data bank and to provide input data for the models. In addition to the potential impact of the intrinsic characteristics of the materials (water-cement ratio, cement and aggregates type), [START_REF] Gallé | Behaviour of high performance concrete under high temperature (60-450°C) for surface long-term storage: thermo-hydro-mechanical residual properties[END_REF][START_REF] Gallé | Permeability and Pore Structure Evolution of Silico-calcareous and Hematite High-strength Concretes Submitted to High Temperatures[END_REF] other influencing factors, designed as external characteristics, were investigated: long duration and cyclic thermal loadings, heating rate, concrete water saturation degree [START_REF] Gallé | Permeability and Pore Structure Evolution of Silico-calcareous and Hematite High-strength Concretes Submitted to High Temperatures[END_REF][START_REF] Gallé | Effect of the heating rate on residual thermo-hydromechanical properties of a high-strength concrete in the context of nuclear waste storage[END_REF] ... The majority of the phenomenological studies were carried out with two high-strength concretes (HSC). A standard HSC (w/c 0,43) was prepared with a CEM I-OPC cement (350 kg/m 3 ) and with silico-calcareous aggregates. A second one was manufactured with the same cement but with hematite aggregates (already used for nuclear radiobiological applications). Hematite is characterized by high density (, 5.10 3 kg/m 3 ), high thermal conductivity (11 W/m/°C) and high elastic modulus (209 GPa). Thermal loadings focused mainly on a temperature range between 60 and 250°C; extreme loadings were also performed up to 450°C in order to estimate margins). HS concretes characteristics are resumed in Erreur ! Source du renvoi introuvable.. Example of properties evolutions with temperature of such concretes under severe thermal loadings are provided within the figures 1 and 2. Concerning the effect of thermal treatment on microstructure, it was observed that the different heating generated important capillary porosity. This porosity is represented by a large population of pore access located in the region between 0.02 and 0.2 µm. This phenomenon was not only detected for high temperatures (over 110°C) but also for lower temperatures such as 60°C (Figure 1a). The development of cracking that influenced significantly the permeability was mainly observed over 250°C -the reference limit for accidental scenario -(Figure 1b), above all for the standard HSC. Concerning thermal conductivity, a linear evolution was globally observed for both concretes. Effect of thermal treatments on MIP pore access distribution (a) and effective gas permeability (b) for standard and heavy HSC [START_REF] Kalifa | Spalling and pore pressure in HPC at high temperatures[END_REF][START_REF] Gallé | Permeability and Pore Structure Evolution of Silico-calcareous and Hematite High-strength Concretes Submitted to High Temperatures[END_REF] About the evolution of mechanical properties with the temperature, it was detected that up to 110°C, HSC were characterized by an increasing compressive strength (heavy concrete) or at least a stable one (standard concrete). Over 110°C a monotonous decrease of this parameter was observed (Figure 2a). At 250°C a relative decrease of 16% was estimated for the standard concrete while for the heavy concrete the compressive strength loss was limited to 4%. The evolution of the elastic modulus (residual parameter) with the temperature was assumed to be linear. Similarly to the other parameters, once again it was observed that due to the mechanical characteristics of hematite aggregates, the residual elastic modulus for the heavy concrete (from 85 GPa at 60°C to 69 GPa at 250°C) was higher than the standard concrete ones (from 44 GPa at 60°C to 26 GPa at 250°C) for which the relative property loss -between 60 and 250°C -was about 50 %.

 (kg m -3 )  (%) k g (10 -19 m 2 )  (W m -1 K -1 ) f c (MPa) E (GPa) Standard 23788 10

a b

Figure 2. Evolution of compressive strength with temperature for standard and heavy HSC (a) and effect of heating rate on elastic modulus for standard HSC (b) [START_REF] Gallé | Behaviour of high performance concrete under high temperature (60-450°C) for surface long-term storage: thermo-hydro-mechanical residual properties[END_REF][START_REF] Gallé | Permeability and Pore Structure Evolution of Silico-calcareous and Hematite High-strength Concretes Submitted to High Temperatures[END_REF] Another specific study focused on the potential effect of the heating rate on the mechanical properties of a standard HSC. Numerical simulations showed that in case of the temporary failure of the cooling systems of the storage facility, high temperatures could be reached in the structures -up to 250°C. In that configuration, temperature rate increase should be around 0.1°C/min and not exceed 0.5°C/min. In the framework of this problem, impact of heating rates ranged from 0.1 to 10°C/min on concrete thermo-hydric-mechanical properties was estimated. The results showed that the lowest compressive strength loss was observed when using the lowest heating rate (0.1°C/min) and that this behaviour was consistent with previous studies performed on ordinary and HSC materials. On the contrary, an opposite trend was detected with the elastic modulus (Figure 2b) for which the lowest residual values were obtained with the 10°C/min rate. Nevertheless it was noticed that the impact of the heating rate was quite limited for both parameters (less than 10% relatively).

MODELLING OF CONCRETE BEHAVIOUR

Still in the framework of the same research program, several models have been developed to address different aspects of the behaviour of concrete at high temperature. The TH aspects are treated by a simplified TH model that rests on the assumption that the gaseous phase is only composed of vapour which permits to combine the two initial mass conservation equations of water in liquid and gaseous phase into only one. To deal with TM problems in concrete, a plasticity model has been developed. Finally, a coupled THM model based on the mechanics of partially saturated porous media [START_REF] Coussy | Mechanics of Porous Continua[END_REF] is being developed with the same assumption that the TH model and with a damage model to describe cracking.

The plasticity model used in the thermo-mechanical computations presented below is based on the model proposed by Feenstra and de Borst [START_REF] Feenstra | A composite plasticity model for concrete[END_REF] . The model has been generalized to 3D problems. The main reasons to choose this model were: a good description of the tensile and compressive behaviour, the possibility to take into account temperature dependent material parameters and input parameters allowing a simple identification. The stiffness degradation due to cyclic loading is not modelled. The model uses two yield criterions. A Drucker-Prager yield function and a Rankine yield function are used respectively in compression and in tension with the following non linear and exponential hardening laws ))

) ( 1 1   and ) ( 2 2   :                                 u u m m u m c m m m c if if f k if f k 1 1 1 1 1 2 1 1 1 1 0 1 1 2 1 1 1 0 1 1 0 ) ( ) ( 1 ) ( ) 1 ( 1 ) ( 2 2                          u t f 2 2 2 2 exp ) (    
1 exp( 1 ( 1 2    t f u hf G  , f G is the fracture energy, t
f the tensile strength and h an equivalent length which is a representative dimension of the mesh size. It is assumed that

A h   for a 2D domain and 3 V h  
for a 3D domain. Hardening laws are shown in figure 3. The input parameters of the model are Young's modulus, Poisson's coefficient, uniaxial and biaxial compressive strengths, yield limit, peak strain and ultimate strain in uniaxial compression, tensile strength and fracture energy.

In the example presented in the next paragraph, the parameters are temperature dependent. The evolutions of Young's modulus and compressive uniaxial strength are given by experimental data, as shown in chapter 2. The biaxial compressive strength is assumed to vary as the uniaxial compressive strength and the yield limit is assumed to vary as Young's modulus. Poisson's coefficient, thermal expansion, tensile strength and fracture energy are considered constant. a b

EXPERIMENTAL VALIDATION OF MODELS

All the tests carried out by the CEA share a high degree of attention to instrumentation, not only thermal and mechanical but also hydric, which is relatively unusual on large-scale structures. This is the case in particular for the MAQBETH test mockup, an annular reinforced concrete structure 2.2 m in diameter fitted with pore pressure sensors and relative humidity sensors; the unit can be heated to an internal temperature of 200°C for periods of several months. The instrumentation has been used to clearly identify water transport mechanisms in concrete in the presence of a strong temperature gradient and structural damage. In the hot compression-loaded zone, transfer in porous media corresponded to the mechanism determined using transport models based on capillary suction pressure and on the pressure of the gas mixture. In the cooler zone, however, water transfer was measured in liquid form via cracks under Darcy flow conditions. We were able to measure the amplitude of the fluid pressures, which did not exceed 10 bars. Hydric measurement data were compared with computations based on input data obtained by characterization of centimeter-scale test specimens (Figure 4a). Thermo-mechanical computations of the mockup using the plasticity model presented above have been carried out. A two-dimensional axisymmetric model of one half of the mock-up is considered. First a thermal analysis is carried out to determine the temperature distribution history within the structure. A thermo-mechanical analysis is then performed using these temperature fields as thermal loading. The temperature history applied by the heating device is imposed on the internal surface while a zero thermal flux is imposed at the ends because of the thermal insulation.

Temperatures obtained by a thermal analysis are in good agreement with measured temperatures. In the central area of the mock-up the stress state is close to a biaxial state: radial stresses are small compared to vertical and orthoradial stresses. The temperature gradient in the depth of the structure generates compressive stresses near the heated surface and tensile stresses near the external surface. In the compression zone, plastic deformations occur over a width of 8 cm whereas in the tension zone two-thirds of the depth of the structure experiences plastic deformation. Figure 4b shows the comparison of the vertical strains predicted by the thermo-mechanical model against measured strains given by two different gauges. The strain represented is the mechanical strain (total strain minus thermal strain) on the heated surface.

APPLICATION TO LONG-TERM INTERIM STORAGE

The third research topic of the French 1991 nuclear waste management act covers long-term interim storage. Two options are being investigated for high-level waste packages: surface (Figure 5) and subsurface (Figure 6) storage of two leaktight barrier canisters for spent nuclear fuel or vitrified wasteforms.

The guidelines determined the fundamental design options applicable to the concept and container:  no electromechanical components subject to wear in the storage facility,  safe containment of the radioactive inventory by two barriers constituted by the interim storage containers, the outer barrier being subjected to very slow and predictable corrosion at low relative humidity values,  continuous drainage by gravity of the storage facility,  dry storage and residual thermal power removal by natural air convection flow: transverse air flow in semi-underground vaults for surface storage [START_REF] Lagrave | Dossier de synthèse des études d'entrepôts de longue durée -version définitive -entreposage en surface des combustibles usés[END_REF] , axial air flow in vertical shafts for subsurface storage [START_REF] Lagrave | Dossier de synthèse des études d'entrepôts de longue durée -version définitive -entreposage en subsurface des combustibles usés[END_REF] .

Facilities design

The facility is designed to withstand the design-basis earthquake. The modules are capable of withstanding the impact of a fighter plane and ensure the protection of the containers from fire risks (burning kerosene)

The design takes into account the possible overloads arising in operation: thermal loads resulting from the heat released by the containers and mechanical stresses produced in the structure by the kinetics of concrete shrinkage.

The new approach described in this paper was applied to three strategic zones of the storage: the subsurface vault of the upper gallery, the subsurface shaft's structure and the surface slab. A design study of a shell for the storage shaft has been carried out based on the regulatory framework of Eurocode 2. In order to optimise the excavation dimensions and the durability of the structure, high performance concrete was selected with a mechanical strength exceeding 100 MPa, with high stress relaxation properties and stainless steel reinforcement to minimize the wall thickness Thermomechanical calculations were carried out for the shell to assess the impact of the rebars on the stress distribution (Figure 8) in the concrete, to determine the stress concentration zones versus height, and to examine the mechanical impact of different thermal and mechanical loading scenarios. The reinforcements were found to have a significant effect on the stress distribution, particularly the clips interconnecting the reinforcement layers. In addition, the sequencing was optimised to favour the interpretation of the results by beginning with mechanical loading without thermal loading with, then mechanical loading. This scenario limits damage to the shell and limits concrete cracking. Similar stress states can be obtained with a low heating rate.

DEMONSTRATION PROGRAMS

Galatee

The Galatee gallery is a demonstration tunnel for long-term subsurface storage. The experiments are designed to study the thermo mechanical response over time of the structure under normal operating conditions (corresponding to an interim storage facility loaded to full capacity).

The walls of the vault are heated continuously at 80°C and are instrumented to monitor the temperature and relative humidity of the concrete throughout its thickness in the course of drying. At the end of the test the residual deformation of the arch will compared with its initial geometry by local and topographic measurements.

A comparison wall (Figure 9) has been built to be able to allow destructive sampling of structural concrete representative of the construction of the vault. The industrial production of items such as the concrete shell based on the final design requires qualification of its long-term behaviour under high thermal and mechanical loads. A scale 1 test facility would be designed for this purpose.

The qualification test facility would be capable of reproducing the thermal and mechanical loads encountered in long-term interim storage, in normal and accidental situations.

CONCLUSION

A new approach has been undertaken to specify the thermal, hydric and mechanical history of these structures. This technical approach is based on three areas of research: material characterization, modelling to identify weaknesses in the structure and validation by experimental tests on heavily instrumented structures subjected to representative loads.

The procedure adopted for long-term interim storage facilities can also be applied to other domains: nuclear reactor containments, civil engineering works exposed to severe thermal and mechanical loading like new bridges or high deep tunnels for instance.

Figure 1 .

 1 Figure1. Effect of thermal treatments on MIP pore access distribution (a) and effective gas permeability (b) for standard and heavy HSC[START_REF] Kalifa | Spalling and pore pressure in HPC at high temperatures[END_REF][START_REF] Gallé | Permeability and Pore Structure Evolution of Silico-calcareous and Hematite High-strength Concretes Submitted to High Temperatures[END_REF]

Figure 3 . 1  , m 1  and u 1 

 3111 Figure 3. Hardening law in compression (a) and in tension (b)

Figure 4 .

 4 Figure 4. Comparison of predicted relative humidity in the layer of concrete versus measurement at different moments (a) and model predicted vertical strain versus measured strain on the heated surface (b)

Figure 5 . 6 .Figure 7 .

 567 Figure 5. Example of a surface storage vault Figure 6. Example of a subsurface storage shaft

Figure 8 .

 8 Figure 8. Reinforcement effect on  zz stress (in Pa)

Figure 9 : 6 . 2

 962 Figure 9 : Overall view during construction of the vault of GALATEE. View of the instrumentation embedded and of the comparison wall built during construction 6.2 Perspectives on shaft experimental program

Table 1 .

 1 Initial characteristics of HS concretes (at 20°C)

	HSC type


Chule Yang

Alessandro Renzaglia

Anshul Paigwar

Christian Laugier

Danwei Wang

Driving Behavior Assessment and Anomaly Detection for Intelligent Vehicles

Ensuring safety of both traffic participants and passengers is an important challenge for rapidly growing autonomous vehicle technology. To this purpose, intelligent vehicles not only have to drive safe but must be able to safeguard itself from other abnormally driving vehicles and avoid potential collisions. Anomaly detection is one of the essential abilities in behavior analysis, which can be used to infer the moving intention of other vehicles and provide evidence for collision risk assessment. In this paper, we propose a behavior analysis method based on Hidden Markov Model (HMM) to assess the driving behavior of vehicles on the road and detect anomalous moments. The algorithm uses the real-time velocity and position of the surrounding vehicles provided by the Conditional Monte Carlo Dense Occupancy Tracker (CMCDOT) framework. Next, by associating with the road information, the movement of each vehicle can be classified into several observation states, namely, Approaching, Braking, Lane Changing, and Lane Keeping. Finally, by chaining these observation states using a Markov model, the abnormality of driving behavior can be inferred into Normal, Attention, and Risk. We perform experiments using CARLA simulator environment to simulate abnormal driving behaviors, and we provide results showing the successful detection of abnormal situations.

I. INTRODUCTION

According to the Global status report on road safety 2018 [START_REF]Global status report on road safety 2018[END_REF], traffic accidents caused approximately 1.35 million deaths each year and up to 50 million injuries. With the rapid development of robotics, data science and sensor technology, many intelligent vehicles have appeared in our daily lives, such as Uber, Tesla, Waymo, and so forth. However, even the smartest autonomous vehicles cannot guarantee road safety and can cause lethal death or injuries [START_REF] Lee | The hype around driverless cars came crashing down in 2018[END_REF]. Among them, the inability to detect and identify risky driving behavior in time is an essential cause of accidents.

Although many works have been done to detect the abnormal scene by using surveillance videos [START_REF] Lin | Abandoned object detection via temporal consistency modeling and back-tracing verification for visual surveillance[END_REF], [START_REF] Lin | Visual-attention-based background modeling for detecting infrequently moving objects[END_REF] or autonomous robots [START_REF] Yang | Knowledgebased multimodal information fusion for role recognition and situation assessment by using mobile robot[END_REF], [START_REF] Yang | Probabilistic reasoning for unique role recognition based on the fusion of semanticinteraction and spatio-temporal features[END_REF], it is still challenging to detect the anomalous moments in highly dynamic traffic scenarios. Existing works on traffic scenario mainly involves vehicle trajectory prediction [START_REF] Glaser | Maneuver-based trajectory planning for highly autonomous vehicles on real road with traffic and driver interaction[END_REF], vehicle behavior prediction (lane-change) [START_REF] Galceran | Multipolicy decision-making for autonomous driving via changepoint-based behavior prediction: Theory and experiment[END_REF], and time to collision risk assessment [START_REF] Laugier | Dynamic scene understanding and upcoming collision prediction to improve autonomous driving safety: A bayesian approach[END_REF]. Few works focus on associating the anomalous driving maneuver with the behavior risk in a traffic scenario. The abnormal behavior of the vehicle can indicate the potential danger of the driver so that more attention should be paid to him. Thus, this paper proposes an approach to analyze the driving behavior of nonego vehicles and assess the behavior risk on the road. First, we use the Conditional Monte Carlo Dense Occupancy Tracker (CMCDOT) [START_REF] Rummelhard | Conditional monte carlo dense occupancy tracker[END_REF] to obtain the representative velocity and position of detected vehicles on the road. Next, by analyzing their real-time motion and position, observation states of each vehicle can be classified. The abnormality of driving behavior of each vehicle can then be inferred by chaining the observation states through Hidden Markov Model (HMM). The main contributions of this paper are as follows:

proposed algorithm. The rest of this paper is structured as follows. Section II reviews the related works. Section III introduces the basic idea of our perception framework, the CMCDOT. Section IV gives an overview of the proposed method and describe the methodology in detail. Section V shows the experimental procedures and results. Section VI concludes the paper and discusses future works.

II. RELATED WORK

A. Trajectory Prediction

Trajectory prediction is an essential function for autonomous vehicles to gain a more comprehensive understanding of the surrounding environments, to plan their own moving actions, and to avoid collisions. Existing methods either based on Recurrent Neural Network (RNN) or probabilistic modeling. RNN-based methods usually leverage a Long Short-Term Memory (LSTM) network to predict a smooth trajectory for autonomous vehicles. Quite a few works study highway scenarios by applying LSTM. [START_REF] Altché | An lstm network for highway trajectory prediction[END_REF] predicts future longitudinal and lateral trajectories for vehicles by training in a large dataset in various traffic densities rather than a low number of trajectories collected from a few drivers, while [START_REF] Deo | Multi-modal trajectory prediction of sur-rounding vehicles with maneuver based lstms[END_REF] predicts vehicle maneuvers by learning a model that assigns probabilities for different maneuver classes, and outputs maneuver specific predictions for each maneuver class. Other research using probabilistic modeling to deal with lane-change on the road [START_REF] Woo | Lane-change detection based on vehicletrajectory prediction[END_REF], and meanwhile considering the possibility of a crash with adjacent vehicles to reduce false alarms. [START_REF] Xie | Vehicle trajectory prediction by integrating physics-and maneuver-based approaches using interactive multiple models[END_REF] takes advantage of both physics-and maneuver-based methods and combines them to take into account both short-term and longterm effects for more accurate trajectory prediction.

B. Risk Assessment

The risk in this paper is the level of uncertainty that other vehicles introduce for ego-vehicles. For autonomous vehicles, the risk assessment needs to analyze the impact of the movement of surrounding vehicles on their normal maneuver, the likelihood of collisions, and the possibility of avoiding the danger. According to the cause of the risk, it is divided into collision risk based on future trajectory and risk based on unexpected behaviors [START_REF] Lefèvre | A survey on motion prediction and risk assessment for intelligent vehicles[END_REF]. The first type of risk assessment considers collisions to be the most critical risk and analyzes the risk based on the likelihood of collisions. The second risk assessment focuses on the dangerous situation caused by the driver's accidental operation. In the autonomous vehicle scene, the first case is more widely researched topic. [START_REF] Noh | Risk assessment for automatic lane change maneuvers on highways[END_REF] proposes two complementary risk indicators (collision time and minimum safety margin) in terms of time and space to cover all risk scenarios that may occur with vehicle lane changes on highways. This method reliably assesses the risk of collision. In [START_REF] Chen | Comparison of time to collision and enhanced time to collision at brake application during normal driving[END_REF], the authors present a large-scale data analysis to compare Time to Collision (TTC) and the Enhanced Time to Collision (ETTC) indicator with acceleration under a normal vehicle driving. The distribution of ETTC and TTC at the start of braking in normal following conditions is compared, and the probability model of TTC and ETTC values at braking at a series of vehicle speeds is given, which provides data for subsequent research in related fields. In [START_REF] Chen | A rear-end collision risk evaluation and control scheme using a bayesian network model[END_REF], the three indices of the transient speed of the following vehicle, the transient speed difference between the leading and following vehicles, and the following vehicle's acceleration (FA) are used to estimate the probability of collision risk between vehicles. In fact, when these quantities are relatively large, the vehicle does not have enough time to deal with the emergency. [START_REF] Zavareh | The effects of indicating rear-end collision risk via variable message signs on traffic behaviour[END_REF] considers TTC and safety margins, as well as speed as a risk perception indicator to assess the impact of information containing risk levels. The Stop Distance Index (SDI) is obtained, and the risk level classification is performed using the collision risk index based on the SDI value.

III. PERCEPTION FRAMEWORK: CMCDOT

The CMCDOT [START_REF] Rummelhard | Conditional monte carlo dense occupancy tracker[END_REF] is a perception system able to provide a dense and generic representation of the environment [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF], [START_REF] Moravec | Sensor fusion in certainty grids for mobile robots[END_REF] through probabilistic occupancy grids (see Fig. 2), based on Bayesian fusion, filtering of sensor data and Bayesian inference. In particular, by tracking spatial occupancy, it infers the dynamics of the scene via a hybrid representation of the environment, which is based on static and dynamic occupancy, free spaces and unknown regions (see Fig. 3). This differentiation enables the use of state-specific models, such as classic occupancy grids for static components and sets of moving particles for dynamic occupancy, as well as confidence estimation and management of areas with no information.

Additionally, an important feature of the CMCDOT, and main subject of study of this paper, is the velocity estimation for each cell of the grid. Most of the existing velocity estimation methods consist in detecting and tracking dynamic objects in the scene. Detection of moving objects requires to predefine shape of objects and involves fitting of object models in point clouds to find the location of objects in the environment. As the number and types of objects increases, the detection and tracking becomes a difficult and costly problem to solve. The grid-based approach used in the CMCDOT framework [START_REF] Rummelhard | Conditional monte carlo dense occupancy tracker[END_REF] instead of detecting objects, directly computes estimations of the position of every cell of the grid. By comparing these estimations with the ones obtained in previous time steps, the distinction between static and dynamic cells can be easily made, which in turn helps to estimate the velocity (speed and direction of motion) for each cell.

IV. METHODOLOGY

The overall workflow adopted for this work is shown in Fig. 4. The CMCDOT framework takes a point cloud and odometry data as input, generates an occupancy grid map and estimates velocities for each cell. A 3D object detection algorithm [START_REF] Paigwar | Attentional pointnet for 3d-object detection in point clouds[END_REF] is then used for detecting the position x of the vehicles and the bounding boxes. The representative velocity v for each vehicle is obtained by taking the mean value of the cells within the bounding box of the vehicle. By tracking and analyzing the real-time motion and position between the nonego vehicle and the ego-vehicle for consecutive time-steps, observation states can then be obtained. Finally, the driving behavior of each non-ego vehicle is inferred through a HMM.

In order to validate our methodology comprehensively, we rely on CARLA simulator to generate a large number of scenarios with abnormal driving behavior by parameterizing the movements of the vehicles. Simulation setup and scenario description are explained in detail in section V.

For each vehicle j detected in CMCDOT, we extract its overall velocity v j , acceleration a j and position x j . For the ego-vehicle, the perceptual parameters are denoted as v e , a e , and x e , respectively.

A. Motion and Position Analysis

To determine the movement state (observation state) of each non-ego vehicle, we first need to track and analyze their real-time motion and position and calculate the behavioral risk according to the distance, acceleration, and direction of velocity. After that, we normalize it to [0,1], and the resulting standardized behavioral risk means the degree of stability of driving behavior on the above three indicators. For each detected vehicle j, the corresponding risks are calculated as follows:

r a j = P (r a j |a j , a e ) = e aj -ae -e ae-aj e aj -ae + e ae-aj , r a

j ∈ [-1, 1] (1)
The acceleration risk r a j represents the relative acceleration between the vehicle j and the ego-vehicle. The greater the fluctuation in speed, the larger the absolute value of the risk.

r d j = P (r d j |v j , l) = cos( v j |v j | , l), r d j ∈ [0, 1] (2)
The direction risk r d j indicates the degree to which the vehicle j deviates from the expected road direction. l is a unit vector indicating the expected moving direction of the road. The larger the risk value, the smaller the deviation from the expected direction.

r p j = P (r p j |x j , x e , v j ) = cos( v j |v j | , x e -x j |x e -x j | ), r p j ∈ [-1, 1] (3)
The position risk r p j denotes the relative position between the vehicle j and the ego-vehicle. If the risk is greater than 0, it means that the vehicle j is in front of the ego-vehicle. Otherwise, the vehicle j is located behind the ego-vehicle.

B. Observation Classification

As illustrated in Fig. 5, the observation state O is defined as the movement status of vehicle that can be directly observed by the ego-vehicle. In this study, the possible movements are classified into four states, namely:

Approaching(o 0 ), Brake(o 1 ), Lane Changing(o 2 ), Lane Keeping(o 3 ), O ∈ [o 0 , o 1 , o 2 , o 3 ].
Let α, β, γ ∈ [0, 1] define thresholds for the respective acceleration risk, direction risk, and position risk.

1) Approaching: Approaching means that the position of the detected vehicle is behind the ego-vehicle, accelerating or having a higher speed than the ego-vehicle, and the direction of movement remains almost unchanged. 2) Brake: Brake means that the position of the detected vehicle is in front of the ego-vehicle, decelerating or having a lower speed than the ego-vehicle, and the direction of movement remains almost unchanged.

3) Lane Changing: Lane Changing means that the direction of the velocity of the detected vehicle is deviating from the expected direction (predetermined by the lane).

4) Lane Keeping: Lane Keeping means that the direction of the velocity the detected vehicle is consistent with the expected direction (predetermined by the lane).

Thus, the observation state can be obtained as follows:

O =       
o 0 , if r a j ≥ α and r d j > β and r p j ≥ γ o 1 , if r a j < -α and r d j > β and r p

j < γ o 2 , if r d j ≤ β o 3 , otherwise (4)

C. Driving Behavior Inference

The inferred abnormality of driving behavior B of the detected vehicle determines how much attention the egovehicle needs to allocate to it. In this study, the behavior can be inferred as N ormal(b 0 ), Attention(b 1 ), and

Risk(b 2 ), B ∈ [b 0 , b 1 , b 2 ]
. Thus, for each detected vehicle j, the inferred behavior at timestamp t can be obtained as follows:

P j (B t |O 1:t-1 ) = Bt-1∈B P j (B t |B t-1 ) Behavior Transition • P j (B t-1 |O 1:t-1 ) Previous Belief , (5)

V. SIMULATION & RESULTS

To validate our approach we rely on simulated scenarios generated with the CARLA urban driving simulator [START_REF] Dosovitskiy | CARLA: An open urban driving simulator[END_REF]. This is motivated by the fact that obtaining real data for abnormal driving can be unsafe. Furthermore, our detection pipeline is only based on two sensors: LiDAR and wheel odometer. LiDARs generate point clouds by measuring the time a beam of light takes to return after hitting an obstacle in the path. In comparison to other sensors (e.g. cameras) simulated LiDAR sensor data are very similar to real data [START_REF] Petit | Remote attacks on automated vehicles sensors: Experiments on camera and lidar[END_REF].

CARLA simulation environment consists of complex urban layouts, buildings, and vehicles rendered in high quality, allowing for a realistic representation of real-world scenarios. The ego-vehicle and its sensors, as well as other moving vehicles, as depicted in Fig. 3, are configured in the simulation to match with the actual system. We developed a parameterbased approach which streamlines the process through which the velocities, position and steering angle of non-ego vehicles are specified.

A. Scenario Description

Our simulation scenario focuses on the behavior of cars in a multi-lane highway. For simplification, we bound our simulation to a scenario with one ego vehicle and two nonego vehicles. Also, we fix the lanes and behavior for each simulated vehicles, as shown in Fig. 6. Non-ego vehicle-1 is in left lane moving forward with randomized acceleration and deceleration parameters over time. Non-ego vehicle-2 is in the right lane also moving forward but with randomized minor alterations in steering angle. The ego vehicle is in the center lane moving forward with a constant velocity.

B. Results

The proposed method has been tested on the simulated scenario presented in Fig. 6 and the assessment results of the non-ego vehicle-2 are shown in Fig. 7 The inference result of non-ego vehicle-2 is shown in Fig. 7. In the beginning, the non-ego vehicle-2 is located in front of the ego-vehicle but is caught and overtaken by the egovehicle within the next 20 seconds. The second chart describes the real-time velocity (both longitudinal and lateral) of the ego-vehicle and the non-ego vehicle-2. The ego-vehicle keeps a relatively steady speed in both longitudinal and lateral direction, while the non-ego vehicle-2 has a fluctuation of velocity either on magnitude or direction. The speed variation on the longitudinal direction will cause the observed state to be Approaching or Brake according to their relative position, and the speed fluctuation on the lateral direction will cause the observed steady to be Lane Changing. As a result, the behavior can be inferred as Attention when the motion of non-ego vehicle-2 is fluctuating. But if the fluctuation is beyond a certain level, it will be judged as Risk. As you can see the timestamp at around 780 and 786, the sudden speed change of the non-ego vehicle-2 leads to the observed states to be Lane Changing and Brake, which leads to the inference to be Risk.

Similar to non-ego vehicle-2, the inference result of non-ego vehicle-1 is shown in Fig. 8. The non-ego vehicle-1 keeps a similar motion and relatively stable distance to the ego-vehicle, but have some sudden changes in the direction and magnitude of velocity. The proposed algorithm successfully detect the sudden Brake and sudden Lane Changing, and inferred to Risk after a period of Attention.

VI. CONCLUSIONS

In this paper, we proposed a systematic approach for vehicle driving behavior assessment and detection of anomalous moments. The proposed method takes advantage of CMC-DOT framework to obtain real-time motion and velocity of vehicles. Then, the ego-vehicle can observe the movements of nearby vehicles and infer their driving behavior by tracking and analyzing their motion and position. We used CARLA simulator to generate several typical abnormal movements for two vehicles. The results show that the proposed method can successfully detect the risk moments. For future work, we intent to use machine learning methods, such as learning from demonstrations, to improve behavioral reasoning. We will also add more road users, such as pedestrians or buses, to study the impact of different semantic objects on decision making.

Fig. 1 :

 1 Fig. 1: (Top) Simulation environment with CARLA simulator. The white vehicle is the ego-vehicle and two non-ego vehicles are simulated to perform anomaly movements. (Bottom) Perceptual environment with CMCDOT framework. By analyzing the real-time velocity and position of the vehicle, the state and behavior of vehicles can be inferred.

Fig. 2 :

 2 Fig. 2: Data representations in CMCDOT formulation. The environment is divided into cells, to which static, dynamic, empty and unknown coefficients are associated. The dynamic part is expressed by weighted particles which sample the velocity space.

Fig. 3 :

 3 Fig. 3: (top) The environment around the vehicle as observed using a 64 layered LiDAR. (Bottom) Representation of the environment as dense probabilistic occupancy grid output from CMCDOT framework. Colors represents different states of each cell in the grid: occupied static (blue), occupied dynamic (green), empty (black), unknown area (red).

Fig. 4 :

 4 Fig. 4: Abnormal driving behavior detection pipeline.

Fig. 5 :

 5 Fig. 5: Graphical model of the proposed algorithm. v and x are the velocity and position of each vehicle, O is the observation state, and B is the inferred vehicle driving behavior.

  and Fig. 8. In these figures, the first two charts show the real-time relative position between the ego-vehicle and the detected vehicle and their velocities. The third chart shows the observation sequence of the detected vehicle, where o 0 , o 1 , o 2 , o 3 represent Approaching, Brake, Lane Changing, and Lane Keeping. The last chart shows the inferred driving behavior of the detected vehicle, where b 0 , b 1 , b 2 represent Normal, Attention, and Risk.

Fig. 6 :

 6 Fig. 6: Simulation environment and scenario in this study. In the upper part of this figure, the white car is the ego-vehicle, and the red car on the right is the detected vehicle with some risk behaviors. The lower part of his figure shows the perceptual results by CMCDOT.

Fig. 7 :

 7 Fig. 7: Experimental results of the non-ego vehicle-2. Sudden lane changes toward the ego-vehicle are detected and anomalous moments are inferred.

Fig. 8 :

 8 Fig. 8: Experimental results of the non-ego vehicle-1. Sudden brake and lane change toward the ego-vehicle are detected and anomalous moments are inferred.

School of Electrical and Electronic Engineering, Nanyang Technological University, 639798 Singapore; e-mail: yang0438@e.ntu.edu.sg, edwwang@ntu.edu.sg

2Univ. Grenoble Alpes, Inria, Chroma, F-38000 Grenoble, France; e-mail: alessandro.renzaglia@inria.fr, anshul.paigwar@inria.fr, christian.laugier@inria.fr

• Four observation states are defined and classified by analyzing the motion and position information obtained from CMCDOT. • Driving behaviors are linked to observation states based on an HMM, to infer anomalous behaviors. • We demonstrate the use of CARLA simulator[START_REF] Dosovitskiy | CARLA: An open urban driving simulator[END_REF] to generate abnormal driving behaviors and validate the

ACKNOWLEDGMENT

This work was supported in part by the ST Engineering-NTU Corporate Laboratory through the NRF Corporate Laboratory@University Scheme. This work has also been conducted within the ENABLE-S3 project that has received funding from the ECSEL joint undertaking under grant agreement NO 692455. This joint undertaking receives support from the European Union's Horizon 2020 research and innovation programme and Austria, Denmark, Germany, Finland, Czech Republic, Italy, Spain, Portugal, Poland, Ireland, Belgium, France, Netherlands, United Kingdom, Slovakia, Norway.