Ulrich Bieder
email: ulrich.bieder@cea.fr

Clarisse Genrault

Pierre Ledac

Hydraulic Forces Acting on Full Cross Section Fuel Assemblies with 17×17 Fuel Rods

The pressure and velocity distribution within and between full cross-section fuel assemblies with generic but realistic mixing grids are analyzed by CFD using high performance computing (HPC). The reference fuel assembly consists of a bundle of 17×17 rods, including fuel rods, control rod guide tubes and a mixing grid with split type mixing vanes. The axial length of the analyzed domain spans approximately the interval between two successive mixing grids. Only hydraulic effects of the mixing grids on the behavior of the flow are analyzed, power input and heat transfer are not considered

The presented results are based on a well-validated calculation procedure. The intra-assembly exchanges and pressure fields were analyzed by modelling one single assembly, whereas the inter-assembly exchanges and pressure fields were analyzed by modelling a group of three horizontally aligned fuel assemblies, each build of 17×17 rods and a mixing grid. In the multi-assembly test cases, the water gaps that separate the central assembly from the neighboring fuel assembles vary between 2 mm and 10 mmm. In one additional test case, the rods of the central fuel assembly were slightly inclined in order to treat deformed assemblies.

For a fuel assembly that is surrounded by water gaps of 2 mm, the total force (pressure and shear force) acting in flow direction upon one mixing grid is about 230 N. The total force working on the grid in cross flow direction is small. The prediction of this horizontal force can be falsified by inappropriate boundary conditions as symmetry. The total force acting upon the rods is about 160 N in flow direction and small in crossflow direction.

Variations in the width of two oppositely located water gaps significantly influence the pressure force that is acting on the mixing grid in horizontal direction. Depending on the difference in the width of the water gaps, pressure forces of 39 N and 24 N were calculated for gap couples of 2 mm and 5 mm, and of 3 mm and 7 mm, respectively. This horizontally acting pressure force (Venturi force) tries to push the assembly in the direction to homogenize the water gap widths (self-stabilizing effect).

The inclination of the rods of the central assembly does not lead to significant forces in cross flow direction. In order to assure this result, a calculation with very fine meshing (1000 million tetrahedrons) was performed.

INTRODUCTION

The reactor core of a PWR (Pressurized Water Reactor) is assembled from arrays of nuclear fuel bundles, which in turn are typically composed by a group of 15×15 or 17×17 bars, more precisely of fuel rods and guide tubes. Support and mixing grids are positioned within the assembly at specific axial distances to maintain the lateral spacing of fuel rods. Since turbulent, forced convection dominates the transport of the thermal energy from the surfaces of the fuel rods to the bulk of the fluid, mixing vanes are placed on the top edges of the mixing grid lattice in order to enhance turbulence and thus heat transfer. The cooling water flows mainly in the direction parallel to the rods. Due to the mixing vanes, cross flow occurs in planes orthonormal to the axial flow direction. This cross flow influences the horizontal fluid mixing as well as the pressure distribution in fuel assemblies and therefore the pressure forces on the grids and the fuel rods. [START_REF] Karoutas | 3-D flow analyses for design of nuclear fuel spacers[END_REF] realized one of the first CFD simulations of turbulent flow in a fuel bundle with mixing vanes by using Reynolds Averaged Navier Stokes Equations (RANS). Only half of a single sub channel was simulated using symmetry boundary conditions. The calculations were validated by comparison to experiments. W.K. [START_REF] In | Flow analysis for optimum design of mixing vane in PWR fuel assembly[END_REF] have optimized vane geometries by using RANS calculation (k-ε). [START_REF] Navarro | Evaluation of a numeric procedure for flow simulation of a 5×5 PWR rod bundle with a mixing vane spacer[END_REF] improved these results using a computational domain of a 5×5 rods. Such 5×5 rod bundle domains have been extensively analyzed in various studies. As an example, the MATIS_H experiment [START_REF] Smith | Report of the OECD/NEA KAERI rod bundle CFD benchmark exercise[END_REF] are cited here that was analyzed among others by [START_REF] Rashkovan | Examination of Separate-Effect and Integral Phenomena Within a Grid Spacer with Mixing Vanes: Results for the MATiS-H OECD_NEA Benchmark Exercise[END_REF] and [START_REF] Mikuž | URANS prediction of flow fluctuations in rod bundle with split-type spacer grid[END_REF]. [START_REF] Bieder | LES analysis of the flow in a simplified PWR assembly with mixing Grid[END_REF] analyzed the AGATE experiments, [START_REF] Wells | Computational Fluid Dynamics Benchmark of High Fidelity Rod Bundle Experiments Industry Round Robin Phase 2 -Rod Bundle with Mixing Vane Grids[END_REF] and [START_REF] Kang | Computational fluid dynamics (CFD) round robin benchmark for a pressurized water reactor (PWR) rod bundle[END_REF] studied the NESTOR experiments. Other 5x5 rods experiments were investigated for example by [START_REF] Conner | Hydraulic benchmark data for PWR mixing vane grid[END_REF] or [START_REF] Nguyen | Stereoscopic particle image velocimetry measurements of flow in a rod bundle with a spacer grid and mixing vanes at a low Reynolds number[END_REF].

Due to the important increase of computational power in recent years, first CFD analysis of full size fuel assemblies have been published. [START_REF] Li | Methods of simulating large-scale rod bundle and application to a 17×17 fuel assembly with mixing vane spacer grid[END_REF] simulated a 17×17 rod bundle with a single mixing grid. [START_REF] Yan | Influence of Spacer Grid Outer Strap on Fuel Assembly Thermal Hydraulic Performance[END_REF] have analyzed by CFD the influence of the mixing grid outer straps (protective tabs) on the flow within adjacent assemblies. [START_REF] Mikuž | Low resolution modelling of mixing phenomena in PWR fuel assembly with split type mixing grid[END_REF] analyzed mixing phenomena in a heated PWR fuel assembly of 15×15 rods with split type mixing grids. They assumed all rods being similar and used a so-called low-resolution approach with porous media and momentum models to mimic the effect of spacers. [START_REF] Mikuž | Low resolution modelling of mixing phenomena in PWR fuel assemblies[END_REF] also presented the low-resolution CFD model. The model reproduces similar pressure and temperature distribution in the PWR rod bundle as it has been observed in reference CFD simulations with reduced resolution. [START_REF] Bieder | CFD analysis of intra and inter fuel assembly mixing[END_REF] have analyzed the cross flow in an assembly of 17×17 rods (intraassembly flow) as well as in a cluster of 2×2 fuel assemblies (inter-assembly flow), where each assembly consists of 17×17 rods.

To the knowledge of the authors, hydraulic forces acting on real cross-section mixing grids in 17×17 rods assemblies have not been published to date. The intra-assembly pressure field is analyzed here by modelling one single assembly. Inter-assembly pressure fields are analyzed by modelling a bundle of three horizontally aligned fuel assemblies, each build of 17×17 rods and a mixing grid. Finally, the influence of inclined rods on flow field and hydraulic forces is discussed. Only hydraulic effects of the mixing grids on the behavior of the flow are analyzed, power input and heat transfer are not considered yet.

THE GENERIC FUEL ASSEMBLY

The calculation domain comprises a part of a PWR fuel assembly (FA) of 17×17 fuel rods including a mixing grid. Mixing grids maintain the lateral spacing between fuel rods and support the fuel rods at intervals along their length. The mixing grid consists of an egg-crate arrangement of interlocking straps. The straps contain spring fingers and dimples for fuel rod support as well as coolant mixing vanes [START_REF] Weihermiller | Batelle Topical Report[END_REF].

Guide taps on the outer straps assure the unhindered axial movement of the fuel assembly during charge and discharge of the reactor core. A typical fuel assembly as well as the used mixing grid are given in Fig. 1 The resulting CAD model of the generic mixing grid is shown in more detail in Fig. 2. This CAD model was developed in SALOME1 platform; it is parametrized by PYTHON scripts in order to have the possibility to modify easily the geometry of the mixing grid for further studies. A zoom on the CAD-model of fuel rods, guide tube, springs, dimples, lattice strap and protective tabs is given in Fig. 2a. Guide tubes have a larger diameter than fuel rods. They are not fastened by springs and dimples but by guide thimbles that are joined to the straps. The protective tabs are located on the outer lattice straps, which have a slightly higher width than the inner straps [START_REF] Weihermiller | Batelle Topical Report[END_REF]. The geometry of the springs is simplified; they are not modelled as bands. The tangential contact between dimple and fuel rod was avoided by increasing their contact zone. Fig. 2b shows a generic form of split type mixing vanes with a triangular central hole. Drop wise holes as used by [START_REF] Delafontaine | Simulation of unsteady fluid forces on a single rod downstream of mixing grid cell[END_REF] are used for the calculations presented here. The angle between grid strap and vane is 33°. The mixing vanes are arranged on the straps to initiate diagonal cross flow in order to increase the flow mixing between sub-channels. Protective tabs are located at the lateral grid straps, all oriented inside the mixing grid. A summary of the geometry of the FA is given below;

