M Belliard
email: michel.belliard@cea.fr

A Gerschenfeld
email: antoine.gerschenfeld@cea.fr

S Li

A LOCAL DEFECT CORRECTION SCHEME FOR THE COUPLING OF THERMAL-HYDRAULICS SAFETY-SYSTEM CODES WITH COMPUTATIONAL FLUID DYNAMICS CODES

Keywords: Sodium fast reactor, Multi-scale coupling, Thermal-hydraulics safety-system code, CFD, Local Defect Correction method

In the framework of the 4th-generation sodium fast reactor ASTRID, numerical simulations are widely used to assess the design and the safety of the reactor. In particular, the need of tools with an increased precision is growing, including multi-scale coupling. For instance, the pipes and the exchangers of the primary or secondary circuit are well described by thermal-hydraulics safety-system codes, but large volumes such as collectors need to be described by Computational Fluid Dynamics (CFD) codes. It is of prime importance to have robust and powerful coupling tools for increasingly complex studies in order to get the most benefit of these studies in term of reactor performance and of safety margin. This paper concerns the resolution methodology of the multi-scale coupled problems appearing when simultaneously considering the systemscale thermal-hydraulics description and the local-scale CFD description. We present the implementation of the Local Defect Correction method (LDC) -a coupling method for an enhanced local description by overlap of a local domain with a global one -into the CEA CATHARE3 system code coupled with the CEA TrioCFD, TrioMC (core thermal hydraulics) codes via the specifically developed MATHYS coupling interface. After the presentation of the LDC method and its application to the CATHARE's balance equations, we illustrate the validity of this approach by the presentation of results concerning di-phasic problems obtained by the coupling CATHARE3/CATHARE3. Finally, we list the perspectives emerging for the application of the methodology to monophasic problems concerning the simulation of ASTRID reactor through the CATHARE3/TrioCFD coupling.

introduction

In the context of the future ASTRID 4th-generation reactor, the CEA aims to carry out detailed thermohydraulic simulations of the vessel. For this goal, system-scale and local-scale coupling is a necessary pathway. It is of prime importance for the CEA to have robust and efficient tools for ever more complex studies and to master them in their globalities (physics, numeric, computing) in order to benefit from them. This work is done in this context. It presents the potential application of the Local Defect Correction (LDC) method to the CATHARE system thermodynamics code serving the CATHARE/TrioCFD (or TrioMC) coupling.

The objective is to increase the robustness of the coupled calculations.

The French reference thermal-hydraulics system safety code CATHARE [START_REF] Barré | The CATHARE code strategy and assessment[END_REF], originally devoted to the study of water-cooled reactor transients (standard operations or accidental transients from any kind of failures or size and location of breaks), is based on 0D/1D and 3D modules using six-equation (mass, momentum and energy) two-fluid models on a staggered space discretization and solved implicitly in time. The time step is dynamically computed. The CATHARE code is used to model the compressible diphasic flow in the core, primary pumps, Intermediate Exchangers (EI), secondary circuits, . . . . TrioCFD (previously named "Trio_U") is a programmable software developed by CEA for monophasic Computational Fluid Dynamics (CFD) code, based on the open-sourced TRUST platform (TRio_U Software for Thermo-hydraulics) [START_REF]trust-platform[END_REF]. The TrioCFD code is used to model the incompressible monophasic flow in the large hot and cold collectors of pool-type, as well as in the EI primary side and in the in-core inter-wrapper gap regions. The TrioMC code (Modèle Coeur in French) is a TRUST-based specific diphasic flow application dedicated to the thermal-hydraulics computation of the sodium fast reactor's core sub-channels [START_REF] Gerschenfeld | Development and Validation of Multi-Scale Thermal-Hydraulics Calculation Schemes for SFR Applications at CEA[END_REF]. Each one of these codes link the ICOCO library as coupling API. There are integrated into the new code MATHYS (Multiscale ASTRID Thermal-HYdraulics Simulation) [START_REF] Gerschenfeld | Development and Validation of Multi-Scale Thermal-Hydraulics Calculation Schemes for SFR Applications at CEA[END_REF].

