Michel Belliard
email: michel.belliard@cea.fr

NUMERICAL MODELING OF AN IN-VESSEL FLOW LIMITER USING AN IMMERSED BOUNDARY APPROACH

Keywords: PWR T-H system, In-vessel flow limiter, CFD, Immersed Boundary Methods

whether they are published or not. The documents may come from teaching and research institutions in France or

INTRODUCTION

The context of this work is set in the domain of Generation III nuclear power plants. There are the innovative designs that are under construction or still in design phase [START_REF]Advanced Reactors Information System (ARIS)[END_REF]. More specifically, we focused on the light-water Pressurized Water Reactors (PWRs), which are the main type of reactors built and exploited in France. Nowadays passive safety systems are more and more included in the nuclear-reactor safety strategy to mitigate design basis accidents. A passive safety system is a system that activates itself without the need of mechanical or electrical actuation. At CEA, some studies on passive safety systems have been done in the past years, notably for the in-vessel flow limiter (hydraulic diode) patented by the CEA [START_REF] Gautier | Dispositif limiteur de débit inverse de fluide[END_REF] designed to limit the amount of water lost during the short-term sequence of a Large-Break (LB) Loss Of Coolant Accident (LOCA), cf. Fig. 1. The goal is to interfere as little as possible with the nominal operation flow and maximize the fluid retained in the event of this kind of accident.

Numerical investigations of the benefices induced by the in-vessel flow limiter during LOCA transients have been done using thermal-hydraulic system safety codes. But, the relevance of these system-scale studies depends on the level of realism of the data introduced to take into account the hydraulic diodes. In order to take into account the large-scale effect of hydraulic diodes in safety-system studies, we need information coming from small-scale experimental or numerical experiments. For instance, results from 3D Computational Fluid Dynamic (CFD) simulations can be used to up-scale relevant characteristics as the global pressure drop induced by the hydraulic diode. Moreover, CFD studies can be involved in the optimization process of the flow-limiter geometry to minimize the global pressure drop during nominal operations and to maximize it during a cold-leg LB-LOCA. This optimization process usually needed a big number of simulations. But, as many geometric scales are presented at the same time (typically several meters for the downcomer radial scale and one centimeter for the fin thickness), a CFD simulation of the two-phase flow inside the flow-limiter can be time consuming. Instead of this conventional approach (i.e. body-fitted approach), we are motivated by fictitious domain approach [START_REF] Saul'ev | On the solution of some boundary value problems on high performance computers by fictitious domain method (in Russian)[END_REF][START_REF] Rukhovets | A remark on the method of fictive domains[END_REF] allowing a less precise but faster estimation of the pressure-drop coefficient. Following this way, we consider simulations over a full computational domain including the in-vessel flow-limiter fins and re-introduce their presence adding local external forces on the immersed boundaries. Moreover, a homogeneous relaxed equilibrium model of a liquid-vapour mixture [START_REF] Clerc | Numerical simulation of the Homogeneous Equilibrium Model for two-phase flows[END_REF] can be considered. For same space discretization, this kind of three-balance-equation model generally run faster than a six-balance-equation model. Once a particular geometry exhibited as a good candidate, a reduced number of body-fitted CFD computations can be done to refine the design.

In this paper, we present a methodology concerning the design of a fast-running twophase CFD model of the in-vessel flow-limiter device illustrated by CFD results. The numerical/experimental validation of this simulation tool is not the goal of this paper and only brief elements of validation are given here. The paper is structured as follow. The two-phase fluid CFD model and the first-order in space immersed boundary models are first presented in Sections 2 and 3. Then, the in-vessel flow limiter study, with the GENEPI code as CFD tool, is discussed in Section 4. Results provide a range of values that may be used in input of safety-system codes. Finally, some words concerning the perspectives of this work are given in Section 5.

The two-phase CFD model

The considered two-phase CFD model is the GENEPI one [START_REF] Grandotto | Calculs des écoulements diphasiques dans les échangeurs par une méthode aux éléments finis (in French)[END_REF][START_REF] Grandotto | Steam generator two-phase-flow numerical simulation with liquid and gas momentum equations[END_REF], designed for the steamgenerator two-phase flow steady-state 3D computations through the resolution of three balance equations for a water liquid/steam mixture. It is based on a homogeneous relaxed equilibrium model with thermodynamic equilibrium of the two phases. But closure laws take into account the liquid/steam momentum disequilibrium. This code incorporates the possibility to model thin no-penetration obstacles using Immersed Boundary Conditions (IBCs) [START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF]. Provided that the following assumptions hold: (i) surface tension, viscous and turbulent dissipation are neglected and pressure terms are neglected in the enthalpy balance equation, (ii) steam and liquid have same pressure, (iii) an eddy viscosity model is considered, one obtains for the mixture description of the two-phase flows the following mass, momentum and enthalpy balance equations:

∇.G = 0, ( 1
)
ρ∂ t V + G. ∇V + ∇.(x(1 -x)ρV R ⊗ V R ) = ρg -∇P -ΛV + ∇.µ T ( ∇V + ∇T V), (2)
ρ∂ t H + G.∇H + ∇.(x(1 -x)ρLV R ) = ∇.(χ T ∇H), (3)
with Λ the singular-obstacle tensor and V R the relative velocity given by the drift-flux Lellouche-Zolotar model [START_REF] Lellouche | Mechanistic model for predicting two-phase void fraction for water in vertical tubes, channels and rod bundles[END_REF] and based on the Zuber-Findlay approach [START_REF] Zuber | Average volumetric concentration in two-phase flow systems[END_REF]. The turbulent dynamic viscosity µ T is given by the local scalar Schlichting model [START_REF] Schlichting | Boundary layer theory[END_REF]:

µ T = a S |G|L T (4)
where L T is a characteristic length and a S a coefficient. The turbulent diffusion coefficient χ T for the enthalpy balance equation is defined via the Prandtl number P r = µ T χ T . The density ρ, the static quality x and the latent heat L are determined through the equation of state of the water as a function of the pressure and mixture specific enthalpy. We solve in H, P and V variables. The time term presence in Eq. ( 2) allows to search the steadystate regime through a transient computation of a thermally dilatable fluid, ∇.G = 0, cf.

Eq. ( 1). For that, the Chorin-Gresho method [START_REF] Gresho | On the theory of semi implicit projection methods for viscous incompressible flow and its implementation via finite element method that also introduces a nearly consistent matrix. i, theory[END_REF] (a fractional-step method) is used to solve the coupled mass-momentum equations. The non linearity are dealt by the Picard iterative process. The time discretization is based on a semi-implicit Crank-Nicholson scheme. The spatial discretization is based on the unstructured hexahedral finite elements (constant pressure by element and tri-linear velocity by node). The physical data ρ, µ T and Λ are constant by element. The Streamline Upwind Petrov-Galerkin method is applied to correct the convective term [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]. A conjugated gradient method, preconditioned by the diagonal, is used to solve the arising linear systems.

According to the hyperbolic nature of the flow equations, Dirichlet boundary conditions are used at the inlets of the domain (mass flux and enthalpy) and Neumann boundary conditions at the outlets (pressure). The other boundaries of the domain are impermeable walls. Generally, these are considered adiabatic and with no shear stress.

The Immersed Boundary Model

In the fictitious domain approach, introduced in the fifties by Hyman [START_REF] Hyman | Non-iterative numerical solution of boundary-value problems[END_REF] and the Russian's school [START_REF] Saul'ev | On the solution of some boundary value problems on high performance computers by fictitious domain method (in Russian)[END_REF][START_REF] Rukhovets | A remark on the method of fictive domains[END_REF], the original domain Ω is embedded in a fictitious domain Ω which is geometrically bigger and generally simpler-shaped. Doing this, some immersed boundary Σ appears such that Ω = Ω ∪ Σ ∪ Ω e , where Ω e is the complementary or 'exterior' domain (as the fins). The spatial discretization is now performed in Ω, independently of the shape of the original domain Ω. Then, the resolution of the new problem in Ω will be faster and simpler. The main issue is to enforce the original boundary conditions on the immersed interface Σ which is non-aligned with the mesh.

The ISI method

In this work, among the numerous fictitious domain methods (see [START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF] for a short introduction), we consider an element of the set of Immersed Boundary Methods: the fictitious domain method with Immersed Spread Interface (ISI) [START_REF] Ramière | A fictitious domain approach with spread interface for elliptic problems with general boundary conditions[END_REF][START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF]. The fictitious problem to be solved in Ω is built from the original problem in Ω, but an additional term takes into account the immersed boundary conditions. For velocity Dirichlet boundary conditions, the singular-obstacle tensor Λ of Eq. (2) will play this role. It allows us to take into account the no-penetration condition of the flow limiter fins in an implicit way during the first step of the Chorin-Gresho method. Let us (u, v, w) be the local basis linked to a given obstacle (i.e. fin). The vectors u and v are tangential to the obstacle and the vector w is normal to the obstacle. For the element e, we define the singular-obstacle tensor by:

Λe = ρ e A e Ω e   Λ u 0 0 0 Λ v 0 0 0 Λ w   (5)
with A e the measure of this obstacle (area, m 2 ) intercepted by the element e and Ω e the measure (volume; m 3 ) of this element. Λ u , Λ v and Λ w are the tensor coefficients in, respectively, the directions u, v and w. Here, we consider no-penetration obstacles in the normal direction, Λ w = 1/ with 0 < << 1, and slip conditions in the tangential directions, Λ u = Λ v = 0. Immersed interfaces Σ, as the flow-limiter fins, are modeled by a collection of linear plane surfaces intercepting elements of Ω. In each intercepted element, the measure of the intercepted surface and the external normal vector are known. As the singular-obstacle tensor is defined by element, all the nodes belonging to this element are concerned: i.e. we have a spread interface approximation of Σ.

Convergence order and elements of validation

As for the L 2 -penalty methods [START_REF] Saul'ev | On the solution of some boundary value problems on high performance computers by fictitious domain method (in Russian)[END_REF][START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], that the ISI method generalizes, the parameter is known as the penalty parameter. Regardless to the Navier-Stokes solving method, the theoretical rate of convergence of the penalized solution toward the body-fitted one is comprised in the range [O( 1/4 ); O( 1 )] in L 2 (Ω) norm [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. Let us notice that for Dirichlet boundary conditions and elliptic problems, the theoretical rate of convergence in space of the Q 1 -finite element method with non-boundary-fitted meshes is O(h 1 ) in L 2 (Ω) norm, with h the space step [START_REF] Ramière | Convergence analysis of the Q1-finite element method for elliptic problems with non-boundary-fitted meshes[END_REF]. As a whole, contributions to the validation of this IB approach can be found in [START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF] and [START_REF] Introïni | A second order Penalized Direct Forcing for hybrid Cartesian/Immersed Boundary flow simulations[END_REF]. On one side, in the context of dilatable two-phase flow elliptic problems, the work mentioned in [START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF] validates the ISI method with respect to body-fitted finite-element computations and to the JEBC method (an IB method using a finite-volume discretization). A first-order rate of convergence in space is numerically reached. On the other side, in the context of incompressible one-phase flow Navier-Stokes equation, the work mentioned in [START_REF] Introïni | A second order Penalized Direct Forcing for hybrid Cartesian/Immersed Boundary flow simulations[END_REF] gives elements of validation for a finite-volume first-order penalty method very similar to the ISI method. Again, a first-order rate of convergence in space is numerically reached on the test case of a laminar flow around a static cylinder of diameter D (Reynolds number = 20). These results are in very good agreement with those proposed in the literature, cf. Table 1. In a lesser degree, it is also true for our own GENEPI results using the ISI method (about 10% on the drag coefficient and 30% on the recirculation length) giving confidence in the ability to catch the magnitude of an obstacle's drag coefficient.

Table 1: Hydrodynamic coefficients associated with the problem of steady flow around a static cylinder of diameter D (Reynolds number 20). B.F.: Body-fitted. C d : drag coefficient. L w : recirculation length. [START_REF] Introïni | A second order Penalized Direct Forcing for hybrid Cartesian/Immersed Boundary flow simulations[END_REF] (base) refers to the first-order penalty method of [START_REF] Introïni | A second order Penalized Direct Forcing for hybrid Cartesian/Immersed Boundary flow simulations[END_REF]. ISI: 18 cells in the diameter D and = 10 -5 .

GENEPI

References

CFD studies of the in-vessel flow limiter