• Rod diameter: 9.5 mm • Rod-to-rod pitch:

13.0 mm • Rod-to-rod gap:

3.5 mm • Hydraulic diameter of a sub channel:

13.15 mm Three arrangements of the generic, full cross-section fuel assembly of 17×17 rods were investigated as shown in Fig. 3: one single FA (Fig. 3a) and two groups of three FAs that are aligned in one row. In the case of groups of FAs, two configurations are modeled; three FAs, each with straight rods (Fig. 3b) where the width of the water gaps d1 and d2 can vary between 2 mm and 10 mm, and three FAs where the central FA has partly inclined rods. In the case of a FA with inclined rods, the variation of the water gaps is shown in Fig. 3c. The assemblies are numbered from left to right according to A, B and C assembly. The orientation of FA's mixing grid is always identical; hence, the grids are never rotated horizontally.

THE NUMERICAL MODEL

The TrioCFD code [START_REF] Angeli | Overview of the Trio_U code: Main features, V&V procedures and typical applications to engineering[END_REF] is used for all calculations presented here. TrioCFD is a general CFD code, developed at CEA, dedicated to simulate incompressible or dilatable flows in nuclear engineering applications. The code is open source and uses efficiently the capabilities of massively parallel HPC machines. Bieder and Genrault (2020) have already described the numerical model and give information on the meshing procedure, mesh sensitivity studies and on mesh convergence as well as on the validation of the model. This information is not repeated here; rather the focus is directed to the application of the model.

In the calculations presented in this study, the fluid is assumed Newtonian and incompressible. As only isothermal flow is treated, buoyancy effects are not taken into account. The instantaneous velocity of such a fluid is expressed by the equations of mass-and momentum conservation. The equations can be found in in standard CFD books as e.g. [START_REF] Pope | Turbulent Flows[END_REF].

Turbulence modelling

The Reynolds number of the flow in the analyzed fuel assembly is about 5×10 5 . Turbulence is modelled by unsteady Reynolds averaged Navier-Stokes equations (U-RANS). In Reynolds-averaged turbulence modelling, an averaging operator is applied to the instantaneous velocity. The non-linearity of the Navier-Stokes equations gives rise to Reynolds stress terms that are modeled by turbulence models. In this study, the Reynolds stresses are modeled according to Boussinesq's concept of eddy-viscosity, which assumes that the Reynolds stresses are aligned with the main strain rates. The high Reynolds number k-ε model [START_REF] Pope | Turbulent Flows[END_REF] is used to calculate the turbulent viscosity. This turbulence model has been selected because of the profound validation effort that has been done at CEA. In fact, for the simulation of the flow in fuel assemblies with mixing grids the use of the k-ε model has been validated by means of the AGATE experiments of CEA [START_REF] Bieder | LES analysis of the flow in a simplified PWR assembly with mixing Grid[END_REF] and the MATHIS_H experiment of KAERI [START_REF] Bieder | Analysis of the Flow Down-and Upwind of Split-type Mixing Vanes[END_REF].

Discretization method

TrioCFD [START_REF] Angeli | Overview of the Trio_U code: Main features, V&V procedures and typical applications to engineering[END_REF] uses a finite volume based finite element approach on tetrahedral cells to integrate in conservative form all conservation equations over the control volumes belonging to the calculation domain. As in the classical Crouzeix-Raviart element, both vector and scalar quantities are located in the centers of the faces [START_REF] Angeli | FVCA8 Benchmark for the Stokes and Navier-Stokes Equations with the TrioCFD Code-Benchmark Session[END_REF]. The pressure, however, is located in the vertices and at the center of gravity of a tetrahedral element [START_REF] Angeli | FVCA8 Benchmark for the Stokes and Navier-Stokes Equations with the TrioCFD Code-Benchmark Session[END_REF]. This discretization leads to very good pressure/velocity coupling. Along this staggered mesh arrangement, the unknowns, that is the vector and scalar values, are expressed using nonconforming linear shape-functions (P1-non-conforming). The shape function for the pressure is constant for the center of the element (P0) and linear for the vertices (P1). [START_REF] Angeli | FVCA8 Benchmark for the Stokes and Navier-Stokes Equations with the TrioCFD Code-Benchmark Session[END_REF] give more information on the discretization method.

Numerical scheme, solution method and boundary conditions

The 1 st order Euler backward time marching scheme was used to integrate the momentum equations. A fully implicit predictor corrector method assures mass and momentum conservation. [START_REF] Guermond | On the Stability and Convergence of Projection Methods Based on Pressure Poisson Equation[END_REF] have demonstrated the stability of this scheme. The generalized minimal residual method (GMRES) with Jacobi preconditioning is used to predict the intermediate velocity that is not divergence free. In the corrector step, this velocity is projected with a conjugated gradient method (CG) with symmetric successive over-relaxation (SSOR) preconditioning in the divergence free space by solving Poisson's equation of the pressure. The CG algorithm of the PETSc2 library is used.

The discretization in space of the Navier-Stokes equations is realized by a second order upwind scheme for the convection term and a second order centered scheme for the diffusion term. The pressure gradient follows a first order approximation.

Dirichlet boundary conditions were used at inflow faces. Fully developed turbulent flow in a bare rod bundle without mixing grids is assumed at the inlet [START_REF] Bieder | LES analysis of the flow in a simplified PWR assembly with mixing Grid[END_REF]. Neumann conditions with an imposed pressure are applied at outflow faces. The general wall function of [START_REF] Reichardt | Vollständige Darstelling der turbulenten Geschwindigkeitsverteilung in glatten Leitungen[END_REF] is used to model momentum exchange between walls and fluid. A slip wall boundary condition (symmetry) is applied at the boundaries that surround laterally the rod bundle in order to mimic free flow conditions in the water gaps between adjacent FAs. Table 1 summarizes the complete numerical scheme.

Meshing of the calculation domains

The CAD model of the mixing grid was transferred in STEP format to the commercial mesh generator ICEMCFD ® . A pure tetrahedral meshing was created around the mixing grid in the region called Meshed Zone that has an axial extension of 43 mm (see Fig. 4). One mm of the surrounding, water gaps are also meshed. As discussed in more detail by [START_REF] Bieder | CFD analysis of intra and inter fuel assembly mixing[END_REF], the best performance concerning the resolution time of one time step is achieved with the Delaunay meshing method and when using a homogeneous mesh size in the whole Meshed Zone; local mesh refinement was not applied. The meshing used in this study has a reference size of 1 mm, what led to about 50 million tetrahedral cells in the Meshed Zone. The height of a tetrahedron is about 0.5 mm with this reference size. As the velocity is discretized in the center of the faces of the tetrahedrons, the distance between two calculation points is about four times smaller than the reference size, that is about 0.25 mm. The meshing of the lateral faces of the Meshed Zone is periodic in x and y direction. The two-step meshing procedure (initial surface meshing followed by volume meshing) leads to tetrahedrons on the walls that are very similar in size. Therefore, the y + values of the wall meshes are very similar and the mesh type will not affect the near wall flow significantly as long as wall functions are used.

Reference meshing of a single fuel assembly

The surface mesh of the axial boundaries of the Meshed Zone, which is represented by triangles, was extruded uniformly in z-direction, more precisely -30 mm in upstream direction and +300 mm in downstream direction. The resulting extruded triangles form prismatic cells. These zones made of prisms are called Extruded Zones. In order to create a meshing made exclusively of tetrahedral cells, each created prismatic cell is divided into 3 tetrahedrons. In the reference meshing, each extruded tetrahedron has an axial length of 3 mm. The resulting calculation domain of one assembly of 17×17 rods is shown in Fig. 4, including Meshed Zone, the Extruded Zones and the inlet and outlet plane.

The total number of cells of the complete calculation domain of one FA is about 130 million tetrahedrons. The number of velocity control volumes is twice as high as the number of tetrahedrons, namely about 260 million. As the pressure is discretized in the center and the vertices of an element, the number of pressure control volumes, that is the control volumes in which the divergence free condition is assured, is about 1.125 times higher than the number of tetrahedrons, namely 146 million. Mesh independency of the velocity field was studied as described in more detail by [START_REF] Bieder | LES analysis of the flow in a simplified PWR assembly with mixing Grid[END_REF] and [START_REF] Bieder | CFD analysis of intra and inter fuel assembly mixing[END_REF].

Reference meshing of three straight fuel assemblies in line