Coupled CATHARE/TrioCFD modeling has already been carried out for the simulation of the Phenix reactor [START_REF] Bavière | A first system/CFD coupled simulation of a complete nuclear reactor transient using CATHARE2 and TRIO_U. Preliminary validation on the Phénix Reactor Natural Circulation Test[END_REF][START_REF] Bandini | Assessment of systems codes and their coupling with {CFD} codes in thermal-hydraulic applications to innovative reactors[END_REF][START_REF] Pialla | Overview of the system alone and system/cfd coupled calculations of the {PHENIX} natural circulation test within the {THINS} project[END_REF] or the ASTRID reactor. The used CATHARE/TrioCFD coupling method relies on a partial overlap of the computation domain [START_REF] Bavière | A first system/CFD coupled simulation of a complete nuclear reactor transient using CATHARE2 and TRIO_U. Preliminary validation on the Phénix Reactor Natural Circulation Test[END_REF][START_REF] Pialla | Overview of the system alone and system/cfd coupled calculations of the {PHENIX} natural circulation test within the {THINS} project[END_REF]. The CATHARE code calculates the entire reactor system (complete domain) and TrioCFD code calculates only a portion of this domain (the overlapping sub-domain), cf. Fig. 1. The computation of the sub-domain requires coupled Boundary Conditions (BCs) provided by CATHARE. The computation of the complete domain requires a feedback method to take account of the relevant information resulting from the computation of the sub-domain. These exchanges of information are performed sequentially according to a coupling period to be defined, cf. Fig. 2. Generally speaking, temperature and mass flow are imposed on the TrioCFD input 2D-surfaces and the pressure on the exit 2D-surfaces. Since the CATHARE code has in 1D-variables as the phase enthalpies, densities and velocities, it is necessary to apply some transformations before fixing the CFD sub-domain BCs (for instance, we can cite the conversions of the enthalpy into temperature and the 1D fields into 2D fields).

Figure 1: Schematic diagram of the system thermal-hydraulics/CFD coupling by partial domain overlap [START_REF] Bavière | A first system/CFD coupled simulation of a complete nuclear reactor transient using CATHARE2 and TRIO_U. Preliminary validation on the Phénix Reactor Natural Circulation Test[END_REF].

Concerning the today feedback method, two terms are locally inserted in the CATHARE computations, cf. Fig. 2. For the enthalpy balance, the convected enthalpies (according to the flow direction) are specified on the scalar nodes of the CATHARE 1D-element located at the borders of the overlap sub-domain. It is done using the values of the CFD code. To obtain these enthalpy values, it is also necessary to transform the temperature 2D-fields into enthalpy 1D-fields. For the momentum balance, a pressure feedback method is used. A source/sink term is added to the Right Hand Side (RHS) of the CATHARE momentum balance equation at the vector nodes located at the boundaries of the overlap sub-domain. This correction is build from the pressure difference between TrioCFD and CATHARE computations at several nodes. Its purpose is to locally impose the pressure gradient calculated with the CFD code.

[P (O ref ) -P (O i )] CAT H. -[P (O ref ) -P (O i )] CF D = 0 (1)
with P (O ref ) the pressure at a reference point in the overlap sub-domain and P (O i ) the pressure at the boundary i. For this purpose, a force is modeled, function of this pressure difference, and involved in a iterative process. In the case of a stationary simulation, this force is evaluated at each CATHARE time step. While, in the case of a transient simulation, the iterative process is internal to the CATHARE time step.

In the next sections, we will present a new definition of the feedback method based on the local defect correction method, all the rest being unchanged.

the local correction method in cathare

In this section we present the local defect correction method and its application to the CATHARE code.

The Local Correction Method

The Local Defect Correction (LDC) method is an iterative coupling method due to Hackbusch [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF], using an overlap of local and global domains in a in a local multigrid architecture, cf. Fig. 3.

As an illustration of its principle, let us suppose that we have to solve the following system of linear equations:

A c X c = B c (2)
in the global computation domain (here the CATHARE's one; named coarse grid in Fig. 3) and that there is a solution XCF D in the overlapping local computation sub-domain (here the CFD's one; named fine grid in Fig. 3). The LDC method introduces a local correction of this system of equations via the truncation of the global operator A c : A c XCF D . The corrected system reads: Generally speaking, there already exists several applications of the LDC approach giving satisfactory results in the field of the nuclear power plant simulation. We can cite works concerning the computation of the diphasic-fluid flows in the steam generators [START_REF] Belliard | Local Zoom Computation of Two Phase Flows in Steam Generators using a Local Defect Correction Method[END_REF] or of the thermo-mechanics behavior [START_REF] Barbié | Strategies involving the local defect correction multi-level refinement method for solving three-dimensional linear elastic problems[END_REF] of the nuclear fuel. Moreover, its implementation is relatively simple and involves few intrusions into pre-existing codes.