This section is devoted to a first insight into the hydraulic of the in-vessel flow limiter through CFD simulations using the GENEPI code. We have set-up a preliminary design of the in-vessel flow limiter, cf. Fig. 3. In this preliminary work, we mainly restrict ourselves to liquid one-phase flows. But it is not a limitation; cf. [START_REF] Belliard | Numerical modeling of an in-vessel flow limiter using an Immersed Boundary Approach[END_REF] for illustrations of simulations with two-phase flows. For simplicity reason, we do not give attention to the gravitational acceleration term in the computations (around the flow limiter, gravitational effects are negligible in comparison with the inertial ones during the fast-depressurization phase). As well, the pressure range considered here is only [40-60] bar. Concerning the turbulence model, cf. Eq (4), the GENEPI-code standard value for the Schlichting coefficient is a S = 0.047 and the turbulence characteristic length L T is related to the biggest eddy structures. As reference, we choose L T ≈ 1 m for the nominal-operation flow direction (azimuth scale in the down-comer) and L T ≈ 0.3 m for reverse flow direction toward the broken cold leg (≈ radial scale in the down-comer). As the GENEPI's turbulence model is quite rough, a parametric study is performed on the coefficient a S and the turbulence characteristic length L T , cf. Eq (4): a S ∈ {a=0.047, a/10} and L T ∈ {0.3, 1.0, 2.0}. The value of the penalty parameter is set to = 10 -5 . We consider that the GENEPI steady state is reached when the relative L 2 -norm difference of the variables (pressure, mass flux and enthalpy) between two consecutive time iterations is less than 5.10 -3 δt with δt the time step. Usually a CPU time of 5 to 6 hours is need to reach the steady state on the finest mesh M 3 . It is compatible with a big number of simulations needed to optimize the flow limiter geometry.

Computational domain and meshing

The CFD computational domain Ω is a simplified rectangular geometry of 4.3 m x 4 m x 0.2 m to which was added the broken cold-leg nozzle starting about 1 m before the down-comer. It extends up to 2 m below the cold-leg axis. Considering a sectorized down-comer as in Fig. 1, we only model one third of the down-comer, including one coldleg entry in the vessel and one hot-leg pipe, cf. Fig. 2. The measures of the computational domain volume, of the down-comer and the cold-leg nozzle surfaces are equal to 3.8 m 3 , 0.9 m 2 and 0.4 m 2 respectively. For the purpose of a mesh convergence study, three meshes M 1 , M 2 and M 3 were built involving N 1 = 6,080, N 2 = 48,640 and N 3 = 164,160 elements respectively. The mean space step ranges from 0.2 m (M 1 ) to 0.07 m (M 3 ) and the ratio between two consecutive space-step indexes is 2.0 (M 1 → M 2 ) and 1.5 (M 2 → M 3 ). The flow-limiter fins are modeled by singular-obstacle surfaces through a collection of plane surfaces, cf. Fig. 3.

Boundary conditions

Slip-wall boundary conditions are considered on the walls, mass flux Q in is imposed on the in-flow boundary and ad-hoc pressure P out = 50 bar on the out-flow boundary. The in-flow and out-flow surfaces are the down-comer and the cold-leg nozzle surfaces depending on the considered main-flow direction. The inlet mass-flux values are Q in = 4,690 kg/s in nominal-operation condition (default flow direction; in ↔ cold leg and out ↔ downcomer) and Q in = 5,200 kg/s in LB-LOCA condition (reverse flow direction; in ↔ downcomer and out ↔ cold leg).

Results