In order to gather clusters of several fuel assemblies of 17×17 rods, the meshing of the lateral faces of the Meshed Zone is periodic in x and y direction. Hence, clusters of three FAs that are arranged in a line as shown in Fig. 3 can be compounded easily by using the mesh manipulation capacities of TrioCFD. The periodic faces of two domains of the Meshed Zone can be merged directly, what leads to a water gap of 2 mm between the mixing grids. For wider water gaps, the meshing of the lateral face of one Meshed Zone can be extruded befor merging two adjacent mixing grids. In this way, several Meshed Zones can be merged with different water gaps between the mixing grids. As in the case of a single FA, the surface mesh of the axial boundaries of the complete Meshed Zone was extruded uniformly in z-direction for -30 mm in upstream direction and +300 mm in downstream direction. Each prisms layer has an axial length of 5 mm. After splitting the prismatic cells in tetrahedrons, the complete calculation domain of three aligned FA with water gaps of different width include about 390 million tetrahedrons.

Fine mesh of three straight fuel assemblies in line

In the reference meshing, the tetrahedrons in the Extruded Zones have an axial length of 5 mm. This leads in these zones to stretched tetrahedrons that are longer in axial flow direction and shorter in cross flow direction. Such stretched meshes might overestimate flow phenomena that are aligned with the fuel rods and might underestimate flow effects in cross flow direction. To assure that the meshing does not favorite a certain flow direction, tetrahedral meshes that are not stretched were created in the Extruded Zones. Therefore, the number of extruded layers was increased in these zones by a factor 7.15. The height of a tetrahedron is then about 0.7 mm, what is slightly above the tetrahedron's height of 0.5 mm in the Meshed Zone. This fine mesh of 1.030 billion tetrahedrons was used for mesh convergence tests.

Fine mesh of one inclined assembly between two straight assemblies

Fig. 3c illustrates the case of slightly inclined fuel rods in the center assembly. In this case, it is very important that effects related to anisotropic meshing not dominate small forces in cross flow direction. Thus, for this configuration, it is mandatory to use the fine mesh for the analysis since meshing effects cannot be excluded a priori for the anisotropic reference mesh. The inclination is realized by shearing the mesh of the Extruded Zone of the center assembly. The resulting slightly elliptical form of the rods is assumed negligible. Before shearing the mesh, a fine mesh with upward rods was build, with water gaps of 3 mm between assemblies A and B and 7 mm between B and C. The shearing process opened in an axial length of 450 mm gap A-B from 3 mm to 7 mm and closed accordingly gap B-C from 7 mm to 3 mm (see Fig. 3c).

It was tested to create the meshing from a CAD model of the full calculation domain with more realistically inclined rods that is with inclined rods that pass across the mixing grid. Modifying the CAD model accordingly was not difficult due to the applied modelling procedures with python scripts. However, meshing the whole domain with a reference size of 1 mm was not possible with ICEMCFD due to memory limits of the used workstations. The authors did also not succeed to do the meshing via the creation of several sub-meshes, which can be merged to the compete domain. However, our failure does not exclude that such a stepwise meshing procedure is not possible. Finally, it was tested to create the fine mesh by homogeneous mesh refinement. Therefore, as a first step, a coarse mesh of the whole domain with about 125 million tetrahedrons was created. This mesh was then refined in a second step by splitting each tetrahedral cell in eight small tetrahedrons. As it was not possible to project the refined mesh to the initial geometry, the meshing did not represent sufficiently precise the complex geometry of the mixing grid and the curvature of the rods.

Parallelization and used HPC machines

The presented CFD calculations have been carried out by exploiting the parallel calculation capabilities of TrioCFD. The key to run large models is that the code is distributed on the different processor cores, where each core is handling only a part of the data. A master processor core that handles the complete calculation domain is not needed. Using METIS library, each domain is decomposed into several overlapping sub-domains by minimizing the load imbalance of the domain partition. The load imbalance is defined by the product of the maximum cells among sub-domains and the number of sub-domains divided by the total number of cells. Then, all sub-domains were distributed equally among different processor cores, which, by using message passing interface libraries (MPI), communicate mutually only when data transfer is needed. All I/O processes are parallelized with the possibility to read and write from a single file (with HDF5 library) or from distributed files. Post-processing with Visit3 is also available in a parallelized version. When a calculation domain comprises more that about 500 million tetrahedral cells, the integers need to be represented on 64 bits. This feature is included as compiler option. The performance of TrioCFD was evaluated in a preparatory numerical study, which concerns the laminar flow around a circular cylinder. The calculations were performed on Irène-Joliot, the BULL Sequana X1000 supercomputer of the Très Grand Centre de calcul (TGCC) of CEA. The following partition was used:

• Partition AMD with 2292 AMD Rome Epyc 7502 bi-processors nodes at 2.66 GHz with 64 cores per processor. Each node has 256 GB of DDR4 memory. The nodes are interconnected by InfiniBand HDR100 network.

Strong scaling was evaluated on 1.335 billion cells with an Euler explicit time scheme. The pressure solver used is the same as that of the current study (Conjugate Gradient with SSOR preconditioning). The CPU time per time step (dt TS) and the CPU time per conjugate gradient iteration (dt CG) are reported in Table 2. The load imbalance of the different partitions is added to Table 2 as it correlates well with the decrease of parallel efficiency above 12800 cores shown in Fig. 5. Since TrioCFD is a memory bandwidth limited code and in order to spend minimal CPU hours, we run the final calculations on the fine mesh on the SKL partition of Irène-Joliot:

• Partition SKL with 1656 Intel Skylake 8168 bi-processors nodes at 2.7 GHz with 24 cores per processor. Each node has 192 GB of DDR4 memory. The nodes are interconnected by InfiniBand EDR network.

In SKL partition, the network and memory bandwidth per core are 2-3 times higher than in the AMD partition, although fewer cores are available and are slightly slower. The calculations with the reference meshes were performed on Occigen, the BULL B720 supercomputer of the Centre Informatique National de l'Enseignement Supérieur (CINES):

• Partition HSW with 2106 Haswell bi-processor nodes at 2.6 GHz with 24 cores per processor. Each node has 128 GB of DDR4 memory. The nodes are interconnected by InfiniBand FDR network.

The performance of TrioCFD on these HPC computers is summarized in Table 3. The load imbalance was 1.054.

ANALYSIS OF THE HYDRAULIC FORCES IN FUEL ASSEMBLIES

Typical thermal hydraulic operating conditions in a PWR core at nominal power form the basis of this study. These thermal hydraulic conditions close to the core outlet are:

• Total pressure, 16 MPa • Mean coolant temperature, 330°C

The resulting physical properties of the coolant water are:

• Mean coolant density, 648 kg/m 3 • Dynamic viscosity, 7.75×10 -5 Pa s • Mean coolant velocity, 5.35 m/s • Sub-channel Reynolds number, 605000

The physical properties are taken constant in all calculations.

First, the forces acting on one assembly with equal water gap widths are presented. Then the width of the water gaps are varied, and finally the rods of the center fuel assembly were inclined. The pressure force Fp and shear force Fµ that act on mixing grid and rods, respectively, are calculated by integrating the local pressure P and shear force over all solid walls according to:

= -• , (1) = .
(2)

is the unit vector normal to the considered wall, µ the dynamic viscosity and the velocity vector.

Fuel assemblies with identical water gap widths

The main objective of this first study is to provide a better insight into the pressure distribution within a realistic fuel assembly of 17×17 fuel rods with mixing grid.

Analysis of one assembly

It is assumed that a 2 mm water gap exists around the analyzed assembly. In order to simulate these water gaps, symmetry boundary conditions are used at the outer vertical faces of the calculation domain, at a distance of 1 mm from the outer lattice straps of the mixing grids. The whole assembly is shown in a perspective view in Fig. 6a. The spatial distribution of the velocity magnitude is given in color scale. About one-half of the assembly (0.0<y<0.1185) is hidden visually to show the velocity in the rod-to-rod gaps. A predominant axial flow exists upstream of the mixing grid. Downstream of the grid the mixing vanes add significant cross flow components on the main flow. This crossflow disappears about 10 hydraulic diameters downstream of the grid that is after about 0.13 m. Since guide tubes have larger diameters than fuel rods, the velocity in the guide tube-to-rod gaps is reduced significantly compared to the rod-to-rod gaps. The guide tubes also constrain the development of crossflow. The reduced velocity along the guide tubes is well visible in Fig. 6a. Fig. 6b shows the cross-flow velocity in a horizontal plane, cut 0.05 m above the upper edge of the mixing grid straps (z = 0.05 m). The color scale denotes the velocity magnitude and the vectors show the velocity direction. The formation of the diagonal cross flow direction and the reduction of the velocity close to the guide tubes are readily visible. [START_REF] Bieder | CFD analysis of intra and inter fuel assembly mixing[END_REF] The spatial distribution of the pressure is displayed in Fig. 7 in form of 2D plots. Fig. 7a shows the vertical plane at y=0.1185 m and Fig. 7b a horizontal plane at z = 0.0 m. The vertical plane is located in the center of the rod-to-rod gap and traverses a lattice strap. The highest pressure is located upstream of the mixing grid below the water gaps of 2 mm width, which are located in the periphery of the calculation domain. Here the pressure reached a value of approximately 12.9 kPa. The square vertical vanes that are located at the top of the lattice strap are guide tube thimbles. They mark the location of the control rod guide tubes. The pressure gradient in flow direction exceeds largely horizontal pressure variations. Nevertheless, the pressure peaks upstream of the water gaps, which are located in the periphery of the assembly, are readily visible. The low-pressure streaks further downstream the mixing vanes can also be observed easily. They persist along the complete vertical calculation domain.