A c X c = B c -χ(B c -A c XCF D ). (3

Application to the CATHARE Code

We apply the LDC method to the six balance equations of the CATHARE code: phasic mass, phasic momentum and phasic energy It can be seen as a generalization to the whole CATHARE unknown vector (phasic velocities, phasic enthalpies, void fraction, pressure) of the pressure correction previously described. This can potentially be more robust.

Let us begin with some recalls concerning the numerical methods in CATHARE. At each time step the non-linear balance equation are implicitly solved by an iterative Newton algorithm:

J n c δX n+1 c = R c (X n c ) (4)
where

J n c = [ ∂Ac ∂Xc ] n is the Jacobian matrix, δX n+1 c the vector of increments X n+1 c = X n c + δX n+1 c , R c (X n c ) = B c -A c X n
c the residual and n the iteration index. In practice, the solved system solely concerns the unknowns at junctions instead of the whole unknown set:

• the Jacobian matrix and the residue associated with each element are computed,

• the internal degrees of freedom X 1 are eliminated for each element (Schur's complement method), cf.Fig. 4 where the Jacobian matrix J n c entries are denoted by A ii and the residue vector R c entries by

S i , A 11 A 12 A 21 A 22 X 1 X 2 = S 1 S 2 (5)
in order to obtain a system of equations reduced to the unknowns at junctions X 2 ,

(A 22 -A 21 A -1 11 A 12 )X 2 = S 2 -A 21 A -1 11 S 1 . (6)
Let us recall that at each junction j, (X 2 ) j gather the junction variables presented in Fig. 5.

• After the assembly of the junction-reduced Jacobian for the whole set of elements, the arrising system is solved for the junction unknowns.

• The internal unknowns are regenerated from the junction unknowns.

• The algorithm convergence is tested. Implementing the LDC method in CATHARE is easy as only adding of RHS terms χR c ( XCF D ) is required in the Newton equation. Indeed the Newton's version of the LDC algorithm reads:

J n c δX n+1 c = R c (X n c ) -χR c ( XCF D ). (7)
Once the CATHARE code is equipped with the LDC method, we can face the coupling with any external code provided that the ICOCO API is linked and that restricted unknowns XCF D are computed. The coupling algorithm is quite unchanged: The place of these local corrections (defined by the χ function) depends on the type of the CATHARE elements and the nature of the BCs imposed on the CFD sub-domains. There can be at internal degrees of freedom or junctions, cf. Fig. 5 for a CATHARE AXIAL element (1D-element). On Fig. 5, vector nodes are denoted by and scalar nodes by . Only two internal nodes exist for a CATHARE VOLUME element (0D-element with two sub-volumes), but several ones can be addressed for a CATHARE AXIAL element allowing a better coherency between CFD results and CATHARE ones. It can be viewed as a two-grid one-V -cycle algorithm. A full multigrid algorithm involving several V -cycles would need to introduce the local problem resolution inside the Newton algorithm. . .

verification/validation tests on the Vertical-Canon experiment

We present the verification (is the code doing what we expect ?) and the validation (is the code really computing the CATHARE coarse-mesh overlapped domain as the CATHARE fine-mesh one ?) on the diphasic-fluid fast depressurization test known as Vertical Canon test [START_REF] Guiot | CATHARE 2 v1.5qr6 qualification on vertical CANON with internal structures[END_REF].

The Vertical-Canon blowdown tests were carried out to study the blowdown of a vertical pipe initially full of under-saturated water under a strong pressure (here about 5.7 Mpa for test #22; noted CNV22). The blowdown was the consequence of the opening of a small break at the top of the pipe. The test section consists of a stainless steel (type 316L) pipe of 0.100 m internal diameter and 4.483 m total length. The pipe is vertical and the wall thickness is equal to 12.5 mm. It is insulated to prevent the heat exchange with the environment. The bottom end of the test section is closed. The upper head is followed by a converging nozzle of 7 degrees and 0.098 m length. The junction between the pipe and the nozzle is made by a toric piece with a radius depending on the break diameter. The nozzle is followed by a straight pipe of 4 mm length with an internal diameter equal to the wished break internal diameter. The closed test section is heated and pressurized up to the wanted initial fluid pressure and temperature. Then, the break is opened and blowdown occurs. It is a good test for the LDC algorithm as all the variables in the unknown vector play an important role. From a physical point of view (cf. Fig. 7 to 13), in the short term (about the 5-10 first seconds), a mass boiling phase is induced by the capacity depressurization, followed by a gravity-driven phase during which a stratification occurs, the steam moving upward and the liquid downward.

Verification on CATHARE unitary tests

For verification purpose, we consider what we call unitary tests specifying fixed user-given LDC restricted variables XCF D,user on parts of the circuit. The goals are to asses the LDC method's coding and to test its robustness. Among the several unitary tests carried out on the CNV22 simulation, we present here a particular one for which the following vector XCF D,user is fixed on the scalar and vector nodes of the element CANAL1:

• P = 4.3 MPa, • H l = 1000.514 kJ/kg, • H g = 2798.376 kJ/kg, • V l = 10 -7 m/s, • V g = 10 -7 m/s, • α = 10 -5 ,
and on the first vector node of the element CANAL2:

• V l = 10 -7 m/s, • V g = 10 -7 m/s.
Let us notice that the values of the saturated liquid and gas enthalpies at 4.3 MPa are 1109.276 kJ/kg and 2798.376 kJ/kg respectively. Then, the restricted vector defines a static monophasic fluid with undersaturated liquid and residual saturated vapor. In this test, the corrected residuals are located at the same nodes that the restricted vectors. Fig. 7 gathers various pressure and void fraction profiles on the whole computation domain (CANAL1 + CANAL2 + CANAL3) at times: 0, 15, 40 and 100 sec. after the break opening. Indeed as the phasic velocities are null on the CANAL1/CANAL2 junction, the local defect corrected case can be viewed as the standard simulation of a short version of the straight pipe (without CANAL1). The full comparison of the LDC's results (solid lines) to the standard-code short-pipe results (dashed line) gives confidence in the implementation of the method, cf. Figs 8 to 9.

Validation on Coupled CATHARE models

For validation purpose, we consider a coupled model for which the central part of the global circuit describing the straight pipe is locally refined using an other CATHARE circuit. The initial 23-cells mesh of the global-circuit central part (CANAL2) is refined using a ratio of 5 (CANAL2_REF defining the local circuit). Each circuit is sequentially computed following the coupled algorithm 1-2-3-4 presented in Section 2.2. Let us recall that the time step is dynamically computed in CATHARE for each circuit. Here, we retain the smallest one as common time step. Figs 10 to 12 present various profiles of the main variables for several times equal to 0, 1, 2 and 5 seconds. Results from the global circuit (symbol X) and from the local circuit (symbol +) are greatly in coherence. Also it is the case comparing with standard-code reference solution refined in space (on CANAL2) and time (solid line). But on the overlapped sub-domain, the global circuit's results are the closest with the reference ones. In comparison the results obtained using the pressure feedback method (1) solely converge toward the refined in space approximation (i.e. not in time), cf. Fig. 13. This is due to the time evolution of the added LDC term in the Newton equation leading to a steeper RHS and smaller time step (mean value: 0.02 instead of 0.3 second). In the case of the pressure feedback approach, we should impose a smaller CATHARE time step to obtained similar results to the LDC ones. Let us notice that the CPU time of the LDC computation is about 10 times greater than the pressurefeedback computation time. But it is twice as small as the computation time of the standard-code solution refined in space and time.

conclusion and perspectives

In the framework of multi-scale computations of the thermal-hydraulics behavior of the sodium fast reactors as the ASTRID project one, we have presented the implementation of the local defect correction method in the French system-scale safety reference code CATHARE. Some illustrations of the coupling capacities of the LDC method have been produced on verification/validation tests concerning compressible diphasic fluid behavior during a fast depressurization transient. On-going work concerns the realization of CFD/system-scale coupling for monophasic fluid using TrioCFD and CATHARE. But the performances shown by coupling method for the diphasic tests give confidence for these incompressible monophasic fluid tests. As a first result, Fig. 14 shows the time evolution of an enthalpy step (computed by TrioCFD) traveling outside the CATHARE computational domain.

Figure 2 :

 2 Figure 2: Schematic diagram of the CATHARE/TrioCFD coupling: CFD boundary conditions and feedback method in CATHARE [4].

Figure 3 :

 3 Figure 3: Schematic diagram of a local two-grid architecture [9].

)

  In this equation, XCF D is provided by the restriction of the local (CFD) sub-domain unknowns X CF D to the part of the global computation domain (CATHARE's one) overlapped by the local sub-domain. Here χ is the characteristic function of this part of the global computation domain taking a value equals to 1 in the overlapping and 0 elsewhere. Considering a set of embedded locally overlapped sub-domains, iterative V or W -multigrid cycles are used to solve the balance equations inside each current time step.

Figure 4 :

 4 Figure 4: Jacobian matrix composition: splitting of the matrix entries into internal-internal and internaljunction parts.

1 .

 1 coupled BCs are defined by extension of the CATHARE's variables X c into the boundaries of the overlapping CFD's sub-domain, 2. a (partial or not) solving of the CFD problem is performed obtaining the solution X CF D , 3. the CFD's variables X CF D are restricted on the CATHARE's overlapped sub-domain: XCF D , 4. the system-scale residual (B c -A c XCF D ) is computed using the restricted variables XCF D and locally added to the RHS of the Newton equation χ(B c -A c XCF D ).

Figure 5 :

 5 Figure 5: LDC corrections on junction and internal nodes of an CATHARE AXIAL element.

Figure 6 :

 6 Figure 6: Meshing for the LDC tests on the Vertical-Canon experiment.

  P, ALPHA, HL, HG, VG, VL)_ldc given from 0.0 to 2.02m; ----: std (short) P, ALPHA, HL, HG, VG, VL)_ldc given from 0.0 to 2.02m; ----: std (short) (b) Void fraction

Figure 7 :

 7 Figure 7: Pressure and void fraction profiles for verification test "short-pipe Vertical-Canon": LDC applied on CANAL1 (solid lines) versus standard-code short-pipe simulation (dashed line).

  P, ALPHA, HL, HG, VG, VL)_ldc given from 0.0 to 2.02m; ----: std (short) P, ALPHA, HL, HG, VG, VL)_ldc given from 0.0 to 2.02m; ----: std (short) (b) Steam velocity