All quantitative results concerning the pressure-drop coefficients for the various turbulencemodel coefficients and meshes are grouped in Tab. 2. In this table, the bold-typed turbulence parameters are the reference ones. Also the global pressure-drop coefficients computed by GENEPI, with or without in-vessel flow limiter, on the finest mesh M 3 are bold-typed. These can be compared to the global pressure-drop coefficients found in literature (Borda-Carnot law or Idel'cik). Although the mesh convergence is not fully reached, the trend of the evolution of the pressure-drop coefficient versus the space-step index is globally caught. Figs 4 to 6 present some field distributions concerning the mixture velocity, the mixture Table 2: Summary of the global pressure-drop coefficients K gl . The default flow direction is defined as the nominal-operation flow direction (from the cold leg to the down-comer). The reverse flow direction is defined as the opposite direction. (*): unsteady computation.

Without limiter

With limiter Flow direction Turbulence [START_REF] Idel'cik | Mémento des pertes de charge[END_REF] K gl ≈ 1. Finally, considering Fig. 6, we conjuncture that the impact on the flow of the fins located at the bottom of the flow limiter is much greater during LB-LOCA condition than during nominal-operation one.

K M 1 K M 2 K M 3 K M 1 K M 2 K M 3 Default a; L T =2 m 2.

Conclusions and perspectives

In this paper, we have presented a CFD model of an in-vessel flow limiter to mitigate the consequences of a large-break loss of coolant accident in a pressurized-water reactor. The principle of this safety device is based on fins designed to create a strong flow vortex increasing the pressure drop toward the broken cold leg. A rough CFD model using a homogeneous relaxed equilibrium model of a liquid-vapour mixture and an immersed boundary approach has been set-up, allowing a less precise but fast estimation of the pressure drop following the geometry of the fins. Through 3D/1D up-scaling of a global pressure-drop coefficient, local pressure-drop coefficients can be provided to thermal-hydraulic system safety codes, allowing the study of the in-vessel flow limiter effect on the thermal-hydraulic system. One-phase fluid CFD simulations have been run using the GENEPI code with a computation domain defined in coherence with the downcomer of a PWR. Parametric studies on the turbulence model lead to determine the range of the global pressure-drop coefficients depending on the direction of the flow and the presence of the flow limiter, opening the way to take into account the in-vessel flow limiter in a realistic manner. Nevertheless, considering the limitations related to the turbulence model, this conclusion needs to be consolidated by body-fitted CFD studies with more precise turbulence models. Some perspectives can be outlined about the improvement of this numerical model in order to proceed to the geometry optimization of the design of the flow-limiter fins. In particular, we can mentioned the space-interpolation scheme across the boundary interface to reach the second order [START_REF] Introïni | A second order Penalized Direct Forcing for hybrid Cartesian/Immersed Boundary flow simulations[END_REF] and the definition of immersed-wall laws for RANS/large-eddy simulations.

Figure 1 :

 1 Figure 1: Scheme of in-vessel flow limiters (hydraulic diode) located between the cold legs and the downcomer [2].

  (a) Down-comer surface (green).(b) Cold-leg nozzle surface (green).

Figure 2 :

 2 Figure 2: Example of the mesh used for the CFD study of the in-vessel flow limiter (mesh M 2 ; 48640 elements). The walls are colored in blue.

Figure 3 :

 3 Figure 3: Example of preliminary meshing of the flow limiter.

2

 2 Borda-Carnot (ξ = 1) K gl ≈ 1.1 pressure and the local external forces taking into account the flow-limiter fins ΛV, cf. Eq. (2) for the two studied flow configurations (nominal condition and reverse condition). For the nominal-operation configuration (fluid flowing toward the down-comer), the pressure and velocity distributions are shown in Figs 4(a) and 5(a). The local external-force distribution is illustrated in Fig.6(a). The introduction of the flow-limiter device induces an irregular flow path with the emergence of flow channels between the fins. Also, the in/out-flow pressure difference is increased in reference to the case without in-vessel flow limiter:K M 3 ∈ [-0.3; +0.1] → [0.7; 1.6], cf.Table 2.Let us notice that, without flow limiter, the range of CF-computed values K M 3 ∈ [-0.3; +0.1] is compatible with the Idel'cik (output of a rectilinear diffuser on a screen) or Borda-Carnot (sudden flow-section expansion) estimation one K gl ∈ [-0.5; +0.5]. For the reverse-flow configuration (fluid flowing toward the broken cold leg), the pressure and velocity distributions are shown in Figs 4(b) and 5(b). The local external-force distribution is illustrated in Fig. 6(b). The introduction of the flow-limiter device clearly increases the fluid vortex at the entry of the broken cold leg and the in/out-flow pressure drop. The global pressure-drop coefficient is multiplied by almost a factor two in case of flow limiter: K M 3 ∈ [3.6; 3.7] → [5.4; 6.4], cf. Table 2. Obviously the geometry of the fins Michel Belliard (a) Nominal-operation condition (b) LB-LOCA condition

Figure 4 :

 4 Figure 4: Pressure; mesh M 3 (164,160 elements); (a S = a; L T = 0.3 m). Pressure iso-values ranges from 4.8 to 5.2 bar.

Figure 5 :

 5 Figure 5: Velocity; mesh M 3 (164,160 elements); (a S = a; L T = 0.3 m).

  has to be optimized to enhance this effect while limiting the flow-limiter impact during nominal operations. Contrary to the nominal-operation case, the range of CF-computed values without flow limiter K M 3 ∈ [3.6; 3.7] sensibly differs from the Idel'cik (conical collector with front wall and screen) or Borda-Carnot (sudden flow-section reduction) estimation one K gl ∈ [1.1; 1.2]. But the hydraulic path is quite complex and does not reduce to simple configurations.

  (a) Nominal-operation flow; (a S = a; L T = 0.3 m). (b) Reverse flow; (a S = a; L T = 1.0 m).

Figure 6 :

 6 Figure 6: Local external-force distribution; mesh M 3 (164,160 elements).