The horizontal plane is located at the upper edge of the lattice straps. A zoom to the lower left corner of the assembly is presented. Pressure peaks are present on the concave side of the mixing vanes and low-pressure areas are present on their convex side. It is mentioned here that the pressure in the wake of springs is lower than in the wake of dimples. As visible in Fig. 7b, these two design features lead to significant pressure differences within each grid lattice. A pressure loss coefficient Ks was evaluated in the NESTOR round robin benchmark [START_REF] Wells | Computational Fluid Dynamics Benchmark of High Fidelity Rod Bundle Experiments Industry Round Robin Phase 2 -Rod Bundle with Mixing Vane Grids[END_REF] from the MANIVEL experiments. For one mixing grid span and a bundle of 5×5 rods, a pressure loss correlation was derived for Re < 100,000, where l is length of one span: . ) . ( 4) [START_REF] Benhamadouche | Pressure Drop Predictions using CODE_SATURNE in NESTOR CFD Benchmark[END_REF] has shown that this correlation also reproduces the pressure losses of one grid span at higher Reynolds numbers. Correlation and calculations with Reynolds stress turbulence models where compared up to Re=500,000. The calculations underestimate at high Reynolds numbers the correlation by about 10%.

= • • • , ( 3
) = 7.2701 • %& '(
In this study, for Re = 605,000, the pressure loss is ∆P/∆z ≈ 30.5 kPa/m, calculated for an axial span length of ∆z ≈ 0.41 m. The pressure loss of the NESTOR correlation (eqs.( 3) and ( 4)) is ∆P/∆z ≈ 34.2 kPa/m, calculated for an axial span length of about ∆z ≈ 0.34 m. This comparison shows that the axial pressure gradient calculated in this study for full cross-section mixing grids is in accordance to [START_REF] Benhamadouche | Pressure Drop Predictions using CODE_SATURNE in NESTOR CFD Benchmark[END_REF] and the NESTOR benchmark [START_REF] Wells | Computational Fluid Dynamics Benchmark of High Fidelity Rod Bundle Experiments Industry Round Robin Phase 2 -Rod Bundle with Mixing Vane Grids[END_REF].

The calculated pressure and shear forces (eqs(1) and ( 2)), which act on the mixing grid as well as on fuel rods and guide rubes, are summarized in Table 3. Here, Fx and Fy denote the horizontal forces in cross-flow direction and Fz the vertical force in flow direction. The vertical pressure force represents to the force for one axial grid distance. For the mixing grid, the viscous forces are small compared to pressure forces, although a noteworthy viscose force exists in flow direction. Concerning the pressure forces, the component in flow direction is significantly higher than the components in cross-flow direction. Further, the two horizontal pressure forces are not equal due to the unsymmetrical geometrical arrangement of mixing vanes, springs and dimples within the mixing grid. For fuel rods and guide tubes, the viscose force in flow direction dominates all other forces. The pressure forces acting in horizontal direction are somewhat higher than the corresponding viscose forces.

Analysis of three assemblies in line

The validity of the symmetry boundary conditions for the prediction of pressure forces as used in chapter 4.1.1 is reassessed. Therefore, a more general modeling approach was introduced, in which three FAs were gathered. Each assembly consists of 17×17 rods. The assemblies are arranged in line along the x-axis as shown in Fig. 3b. It is assumed that the water gaps, i.e. the gap between two adjacent the FAs, have a width of d1 = d2 = 2 mm. Table 4 summarizes the calculated pressure and shear forces, which act on the mixing grid as well as on fuel rods and guide rubes of the central assembly. Both the shear forces and the pressure forces in flow direction (z-direction) are very similar in the two modelling approaches. The symmetry boundary condition, which is used at the peripheral faces of the single assembly, leads to an overestimation of the pressure force in x-direction compared to the corresponding force acting of the grid of the central assembly. It seems that pressure forces calculated with symmetry boundary conditions can be used only with caution. It is thus expected that the calculated pressure forces in y-direction are also overestimated.

In summary, to give orders of magnitudes for the total force (Fp+Fµ) acting upon the mixing grid in flow direction is about 230 N. Assuming that only the forces in x-direction of the central assembly are representative for cross flow forces, the total force acting in cross flow direction seems to be small. The total force acting upon the rods is 170 N in flow direction and negligible in crossflow direction. Self-evidently, these are only guiding values for a realistic fuel assembly at high Reynolds number. Nevertheless, the estimated forces denote a first insight into a realistic distribution of hydrodynamic forces that are acting upon mixing grids and rods.

Fuel assemblies with different water gap widths

The main objective of this second study is to provide insight into the pressure distribution between realistic fuel assemblies of 17×17 fuel rods with mixing grid. Especially the influence of the water gaps widths on the forces acting on fuel assemblies is addressed. For that reason, three FAs were gathered and arranged in line along the x-axis. Water gaps of different sizes exist between the assemblies. Fig. 8 shows the top view on the three FAs called "A", "B" and "C". As indicated, the water gaps have a width of 2 mm and 5 mm, respectively. At the outer vertical faces of the calculation domain, at a distance of 1 mm from the outer lattice straps of the mixing grids, symmetry boundary conditions are used. This boundary condition simulates a 2 mm water gap. Each fuel assembly is identical to the FA presented in chapter 4.1. A tracer with the concentration 1 is injected from below into assembly "B" in order to mark the flow in the central assembly. Analogously, the tracer with concentration 0 is injected into the assemblies "A" and "C". Fig. 8 shows for the reference mesh, that the tracer injected into assembly "B" penetrates less into assembly "C" (separated by a 5 mm water gap) than into assembly "A" (separated by a 2 mm water gap). This indicates that the wider the water gap, the more the fuel assemblies are separated hydraulically one from another, what reduces inter-assembly exchanges. Test calculations with the fine mesh have confirmed this result. More detail on the pressure distribution within the three assemblies illustrates Fig. 9. Like in Fig. 7a, the pressure is given in a vertical plane at y = 0.1175 m that cuts rod-to-rod gaps in the middle. The whole calculation domain is visualized on the top figure; a section around the mixing grid strap of the central assembly is zoomed in on the lower figure. Upstream of the mixing grid, the highest pressure is located below the three narrow water gaps of 2 mm width, with a maximum value of approximately 12.9 kPa. These peak values exist on one hand below the water gap between assembly A and B and on the other hand below the gaps between the outer symmetry boundary conditions and the adjacent mixing grid. The pressure below the wide gap of 5 mm width is 10.4 kPa and thus 20% lower than that below the narrow gaps. It is interesting to note that the spatial extension of the high-pressure zone below the wide gap is significantly reduced with respect to the high-pressure zone below the narrow gaps. Downstream of the grid, the pressure distribution above each assembly is similar to that given in Fig. 7a for one single assembly.

Velocity vectors were added to the zoomed in section. They show the deflection of the flow by the mixing vanes as already mentioned in Fig. 6a. It is also well visible that the velocity significantly accelerates in the wide water gap, in contrast to the narrow gap, where further wall friction decelerates the flow. In accordance to the Bernoulli equation, the pressure in the wide gap is lower than in the narrow gap. This pressure difference is acting on the external lattice walls of the central mixing grid in direction from the narrow gap to the wide gap. This is shown in more detail in Fig. 10 where the pressure distribution in a horizontal plane at z = 0 m is given. This plane is located at the upper edge of the grid lattice. The figure zooms in the central assembly (B) in order to better visualize the pressure distribution in the water gaps and in the grid lattice. The pressure distribution within the grid was already discussed in Fig. 7b. The difference between the mean pressures in the wide gap and the narrow gap is approximately 3.25 kPa. The pressure differences between the water gaps and within each grid lattice lead to a total mean pressure force of +39 N that is acting in x-direction on the mixing grid of the central fuel assembly. It is important to note, that this pressure force (called Venturi force) has a significant self-stabilizing effect on the alignment of FAs in the reactor core. In fact, as the Venturi force is directed from the narrow to the wide water gap, this force attempts to push the FA in this direction in order to reduce the width of the wide gap. The Venturi force thus tries to even out the widths of all water gaps to similar values. Besides the Venturi forces, local pressure differences within the grid lattice interact with mixing grid and rods (Fig. 10). The resulting pressure forces affecting a mixing grid are summarized in Table 5. Viscous forces are not considered as they are small in cross-wise direction. Three different configurations of the widths of the water gaps are discussed, which confine an assembly in x-direction, namely gaps of 2 mm and 5 mm, 3 mm and 7 mm as well as 2 and 10 mm. In the first two cases, fully developed flow in pure rod bundle is imposed as inlet boundary condition. This hypothesis let in the gap of 10 mm to an unrealistic high axial velocity, since the mixing grids would reduce significantly the velocity in real assemblies. Therefore, a constant velocity of the same mass flow rate is used in this case as inlet boundary condition.