Figure 8 :

 8 Figure8: Liquid and steam velocity profiles for verification test "short-pipe Vertical-Canon": LDC applied on CANAL1 (solid lines) versus standard-code short-pipe simulation (dashed line).

  experiment -LDC -____ : (P, ALPHA, HL, HG, VG, VL)_ldc given from 0.0 to 2.02m; ----: std (short) experiment -LDC -____ : (P, ALPHA, HL, HG, VL, VG)_ldc given from 0.0 to 2.02m; ----: std (short) (b) Steam enthalpy

Figure 9 :

 9 Figure9: Liquid and steam enthalpy profiles for verification test "short-pipe Vertical-Canon": LDC applied on CANAL1 (solid lines) versus standard-code short-pipe simulation (dashed line).

Figure 10 :

 10 Figure 10: Pressure and void fraction profiles for validation on Vertical-Canon: LDC applied on CANAL2 (crosses) versus standard-code reference solution (solid line).

  LDC computation Gas velocity (m/s) vs elevation (m) (b) Steam velocity

Figure 11 :

 11 Figure 11: Liquid and steam velocity profiles for validation on Vertical-Canon: LDC applied on CANAL2 (crosses) versus standard-code reference solution (solid line).

Figure 12 :

 12 Figure 12: Liquid and steam enthalpy profiles for validation on Vertical-Canon: LDC applied on CANAL2 (crosses) versus standard-code reference solution (solid line).

  dpext computation Pressure (1.0D5 Pa) vs elevation (m) dpext computation Void fraction vs elevation (m) (b) Void fraction

Figure 13 :

 13 Figure13: Pressure and void fraction profiles for validation on Vertical-Canon: pressure feedback method applied on CANAL2 (crosses) versus standard-code reference solution (solid line).

* Corresponding author