The pressure forces acting on the mixing grids are summarized in Table 5 for all three grids and all three cases. The pressure forces of the calculation with equal width of the water gaps (2 mm) is added for the sake of completeness. It is interesting to note that the axial pressure forces (Fz) are only slightly dependent on the water gap widths (about 200 N in all cases). The pressure force in y-direction (Fy) are consistent among each other and does show only an unimportant dependency on the water gap width. Nevertheless, a doubling of the force from Grid "A" to Grid "C" is observed. It should be noted however that the pressure force in y-direction might be falsified by the used symmetry boundary condition. -32.4 -10.3 192 Concerning the pressure force in x-direction, the situation is much more complex as the water gap widths influence significantly this pressure forces. A significant net pressure force is acting on the central grid B, which increases from -0.5 N for the gap combination 2-2 mm to 39 N in the gap combination 2-5 mm and decreases

B C A

again to 24 N in the gap combination 3-7 mm. Since the inlet boundary condition is different in the gap combination 2-10 mm, this case cannot be compared directly to the other results. It is important to note that comparative calculations with the fine mash were realized for the water gap widths of 2-5 mm and 3-7 mm. These results with the fine meshes have confirmed the vertical and horizontal forces mentioned above.

Fuel assemblies with inclined rods

The main objective of this third study is to provide further information on the pressure distribution in bundles of 17×17 fuel rods; especially the influence of the inclination of the rods on the flow field and the resulting pressure and shear forces is addressed. The geometry of the three assemblies was given in Fig. 3 and is recalled in Fig. 11. The rods in the central assembly are inclined downstream the mixing grid. The water gap between assemblies A and B opens in 470 mm from 3 mm to 7 mm. Accordingly, the gap between assemblies B and C closes from 7 mm to 3 mm. A tracer marks the flow in the central assembly. The pressure forces that act on the rods of the three assemblies A, B and C are summarized in Table 6. The detected difference between assemblies with only upright rods and assemblies, where the central assembly has inclined rods, seems negligible. The difference of the shear forces that are acting on the rods in upright and inclined configuration are also very small. This is shown in Table 7 for the rods of the three assemblies. The calculated differences are too small to draw robust conclusions. The insensitivity of pressure and shear forces on the inclination of the rods might have two reasons. On the one hand, the inclination angle is small. On the other hand, the distribution of both the cross flow velocity and the pressure is very inhomogeneous within the assemblies due to the presence of mixing grids. These inhomogeneous distributions are visualized in Fig. 13 for the central assembly with inclined rods. The fields of the magnitude of the crossflow velocity (+ , + . ) and of pressure are shown in color scale for a plane at z = 0.35 m. A clear direction of the resulting forces cannot be detected.

CONCLUSION

The velocity and pressure distribution within and between full cross-section fuel assemblies with generic but realistic mixing grids have been analyzed by CFD using high performance computing (HPC). The analyzed assembly consists of bundles of 17×17 rods including fuel rods, control rod guide tubes and a mixing grid with split type mixing vanes. The axial length of the analyzed domain spans approximately the distance between two successive mixing grids.

The flow and pressure distribution in three configurations of assemblies with 17×17 rods are compared: a single assembly test case and two multi-assembly test cases, each with three assemblies arranged in line. Intraassembly flow and pressure distribution are analyzed for the single assembly case; inter-assembly flow and pressure distribution are analyzed for the multi-assembly test cases. In the last-mentioned case, two water gaps of different width separate the central assembly from the neighboring fuel assembles, namely 2 and 5 mm as well as 3 and 7 mm. Finally, it was analyzed if inclined rods can significantly affect pressure and shear forces, which are acting on the rods.

For a fuel assembly surrounded by water gaps of 2 mm, the total force (pressure and shear force) acting upon the mixing grid in flow direction is about 230 N. The total force working in cross flow direction on the grid is small. This horizontal force can be falsified by inappropriate boundary conditions as symmetry. The total force acting upon the rods is about 160 N in flow direction and small in crossflow direction. Self-evidently, these are only guiding values for a realistic fuel assembly at high Reynolds number. Nevertheless, the estimated forces denote a first insight into a realistic distribution of hydrodynamic forces that are acting upon mixing grids and rods.

Variations in the width of two oppositely located water gaps significantly influence the pressure force acting on the mixing grid in horizontal direction. Depending on the difference in the width of the water gaps, pressure forces of 39 N and 24 N were calculated for gap couples of 2 mm and 5 mm as well as of 3 mm and 7 mm, respectively. This horizontally acting pressure force (Venturi force) affects the fuel assemblies in the reactor core in a way that the force tries to adjust automatically the width of water gaps to similar values (self-stabilizing effect).

Very precise calculations with 1 billion tetrahedral cells were realized to calculate the hydraulic forces on inclined rods. The predicted forces on inclined rods do not differ significantly from those upon upright rods. The insensitivity of pressure and shear forces on the inclination of the rods might have two reasons: on the one hand, the inclination angle is small and on the other hand, the distribution of both the cross flow velocity and the pressure is very inhomogeneous within the assemblies due to the presence of mixing grids.

It was shown that symmetry boundary conditions could falsify the calculated hydraulic forces acting on a fuel assembly. In order to overcome this problem, the hydraulic forces on an assembly will be analyzed, which is located in the center of nine assemblies. The forces on this central assembly will not be affected by horizontal boundary conditions of the peripheral faces of the assembly cluster. Further, it is planned to vary more systematically the width of the water gaps and to add thermal effects.

Fig. 1 .

 1 Fig.1. Typical fuel assembly (a) and mixing grid of the generic fuel assembly (b).

  Fig.2 Details of the mixing grid: a) Zoom on fuel rods, guide tubes, grid strap, springs and dimples; b) zoom on the split type mixing vanes with triangular central hole

  Fig.3 Configuration of analyzed assemblies: a) single assembly, b) three assemblies in line with upright rods and c) one assembly with inclined rods in the middle between two assemblies with upright rods

Fig. 4

 4 Fig.4 Calculation domain of the fuel assembly with mixing grid, meshed and extrudes zones as well as with inlet and outlet planes.

Fig. 5

 5 Fig.5 Comparison of the speedup of TrioCFD compared to ideal speedup (laminar test case)

  Fig.6 Spatial distribution of velocity magnitude in the 17×17 rods fuel assembly. a) view on the calculation domain, b) Cross-flow velocity close to a guide tube (horizontal plane at z= 0.05 m)

Fig. 7 .

 7 Fig. 7. Spatial distribution of the pressure in cut planes; a) vertical plane at y = 0.1175 m, b) horizontal plane at z=0 m (zoom to the lower left corner of the assembly).

Fig. 8

 8 Fig.8 Group of 3 fuel assemblies; arrangement of the assemblies and visualization of the inter-assembly flow exchanges by a tracer.

Fig. 9 .

 9 Fig. 9. Pressure in a vertical plane with a zoom to the water gaps and mixing grid lattice

Fig. 11

 11 Fig.11 Concentration of a tracer in test case with inclined central rods; arrangement of the assemblies and visualization of the inter-assembly flow exchangesAs soon as the gap widths change, the sharp frontier between high and low tracer concentration becomes diffuse. This was not observed for constant gap widths (Fig.8). The closing gap pushes the flow more profoundly into the central assembly B whereas the opening gap sucks the tracer from assembly B into the gap. These processes are shown in more detail in Fig.12, where the flow interactions across the water gaps are visualized in a horizontal plane at z = 0.35 m; the opening gap is shown in Fig12aand the closing gap is shown in Fig12b. Velocity vectors visualize the direction of the cross flow velocity. It is interesting to note a certain correlation of the flow in the peripheral rods of assembly B: the flow enters and leaves, respectively, the central assembly at the same y-coordinate, and one out of two columns sustains a similar cross flow behavior. It would be interesting to know if, after discharge of the reactor core, the wear of fuel rods close to the water gaps can be related to the different flow fields observed for upright and slightly deformed assemblies.

  Fig.13 Cross flow velocity (a) and pressure (b) in the assembly with inclined rods at z = 0.35 m

Table 1

 1 Summary of the numerical scheme

	Time scheme	Euler backward	1 st order implicit
		Predictor step	GMRES
		Corrector step	CG with SSOR preconditioning
	Spatial discretisation	Convection term	2 nd order upwind
		Diffusion term	2 nd order centred
		Pressure gradient	1 st order
	Turbulence	U-RANS	k-ε with wall functions

Table 2

 2 Strong scaling evaluation of the parallel performance of TrioCFD

	Number of cores	dt TS [s]	Dt CG [ms]	Cells/core	Load imbalance
	6,400	69.6	88	208594	1.003
	12,800	34	46	104297	1.053
	25,600	18.5	25	52148	1.268
	51,200	10.4	15	26074	1.300

Table 3

 3 Characterization of the parallel performance of TrioCFD for assembly calculations

		Reference mesh	Reference mesh	Fine mesh
		One assembly	Three assemblies	Three assemblies
		Occigen	Occigen	Irène-Joliot
	Number of elements	129,094,176	390,881,109	1,036,437,165
	Number of processor cores	2800	8400	16,680
	Number of elements per core	46,105	46533	62,137
	Mean CPU per time step [s]	4.8	11.0	22.6
	Time spent in solvers [%]	69	84	74
	Network traffic per time step [Mb]	0.11×10 6	6.2×10 6	31.0×10 6
	Machine time to convergence [h]	10.6	24.5	50.2

Table 3

 3 Hydraulic forces acting on the mixing grid and rods.

			Mixing grid			Fuel rods and guide tubes
	Force	Fx	Fy	Fz	Fx	Fy	Fz
	[N]	-9	-6	190	-2	-1	0.2
	[N]	-0.005	-0.004	43	0.7	0.3	167

Table 4

 4 Hydraulic forces acting on the mixing grid and the rods of the central assembly.

			Central mixing grid		Central fuel rods and guide tubes
	Force	Fx	Fy	Fz	Fx	Fy	Fz
	[N]	-0.5	-10.6	204	-0.3	-3.3	0.2
	[N]	-0.006	-0.002	40	0.01	-0.03	169

Table 5

 5 Pressure forces acting on the mixing grid.

	Width of the		A-B 2 mm			A-B 2 mm			A-B 3 mm			A-B 2 mm
	Gap		B-C 2 mm			B-C 5 mm			B-C 7 mm			B-C 10 mm
	[N]	Fx	Fy	Fz	Fx	Fy	Fz	Fx	Fy	Fz	Fx	Fy	Fz
	Grid A	-3.2	-6.2	204	-3.1	-5.6	207	18	-5.0	203	-1.6	-4.3	200
	Grid B	-0.5	-10.6	204	39	-8.9	206	24	-9.8	202	30,6	-7.9	192
	Grid C	-4.8	-14.7	204	-42	-13	206	-51	-13.6	203

Table 6

 6 Pressure forces acting on the rod bundles of the three assemblies

	Width of the Gap		Upright rods			Inclined rods
	[N]	Fx	Fy	Fz	Fx	Fy	Fz
	Rods A	6.0	5.2	0.2	5.2	7.1	0.2
	Rods B	1.1	-4.1	0.2	-1.8	-6.1	0.2
	Rods C	-14.3	-9.3	0.2	-9.8	-7.1	0.2

Table 7

 7 Shear forces acting on the rod bundles of the three assemblies

	Width of the Gap		Upright rods			Inclined rods
	* [N]	Fx	Fy	Fz	Fx	Fy	Fz
	Rods A	0.2	0.7	165	0.3	1.2	171
	Rods B	-1.1	-0.02	168	-1.9	0.2	175
	Rods C	-0.3	-0.6	168	0.5	-1.0	175

https://salome-platform.org/

https://www.mcs.anl.gov/petsc

https://wci.llnl.gov/simulation/computer-codes/visit

ACKNOWLEDGEMENT

This work was granted access to the HPC resources of TGCC and CINES under the allocation A0052A07571 and A0072A07571 made by GENCI.

DES-STMF, CEA, Université Paris-Saclay

F-91191 Gif-sur-Yvette, France ulrich.bieder@cea.fr Hydraulic forces on fuel assemblies with different sizes of water gaps are calculated by using the CFD code TrioCFD. Meshes with up to 1 billion cells were used.

Group of 3 fuel assemblies; arrangement of the assemblies and visualization of the inter-assembly flow exchanges by a tracer.

Forces in flow direction and in crossflow direction were quantified. The influence of boundary conditions on the prediction of forces is discussed. The horizontally acting pressure force (Venturi force) affects the fuel assemblies in the reactor core in a way that the force tries to adjust automatically the width of water gaps to similar values (self-stabilizing effect).


Fatima-Zahra Alem

Samira Ananthabothla

Lise Cougnaud

Cedric Affnar

Issmail Nounah

Bakri Youssef

Alessandra Lopes De Oliveira

Samira Ananthabothla Gita

Zoubida Charrouf

Boutayna Rhourri- Frih
email: boutayna.frih@u-bordeaux.fr

Professor Boutayna Rhourri-Frih

Lawsone quantification in Lawsonia inermis L. by HPLC-MS: How does the temperature and pluviometry affect lawsone concentration?

Keywords: Lawsonia inermis L, Lawsone, HPLC-ESI-QTRAP, Mass spectrometry, MRM, Principal, Component analysis, Classification, Correlation, Morocco

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

Introduction

Morocco, country of North Africa, spans from the Mediterranean Sea and Atlantic Ocean. It includes a rich and large nature including mountains areas in the interior body of the kingdom to the Sahara Desert in the south east and the extreme south. This richness of its nature offers a great biodiversity of medicinal and aromatic plants. This biodiversity has opened the field to the development of phytotherapy and aromatherapy in Morocco but also in other countries using these plants.

additionally, this natural richness has contributed to the development of new therapeutic and or preventive strategies with less toxicity.

Mediterranean plants have been used for treatment of several diseases since time immemorial [START_REF] Conway | The Consultation in Phytotherapy E-Book: The Herbal Practitioner's Approach to the Patient[END_REF]. Artemesia herba alba, Nigella sativa, Thymus vulgaris and others were considered as remedies for diabetes mellitus (Eddouks et al., 2002;[START_REF] Kaleem | >Nigella sativa</i> L seeds in diabetic rats[END_REF]. Moreover, different plants were found to have an anti-cancerous activity , such as Olea europae, Salvia officinalis L, Cicer arietinum and Vicia faba that were used for the treatment of lung cancer [START_REF] Alves-Silva | North African Medicinal Plants Traditionally Used in Cancer Therapy[END_REF][START_REF] Chen | Mulberry anthocyanins, cyanidin 3-rutinoside and cyanidin 3-glucoside, exhibited an inhibitory effect on the migration and invasion of a human lung cancer cell line[END_REF][START_REF] Endrini | Anticarcinogenic properties and antioxidant activity of Henna (Lawsonia inermis)[END_REF][START_REF] Marchand | Cancer preventive effects of flavonoids-a review[END_REF].

Among the large list of plants growing in Morocco, Lawsonia inermis L, has drown our attention by its wide spectrum of biological activities. It's a plant widely used in Moroccan traditions especially for cosmetic purpose in all special occasions such as weddings, birth or circumcisions where women get henna on their hair, hands and feet [START_REF] Staats | Ritual, Strategy, or Convention: Social Meanings in the Traditional Women's Baths in Morocco[END_REF].

The origin of this plant is the Babylon region (the Mesopotamia) 3000-1500 before J-C. It was used as a sign of richness in the noble families and especially for the preparation of the bride [START_REF] Aubaile-Sallenave | Les voyages du Henné[END_REF]. Later it was migrated to south Asia especially India and in Egypt where the pharaons used it as a cosmetic agent for dying their feets and hands. The Egyptian Jewish was also using this plant as perfume mixed with oils from other plants.

After the Islam era, the plant was becoming widely used in the Middle East and North Africa. It has several appellations based on the region: in North Africa, it is Henna or Hinna; in India, it is Mehendi; in ancient Egypt, it was pouker, in Hebrew it is kopher, but the common purpose of it usage is its desired dark color used in cosmetic [START_REF] Aubaile-Sallenave | Les voyages du Henné[END_REF]. As previously reported, the responsible compound of the coloring is Lawsone [START_REF] Lee | Hair Dyes Resorcinol and Lawsone Reduce Production of Melanin in Melanoma Cells by Tyrosinase Activity Inhibition and Decreasing Tyrosinase and Microphthalmia-Associated Transcription Factor (MITF) Expression[END_REF]. The concentration of Lawsone in L. inermis L varies from a region to another depending on the conditions of the culture and weather.

Lawsone (Figure 1) is the major compound present in the leaves of L. inermis L responsible for its specific colour used for dying hair, skin and nails [START_REF] Chaudhary | Lawsonia inermis Linnaeus: A Phytopharmacological Review 2, 8. Chemical Composition and Antioxidant Activity of Lawsonia inermis Seed Extracts from Morocco -Issmail Nounah[END_REF]. In addition to the cosmetic use, purified lawsone is also widely used in the coloration of cotton and silk fabric [START_REF] Hasan | Application of Purified Lawsone as Natural Dye on Cotton and Silk Fabric[END_REF]. Phytochemical studies of L. inermis L extract have shown the presence of hundreds of secondary metabolites, of diverse chemical nature, including coumarins, flavonoids, naphthoquinones, quinoids, naphthalene derivatives, triterpenoids, organic acids, xanthones, phenolics, and phenolic glycosides (Ahmed et al., 2000;[START_REF] Badoni Semwal | Lawsonia inermis L. (henna): Ethnobotanical, phytochemical and pharmacological aspects[END_REF][START_REF] Liou | 1,5-Diphenylpent-3-en-1-ynes and methyl naphthalene carboxylates from Lawsonia inermis and their antiinflammatory activity[END_REF][START_REF] Mengoni | SPME-GC-MS analysis of commercial henna samples (Lawsonia inermis L.)[END_REF][START_REF] Singh | Safety and efficacy of hydroalcoholic extract from Lawsonia inermis leaves on lipid profile in alloxan-induced diabetic rats[END_REF]. Among this wide variety of secondary metabolites, several studies focused on naphthoquinone derivatives. This interest stands for the reported anti-microbial effect of Lawone; 2-hydroxy-1,4-naphthoquinone [START_REF] Sauriasari | Cytotoxicity of lawsone and cytoprotective activity of antioxidants in catalase mutant Escherichia coli[END_REF], as well as the leishmanicidal and antiviral agent of Lapachol; 4-hydroxy-3-(3-méthylbut-2-ényl) naphtaléne-1,2dione (ravelo et al., 2003) and finally the antimetastatic melanoma in vitro effect of Plumbagin; 5hydroxy-2-methyl-1,4-naphthoquinone (Alem et al., 2020).

In parallel, its biological activities have been widely described in the literature [START_REF] Badoni Semwal | Lawsonia inermis L. (henna): Ethnobotanical, phytochemical and pharmacological aspects[END_REF][START_REF] Chaudhary | Lawsonia inermis Linnaeus: A Phytopharmacological Review 2, 8. Chemical Composition and Antioxidant Activity of Lawsonia inermis Seed Extracts from Morocco -Issmail Nounah[END_REF][START_REF] Dasgupta | Modulatory effect of Henna leaf (Lawsonia inermis) on drug metabolising phase I and phase II enzymes, antioxidant enzymes, lipid peroxidation and chemically induced skin and forestomach papillomagenesis in mice[END_REF][START_REF] Endrini | Anticarcinogenic properties and antioxidant activity of Henna (Lawsonia inermis)[END_REF][START_REF] Nounah | Chemical Composition and Antioxidant Activity of Lawsonia inermis Seed Extracts from Morocco[END_REF][START_REF] Pour | Lawsonia inermis L. leaves aqueous extract as a natural antioxidant and antibacterial product[END_REF]. Its antioxidant potential has been found to be higher than the Vitamin E and α-tocopherol [START_REF] Endrini | Anticarcinogenic properties and antioxidant activity of Henna (Lawsonia inermis)[END_REF]. The same study showed that L. inermis L extract displayed a cytotoxic effect against liver cancer cell line. Furthermore, this species constitutes a high antibacterial candidate against several pathogens: Staphylococcus aureus, Bacillus cereus, Corynebacterium pseudotuberculosis, Klebsiella pneumonia, Escherichia coli and Salmonella enterica serovar typhi [START_REF] Pour | Lawsonia inermis L. leaves aqueous extract as a natural antioxidant and antibacterial product[END_REF][START_REF] Sauriasari | Cytotoxicity of lawsone and cytoprotective activity of antioxidants in catalase mutant Escherichia coli[END_REF].

The aim of this work is to analyze and quantify lawsone in Moroccan L. inermis L that grows specifically in the south of Morocco and compare its chemical profile with other reported.

Nevertheless, the difficulty in quantifying lawsone is mainly due to the presence of several quinones concurrently in the extract, which leads to an error in estimating the accurate quantity of this compound. Moreover several works had reported this difficulty and tried to estimate the relative error encountered in the quantification of lawsone using LC-UV-VIS [START_REF] Li | Evaluation and correction on quinones' quantification errors: Derived from the coexistence of different quinone species and pH-sensitive feature[END_REF] as an analytical technique or HTPLC [START_REF] Moharana | Comparative Thin-Layer Chromatographic Studies and Development of a High-Performance Thin-Layer Chromatography Method for the Quantification of Lawsone in Natural and Micropropagated Plant Parts of Lawsonia inermis L. JPC[END_REF][START_REF] Patel | Method development for Lawsone estimation in Trichup herbal hair powder by high-performance thin layer chromatography[END_REF]. Another problem in the literature dealt with the quality of the plant after the bleaching operation carried out before use in cosmetics and which can seriously affect the concentration of lawsone [START_REF] Jyotshna | Validated method for quality assessment of henna (Lawsonia inermis L.) leaves after postharvest blanching and its cosmetic application[END_REF].

In this work, a simple and unique specific method is used for the first time to quantify lawsone in Lawsonia inersmis L from four regions renowned to be the leading ones in the production of this plant in Morocco, namely Foum Zguid, Tazarine, Er Rachidia and Zagoura (Figure 2) and to compare them using principal component analysis. The use of tandem mass spectrometry has enabled the development of a specific and sensitive method for the quantification of lawsone in the different samples from the four regions and the correlation of these data with the geographic data that we were able to collect.

Generally, the humidity and temperature as well as the nature of the soil influence the accumulation of secondary metabolites in a plant. In this work, it was demonstrated for the first time that a small change in temperature alone induces a large change in the concentration of lawsone in henna in geographically very close regions.

Material and methods

Chemical and reagents

Analytical grade ethanol and methanol purchased from Sigma Aldrich were used throughout the extraction and analysis process. Lawsone was purchased from Sigma Aldrich Chemical Corp. (Saint-Quentin-Fallavier, France). Ultra-pure watter was obtained with an Elga apparatus.

Plant material and geographical zone of culture

Lawsonia inermis L was kindly provided fresh from certified Moroccan cooperatives: The Tazarine henna was purchased from Ait Matan cooperative, Foum Zguid, Er Rachidia and Zagora samples were purchased from La Batha cooperative. The harvest and drying process were same for all the samples. The choice of geographic areas for growing henna was made based on the popularity of its quality by locals. It is known that these 4 regions are the most producers of the best quality of Henna in Morocco.

Climate data collection

Climate data (temperature and rainfall) were collected from website climate data .org. and distance from ocean was calculated using the web calculator (distance.org -Ressources et information concernant Distance Calculator Resources and Information.). The collected data are showed in table 1 andfigure 3.

Extraction

The henna leaves were freshly harvested, dried and grounded. The ultrasound-assisted extraction was performed at room temperature in an ultrasonic bath device with 46KH3 transducers annealed to the bottom. However, a slight increase of the temperature to reach 37°C was registered due to the ultrasound's frictions. 10 g of dried henna leaves powder were placed in an Erlenmeyer flask (500 mL) and mixed with 100 mL of ethanol and sonicated for 45 minutes. This process was done simultaneously for all the samples to get the same extraction conditions and to avoid the influence of temperature variation. The suspension was recovered and a second extraction was carried out for 45 additional minutes with 100 ml ethanol. Then, the extracts obtained were filtered and evaporated with a rotavapor for 4 hours. Finally, the extracts were lyophilized. The lyophilizate was recovered with mobile phase prior to LC-MS-MS analysis.

HPLC-ESI-QTrap-MS conditions

A Shimadzu gradient LC system (CBM-20A) interfaced to a Qtrap 5500 mass spectrometer equipped with a turboV electrospray source (SCIEX API 5500) was used for LC-ESI-MS/MS analysis. The separation of the compound was carried out on a C18 Phenomenex 100x1mm column with 3 µm particles, 110 Å. The flow rate was 100 µl/min. Henna extract was dissolved in 10/90 (H2O/MethOH). Afterward, 5 µl of sample volume was injected to the HPLC system.

Standard preparation

Standard stock solutions of lawsone were prepared in the mobile phase at the concentration of 1mg/ml, and then several dilutions were done before the analysis: 2.5-5-10-20 µg/mL.

Data Analysis

The analysis was done by extracting the value of area under curve for m/z (175) analyzed in triplicate, then the standard calibration curves were established.

The final content of lawsone was estimated in each sample using the calibration curve equation after validation.

Method validation

The linearity and precision were evaluated by analyzing three samples at each mass concentration.

The precision of the method at each mass concentration was expressed as a coefficient of variation (CV) by calculating the standard deviation of the mean calculated mass concentration.

Statistical analysis:

A one-way analysis of variance (ANOVA) followed by the Tukey post hoc test were used to analyse statistical significance (p < 0.05) using Xlstat software (Addinsoft 2020. XLSTAT statistical and data analysis solution. NEW York, USA). The analysis was carried out at least in triplicate. Principal component analysis (PCA) was performed on lawsone concentration calculated from every individual sample, temperature of the four regions (°C), rainfall parameters (mm) and distance from ocean (km).

Results and discussion

Geographical analysis:

Morocco is characterized by a very marked spatial and temporal variability in precipitation. The northwest is on average more watered than the rest of the Kingdom. Even in this region, the average annual cumulative precipitation varies considerably. For example, it can reach more than 800 mm on the reliefs, while it does not exceed 300 mm on the surrounding plains.

The Moroccan climate is mostly arid. Two exceptions are recorded: the far north where the climate is subhumid, and the Atlas region where the climate is humid. The further we go to the south and south-east of Morocco, the drier the climate is.

This variability in temperature and precipitation is the main source of plant wealth found there. Additionally, we can notice the direct effect that this variability in temperature and humidity can have on the chemical composition, mainly the secondary metabolites of its plant kingdom.

Secondary metabolites are the main source of the pharmaceuticals and additives in all the aromatic and medicinal plants. They also play an important role in the plant defense from biotic and abiotic stress. The biosynthesis of this metabolites is related to the conditions of the growth. Among the influencers on this production there are high temperature, cold, salinity, light, phytohormones [START_REF] Akula | Influence of abiotic stress signals on secondary metabolites in plants[END_REF] The four regions selected for the quantification of lawsone in Lawsonia inermis L, all belong to the central-east region of Morocco (figure 2). Despite this, they all have different mean temperatures and pluviometry; these data were considered to study the influence of the climate on the lawsone content in the different samples.

The climate data for the 4 regions are summarized in Figure 3 andtable 1.

The four regions show very close temperature in the summer season but a very different precipitation during the year. Tazarine is by far the region with the highest precipitation attempting a mean of 70mm in April while Foum Zguid, Er Rachidia and Zagora the mean precipitation value doesn't exceed 20 mm Foum Zguid is the hottest region with a mean value of 35°C (95 °F) in the summer, while Er Rachidia has registered the lowest temperature around 30°C in the summer.

Stability

Before performing the quantitative analyzes in LC-MS-MS, the lawsone stability study was carried out by preparing standard lawsone solutions in an ethanol / water mixture. The solutions were exposed to daylight at room temperature during four days. Meanwhile, daily analysis using the developed method iin HPLC-UV (244nm) indicated that the recoveries of the target peak of lawsone decreased over the time. We noted a decrease in peak area of 25% after 48 hours and 40% after 4 days ( Figure 4). The validation results indicated that the developed method can provide accurate quantification results within 36 hours post-preparation. This stability test allowed us to develop the quantification analysis conditions afterwards with the use of shaded vials and thermostatically controlled autosampler in the HPLC system.

Liquid chromatography optimization

Plant extracts are always very rich in various compounds, including a large number of isomers.

Despite the use of a very specific mass spectrometry method, the development and optimization of a separation method is essential to fill a quantitative analysis.

To optimize the analytical method, we have worked first on HPLC-UV (244 nm) using Lawsone standard solutions, Reverse Phase Chromatography was the preferred mode because of the low polarity of Lawsone.

Initially a column of 50mm was used for the separation, but the shape and retention time of lawsone peak was unsatisfactory (2.5 min), we decided to use a 100 mm column where retention time of lawsone peak was registered at 4.8 minutes.

Following a comparison of the chromatographic efficiency and asymmetry of the lawsone peak, the best compromise between chromatographic resolution, signal intensity, and time of analysis was achieved using a gradient elution composed of MeOH and Water with a gradient slope of 7.8%/min of organic solvent. This program led to narrower peak of Lawsone with higher efficiency than the one that was obtained with a slope of 2.5% (Figure 5).

The gradient elution was set as the following: 0-1min: 70% MeOH and 30% H2O; 1-8 min:100%

MeOH; 8-14 min; 14-15 min: 70% MeOH and 30% H2O and 15-20 min: 70% MeOH and 30% H2O.

Mass spectrometry optimization ions

Due to its high specificity and sensitivity, liquid chromatography coupled to mass spectrometry was used for the quantitative analysis of lawsone using an ESI-QTrap-MS instrument.

To optimize mass spectrometer source parameters, declustering potential, collision energy, and collision cell exit potential (table 1), Lawsone standard solution was injected in the ESI-Qtrap-MS at the concentration of 1. 25 µg/ml.

A simple mass screening, followed by a MS2 analysis was performed to determine different fragments in the positive ionization mode (figure 6).

The multi-reaction monitoring mode (MRM) was also performed and different transition were followed in order to increase the LOD and LOQ.

The signal to noise ratio was calculated experimentally from both chromatogram in MS2 and MRM mode at the lowest concentration from the calibration curve (1. 25 µg/ml), the S/N of the pic at 4. 8 min was found equal to 16 for MS2 and to 23 for MRM (Figure 7).

Calibration curve

The calibration curves were established to quantify lawsone in the four samples of henna leaves harvested from different regions of Morocco.

Three independent calibrators for each concentration level were prepared (1.25-2.5-5-10-25 µg/ml).

The linearity of response of the LC/MS system across the selected concentration range was determined by linear regression analysis in MS2 and MRM modes.

In all cases linear regression analysis of the curves afforded determination coefficients r 2 > 0.96 indicating excellent linearity of the method within the concentration ranges (Supplementary information, graph A). Also, student test on the slope and the origin was validated.

Fisher test on both curves showed a value higher than F0, proving the significance of regression.

Quantification

The LOD and LOQ were determined experimentally using the S/N Ration, we estimated the LOQ to 0. 76 µg/ml in MS2 and 0. 54 µg/ml in MRM.

The calibration curve obtained with the MS2 mode was used for the calculation of the concentration in the four regions. Foum Zguid has shown the highest (0.6%) amount of Lawsone in the extract analyzed (Figure 7). Our results are in accordance with previous study that quantified Lawsone in Henna leaves in Egypt. Six samples were quantified and the range of Lawsone concentration was found between 0.004 and 0.608g/100g of dried material obtained using high performance thin layer chromatography [START_REF] Charoensup | Pharmacognostic Specifications and Lawsone Content of Lawsonia inermis Leaves[END_REF][START_REF] El-Shaer | Determination of lawsone in henna powders by high performance thin layer chromatography[END_REF].

Principal component analysis

The quantitative data were again analyzed using Principal Component Analysis (PCA) with the aim to differentiate the samples from each other and evaluate which variables among temperature precipitation or distance from ocean, influences their location in the 2D space.

The created dataset was composed by 12 observation (4 samples in biological triplicate) and 4 variables (distance from the Ocean, rainfall, mean temperature and content of lawsone). Score Plot obtained is showed in Figure 9, PC1 represents a variability of 70.9% and PC2: 17.04% consequently this PCA shows 87.94% of variability.

As evidenced, samples are separated along first component in two main groups: first group, (in the highest negative values), is composed by samples collected in Errachidia and Tazarine, the second one, (in the highest positive values), by those collected in Zagora and Foum Zguid (Fig. 10 A andB).

In order to evaluate the influence of each variable on the classification of the samples, the loading plot obtained for the same dataset was then studied and results are shown in Fig. 10 C.

As shown, along the highest positive values of PC1 are located the regions with high mean temperature in the space corresponding to the samples collected in Foum Zguid and Zagora that showed. By contrast, the main information able to discriminate samples from Errachidia and Tazarine is due to the variable pluviometry which is located through the direction of the samples.

Quantitative differences in the content of identified lawsone among samples collected in different geographical areas in south east of Morocco was clearly demonstrated and the variable temperature is a critical component that influences the content of lawsone and its concentration in the leaves of Lawsonia inermis L

Conclusion

In this work, the discrimination of lawsone content in Lawsonia inermis L in four regions from South east of Morocco was presented.

A LC-MS/MS targeted quantitative analysis for lawsone estimation on Lawsonia Inermis L samples collected in the different areas. After statistical analyses, important correlation between temperature and lawsone concentration was demonstrated. In particular two regions showed a very high % of lawsone; Foum Zguid and Zagora that were more in the south. In the opposite there was no significant correlation with the precipitation. Samples with very humid climate and high pluviometry showed lower content of Lawsone.

Compound

Figure 1 :

 1 Figure 1: Lawsone two-dimensional chemical structure and Lawsonia inermis L leaves and powder.

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 7 :Figure 8 :Figure 9 :

 2345789 Figure 2: Map situation of the four South Eastern Moroccan regions where Lawsonia inermis L (Henna) was harvested for this study: Tazzarin, Foum Zguid, Errachidia and Zagora.

Figure 10 :

 10 Figure 10: Loading 2D plot of dataset of 12 observations: 4 samples in biological triplicate and 4 variables (Distance from Ocean, rainfall, mean temperature and quantification of lawsone). In (A) is plotted the variables distribution on the bidimensional plot, in (B) are represented the observations distribution along the two main components F1 and F2 and in (C) is plotted the co-distribution of observations (samples) and variables on the PCA biplot

Table 2 :

 2 Optimized values of mass spectrometer ionization source parameters in Multi Reaction

		Molecular	Ionisation mode	Parent ion	Product ion
		weight		(ESI)
	Lawsone	174.15 g/mol		Positive	175	101-77-70
	2-hydroxy-1,4-naphtho-
	quinone
	Monitoring mode for the rhree main transitions: 175	101; 175 77; 175 70; where DP:
	Declustering Potential, CE: Collision Energy, CxP: Collision Cell exit Potential

Extract