Rémi Gonin
email: remi.gonin@toulouse-inp.fr

Pierre Horgue

Romain Guibert

David Fabre

Rémi Bourguet

Fouad Ammouri

Elena Vyazmina

A computational fluid dynamic study of the filling of a gaseous hydrogen tank under two contrasted scenarios

Keywords: compressed hydrogen, tank filling, thermal gradient, CFD 2010 MSC: 00-01, 99-00 Highlights

The filling of a horizontal hydrogen tank designed for light duty vehicles is investigated by means of multi-physics numerical simulations. The simulation approach, implemented in OpenFOAM, includes compressible Reynolds-Averaged Navier-Stokes (RANS) modeling of the fluid flow and heat transfer in the solid parts. The simulations are carried out for 2D-axisymmetric and 3D configurations. Two filling scenarios of the tank, leading to two distinct thermal behaviors, i.e. homogeneous versus heterogeneous, are simulated and compared to the experimental data issued from the HyTransfer project. In the homogeneous case, where no thermal stratification occurs, the 2D and 3D simulation results are close to the experimental ones. A phenomenon of jet flapping is identified via the 3D simulation. In the heterogeneous case, where thermal stratification occurs, the 3D simulation captures an averaged temperature close to the experimental one, as well as the instant at which the thermal gradients appear. It also captures the deflection of the jet, which is a central element in the emergence of the thermal gradients.

• The kω SST turbulence model underestimates the thermal gradients.

Introduction

Context

Hydrogen is currently considered as an alternative source of energy for light and heavy vehicles compared to regular fossil fuel [START_REF] Nouvelot | Performance metrics for refuelling protocols for heavy duty hydrogen vehicules, fCH-04-2-2019:Refuelling Protocols for Medium and Heavy-Duty Vehicules[END_REF]. For practical convenience, the filling time of the hydrogen tank of light duty vehicule should be comparable to that of a standard fuel tank, about 3-5 minutes [START_REF] Council | A sustainable pathway for the global energy transition[END_REF], for a similar operation range, i.e. 500 km. To match these requirements, the internal pressure reaches 70 M P a at a standard temperature. Composite and plastic materials are employed in order to limit tank weight. These materials are characterized by much lower thermal conductivities [START_REF] Simonovski | Thermal simulations of a hydrogen storage tank during fast filling[END_REF] compared to metallic tanks, leading to reduced heat exchanges with the exterior environment and to a sharp rise of the gas and wall temperatures. The J2601 standard of the Society of Automobile Engineers (SAE) organism [START_REF] Committee | Fueling protocols for light duty gaseous hydrogen surface vehicles[END_REF] advises admissible temperature between -40 • C and 85 • C. Inner temperatures need to be controlled following these specifications.

As shown in previous studies [START_REF] Melideo | CFD analysis of fast filling strategies for hydrogen tanks and their effects on key-parameters[END_REF][START_REF] Ortiz Cebolla | Effect of precooled inlet gas temperature and mass flow rate on final state of charge during hydrogen vehicle refueling[END_REF][START_REF] Kang | Charging strategy through analysis of charging influence factors of ultra-light hydrogen storage with polyethylene terephthalate liner[END_REF], reducing the inlet temperature is a possible way to decrease the maximum temperature reached during the filling. As a result, current filling stations precool the inlet flow down to -40 • C to reduce the final gas temperature in the tank. However, this solution presents some limitations. It is energetically costly and, considering that the temperature field may be heterogeneous during filling due to buoyancy effects, high temperatures may be reached along the upper walls even with precooling. This heterogeneous thermal field phenomenon has been highlighted during the HyTransfer project [START_REF] European Commission | Hytransfer[END_REF], funded by the Fuel Cells and Hydrogen Joint Undertaking (FCH JU) partnership. This program involved experiments and computational fluid dynamics (CFD) in two institutions, Air Liquide Advanced Technology (ALAT) and the Joint Research Centre (JRC). With thermal sensors located inside the hydrogen bulk, the thermal behavior of the fluid has been investigated and showed internal thermal gradients for some tank geometries and under some specific filling conditions. However, the limited number of probes within the fluid region did not allow to clarify the mechanism leading to thermal gradients.

Historically, the first attempts to model tank filling were based on thermodynamic zero dimension (0D) models, which consider the fluid bulk as a unique control volume on which mass and energy balance equations are applied (see e.g [START_REF] Reynolds | Analytical investigation of the blowdown and charging processes in a single gas receiver, including the effects of heat transfer[END_REF][START_REF] Muto | An experimental investigation of the effects of heat transfer during charging and blowdown of single Gas receivers[END_REF][START_REF] Bourgeois | Evaluating the temperature inside a tank during a filling with highly-pressurized gas[END_REF][START_REF] Molkov | Physical model of onboard hydrogen storage tank thermal behaviour during fuelling[END_REF]). These models are suitable to obtain rapid estimations of the averaged temperature of the inner fluid. When reducing the injection velocity, thermal gradients appear in the tank and 0D models cannot capture local temperature variations. Terada et al. [START_REF] Terada | Thermal behavior in hydrogen storage tank for fcv on fast filling[END_REF] observed that this transition occurs below a critical injection velocity of 5 m/s . This criterion was confirmed on Hy-Transfer cases [START_REF] Devilliers | D7.5 hytransfer final synthesis of the project findings for the industry public[END_REF] and gives insights into the occurrence of thermal gradients; however the underlying physical phenomena still need to be explored.

With the increase of computer resources, CFD can now be employed to investigate these configurations. CFD gives access to the entire velocity, pressure and temperature fields, which are inaccessible via experiments. The HyTransfer experimental database has served as a benchmark for validation of the simulations in both involved teams. On one hand, within JRC team, simulations were carried out with the commercial CFD software Ansys CFX, using a Redlich-Kwong [START_REF] Redlich | On the thermodynamics of solutions. v. an equation of state. fugacities of gaseous solutions[END_REF] equation of state as real gas equation of hydrogen. For turbulence modeling, in the first study belonging to the HyTransfer project [START_REF] Melideo | CFD simulations of filling and emptying of hydrogen tanks[END_REF], a modified kϵ turbulence model [START_REF] Ouellette | Turbulent transient gas injections[END_REF] previously used in similar studies from the same team [START_REF] Melideo | CFD analysis of fast filling strategies for hydrogen tanks and their effects on key-parameters[END_REF][START_REF] Galassi | CFD analysis of fast filling scenarios for 70 MPa hydrogen type IV tanks[END_REF][START_REF] Galassi | Assessment of CFD models for hydrogen fast filling simulations[END_REF][START_REF] Melideo | CFD model performance benchmark of fast filling simulations of hydrogen tanks with pre-cooling[END_REF] was used. Subsequently, in Refs. [START_REF] Melideo | Effects of the injector direction on the temperature distribution during filling of hydrogen tanks[END_REF][START_REF] Melideo | Effects of some key-parameters on the thermal stratification in hydrogen tanks during the filling process[END_REF] a kω shear stress transport (SST) model with a Γ -Θ model transitional turbulent-laminar model described by Langtry et Menter [START_REF] Langtry | Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes[END_REF] was used to improve the capture of heterogeneities. These simulations captured the onset of thermal stratification but tended to be less reliable to quantify the stratification level. On the other hand, in the ALAT team, simulations were done with the commercial CFD software Ansys Fluent Inc. using a half domain composed of 640 000 cells, a kω SST turbulence model and real gas model data tabulated from REFPROP© v9.1 database. These simulations were able to predict the onset of stratification phenomenon but not the gradient values.

Case selection from experimental data

To avoid exceeding the temperature limit of 85 • C during the filling, the injected gas is cooled which is energetically costly. The HyTransfer project [START_REF] European Commission | Hytransfer[END_REF] aimed at optimizing fast filling of compressed hydrogen tank using an experimental setup, to find new filling protocols, limiting the pre-cooling, while guaranteeing safety. In this project, three tanks have been used: (i) a type III tank, i.e. composite material with metal liner, capacity of 40 L, provided by Dynetek Industries Ltd, (ii) a type IV tank, i.e. composite material with plastic liner, capacity of 37 L, and (iii) a type IV tank capacity of 531 L provided by ®Hexagon. The first two are considered as small tanks with an aspect ratio L/D (tank length over tank inner diameter) about 2.7 and 2.4. The aspect ratio of the third one is about 5.5. All tanks have been placed in horizontal position, the gravity is normal to the tank axis. The main filling parameters are the inlet mass flowrate and the injection diameter.

The present study will only focus on small tanks for light vehicles. Given the reported cases which lead to different results in terms of thermal gradients, the ®Hexagon 37 L appears to be the best candidate to run CFD validation. For each experimental case, the local temperature has been recorded at selected points. There are 10 probes for ALAT and 6 for JRC. Consequently, local information about the thermal field is limited by the number and location of the probes. To compare the different results, a quantification of the thermal gradients is needed. Based on the available data [START_REF] Ravinel | D4.1 report on the experimental filling test campaign[END_REF], the maximum difference between the highest probe temperature T max and the averaged probe temperature T av is proposed as an indicator of heterogeneity at each time. Table 1 presents ALAT results. Results from JRC [START_REF] De Miguel | Influence of the gas injector configuration on the temperature evolution during refueling of on-board hydrogen tanks[END_REF] showed the same trend.

For the present study, the best case and the worst case scenarios have been be selected: (i) a case with quasi uniform thermal field for 3 mm injection diameter and 8 g/s averaged mass flowrate, which will be named the homogeneous case, 1: Thermal gradients for the 37 L ®Hexagon tank from ALAT [START_REF] Ravinel | D4.1 report on the experimental filling test campaign[END_REF]. Selected cases in bold. *Estimated assuming a uniform velocity profile in the pipe and using the density calculated via the measured inlet temperature and pressure in the tank, the measured mass flowrate and pipe cross section.

and (ii) a case with large thermal gradient for 10 mm injection diameter and 2 g/s averaged mass flowrate, which will be named the heterogeneous case.

Note that for the heterogeneous case, the 10 mm injection diameter corresponds to filling directly into the plug aperture without an injection pipe, while for the homogeneous case, an injection pipe of 3 mm of diameter and 100 mm length is entering into the hydrogen bulk. As more probes are available, the ALAT experimental data were chosen for the present calculations. The experimental absolute pressure in the bulk and local temperatures will be used to set the initial and boundary conditions and to validate the numerical method.

Objectives of this work

The main objective of the present study is to capture, numerically, the different thermal behaviors which occur during tank filling. The different injection scenarios involve physical phenomena related to the three-dimensional nature of the problem so three-dimensional (3D) simulations are carried out. Particular attention is paid to the 3D nature of the underlying physical phenomena and to the impact of a 2D-axisymmetric assumption on their prediction. A specificity of the present study is the use of the open-source software OpenFOAM [START_REF]The OpenFOAM foundation[END_REF]. (version 6.0 from CFD Direct Ltd).

Physical system and its modeling

Physical system

During gas injection in a closed volume, the pressure and temperature gradually rise. Figure 1 represents the physical configuration. The tank is laying along its main axis, the x-axis, in a horizontal position and the incoming gas is injected along the x-direction. The gravity (g) is acting vertically along the y-axis. The tank is axisymmetric around the x-axis. It is composed of 4 or 5 solid regions, depending on the presence of an injection pipe: (i) a composite wrap which undergoes the mechanical constraint due to gas pressure, (ii) a plastic liner which limits the molecular diffusion of hydrogen through the wall, (iii) metallic bosses and (iv) plugs at each tank extremities and (v) a metallic injection pipe (for the homogeneous case). The fluid region is the volume enclosed by the solid regions. The tank is surrounded by air at atmospheric pressure and temperature. The simulations are initialized using the pressure and temperature values from the issued experiments. The temperature data come from 10 probes located in the fluid region. Probe locations are given in the HyTransfer project document [START_REF] Ravinel | D4.1 report on the experimental filling test campaign[END_REF] and shown in figure 1. The inlet temperature and tank pressure were also measured and reported in [START_REF] Ravinel | D4.1 report on the experimental filling test campaign[END_REF]. Figure 2, adapted from [START_REF] Ravinel | D4.1 report on the experimental filling test campaign[END_REF], shows plots of all probe measurement for the homogeneous case (plot A) and heterogeneous case (plot B). The results emphasize the contrasted thermal behaviors between the two cases: (i) in the homogeneous case, all probe temperatures remain close to each other for each part; (ii) for the heterogeneous case, the temperatures diverge at t ≈ 100 s .

Model

The simulation is based on a coupled fluid dynamics and heat transfer model, called conjugate heat transfer (CHT) [START_REF] Perelman | On conjugated problems of heat transfer[END_REF] method. Depending if the region is fluid or solid, different sets of models and equations are considered.

In solid regions, the governing equation is the energy conservation equation, in order to calculate temperature diffusion across walls. As experiments have been carried out outdoors without control of the environmental conditions, the temperature was not uniform at the beginning of the filling. Assuming a negligible impact of the initial temperature disparity on the final fluid and liner temperatures, all tank regions were initialized at the averaged temperature of the fluid region. Thermophysical properties are constant and can be found in table 2.

On the boundaries, a perfect thermal contact is assumed between the different solid regions, continuity of the heat flux and temperature on walls in contact between solid regions. Heat transfers between the outer faces and the atmosphere are modeled by Newton law, The plug and the pipe represent a marginal fraction of the solid tank mass and the impact of this error is therefore negligible.

q = hS∆T , (1)
with q [W ] the thermal flux, h

[W/m 2 /K] the heat transfer coefficient, S [m 2 ]
the surface, and ∆T [K] the difference between the solid temperature and the exterior temperature. As set in the JRC work [START_REF] Melideo | CFD simulations of filling and emptying of hydrogen tanks[END_REF][START_REF] Galassi | CFD analysis of fast filling scenarios for 70 MPa hydrogen type IV tanks[END_REF][START_REF] Melideo | CFD model performance benchmark of fast filling simulations of hydrogen tanks with pre-cooling[END_REF], h = 6 W/m 2 /K and the external temperature is set to 18 • C, based on the experimental data values.

In the fluid region, the governing equations are the Navier-Stokes equations, i.e. the mass, momentum and energy conservation equation. The Reynolds number based on the pipe diameter, injection velocity and viscosity ranges from 15 000 to 45 000 for the heterogeneous case and from 50 000 to 600 000 for the homogeneous case. At these Reynolds numbers, the flow is turbulent [START_REF] Eckhardt | Turbulence transition in pipe flow[END_REF] at the pipe outlet. The turbulence is modeled using the Boussinesq assumption with an unsteady Reynolds-averaged Navier-Stokes (URANS) approach and the k-ω SST equations [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. This model was shown to be accurate in a fully developed turbulent region like a jet flow [START_REF] Kmecova | Circular free jets: CFD simulations with various turbulence models and their comparison with theoretical 25 solutions[END_REF][START_REF] Meslem | Comparison of turbulence models in simulating jet flow from a cross-shaped orifice[END_REF] and in a pipe flow [START_REF] Aliferis | Pipe flow simulation for reynolds number up to[END_REF][START_REF] Hedlund | Evaluation of rans turbulence models for the simulation of channel flow[END_REF]. Assuming that Prandtl number is close to 1, i.e. the thermal boundary layer and the velocity boundary layer have similar thicknesses, it is expected to capture the thin thermal layer, with an adequate mesh.

-100 Concerning the boundary conditions, for the turbulent kinetic energy k and the turbulent dissipation rate ω, no wall law is used in the heterogeneous case and a wall law (in OpenFOAM, kLowReWallFunction) [START_REF] Fangqing | A thorough description of how WallFunctions are implemented InOpenFOAM[END_REF] is employed in the homogeneous case. A no slip condition is used on the walls for the velocity. At the inlet, a uniform mass flowrate is imposed. It is deduced from the experimental pressure input as represented in figure 2 (green curve). For the temperature, a continuity of temperature and heat flux is imposed between the different regions. As the perfect gas approximation is not valid in the considered range of pressure and temperature, real gas data are required to determine the thermodynamical properties. In the literature, several gas models have been used and details can be found in the review [START_REF] Bourgeois | The temperature evolution in compressed gas filling processes: A review[END_REF] and [START_REF] Li | Review on the research of hydrogen storage system fast refueling in fuel cell vehicle[END_REF]. Here, the relevant properties of hydrogen are obtained using the National Institute of Standards and Technology (NIST) tabulated data from [START_REF] Lemmon | Nist standard reference database 23: reference fluid thermodynamic and transport properties-refprop[END_REF]: (i) the density [kg/m 3 ], (ii) the thermal conductivity [W/m/K], (iii) the specific thermal capacity [J/kg/K], (iv) the kinematic viscosity [kg/m/s] and (v) the specific enthalpy [J/kg]. These properties are read from tables. For example, figure 3 represents a plot of the density table as a function of pressure and temperature.

-

Numerical Method

Solver

The most suitable solver in OpenFOAM is called chtMultiRegionFoam [38] as it is required, in the fluid region, to simulate a transient, buoyant and turbulent flow, which exchanges heat with walls. The ChtMultiRegionFoam solver can run a transient simulation with two algorithms. The Pressure Implicit with Splitting of Operators (PISO) algorithm [START_REF]The PISO algorithm in OpenFOAM[END_REF][START_REF] Issa | Solution of the implicitly discretised fluid flow equations by operator-splitting[END_REF], solves the coupled momentum and density equations. The thermodynamic properties (such as density, viscosity, conductivity) are computed explicitly using the last known pressure and temperature values. It involves one velocity predictor step and, by default, two pressure corrector steps at each temporal iteration, to calculate the velocity and pressure fields. To maintain the stability and accuracy of the simulation, it is recommended to use a Courant-Friedrichs-Lewy number (CFL), lower than 1.

The PIMPLE algorithm [START_REF]The PIMPLE algorithm in OpenFOAM[END_REF] is named after a contraction of the PISO and the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) [START_REF] Caretto | Two calculation procedures for steady, three-dimensional flows with recirculation[END_REF] algorithms. It consists in executing, at each time step, several PISO iterations with updated thermodynamic properties of the fluid, using the last known solution. With this iterative process, CF L > 1 is allowed, i.e. larger time steps and faster computations.

No significant difference has been observed in the results by using these two algorithms. The decision of using one instead of the other was motivated by the possible reduction in computation time. The PIMPLE algorithm was used for the homogeneous case and the PISO algorithm for the heterogeneous case. Regarding time and space discretization, first order schemes were used to maintain the stability of the simulations. Both cases ran using constrained maximum CFL values, with time step ranging from 10 -4 s to 10 -6 s .

Mesh

The meshing process is based on a 3D CAD model from ®Hexagon provided during the HyTransfer project [START_REF] European Commission | Hytransfer[END_REF]. The mesh was generated with two opensource tools: (i) Salome [START_REF][END_REF], and (ii) cfMesh [START_REF] Cfmesh | [END_REF]. Figure 4A shows a general view of the 3D meshes in the (x,y) plane. Here, focus is placed on the meshing strategy of the fluid region, which is more challenging than the solid regions. The objective is to generate a mesh mainly composed of hexahedral cells, with a refinement adapted to the expected physics. The difference of scales between the inner pipe area and the main bulk require a variable refinement strategy to reach the relevant cell size in each region. A view at the injection region can be found in figure 4B. A typical cell size in the bulk is to 5mm. Considering that a turbulent jet typically expands with an angle of 12 • [45], the refinement follows a conical shape, with a slightly larger angle. Figure 4C visualizes this region. In the boundary layer regions, the mesh is composed of five layers with a growth rate of 1.4. This design leads to y + < 1 on all walls for the heterogeneous case, and a wall law is not required. For the homogeneous case, y + ≈ 1 could not be reached in the the injection region, a y + < 50 was enforced and a wall law is used.

3D simulations have been carried out: (i) on the full tank geometry, or (ii) on half of the domain using as a cutting plane (x,y). The first ones will be named full 3D and the second ones 3D. In addition to the 3D meshes, 2D axisymmetric meshes, named 2D, were also employed. These meshes represent a slice of the domain with one cell depth. They have the same refinement property as the 3D meshes. The number of cells for each mesh is reported in table 3.

Validation

Validation of the solver and mesh was performed on the 2D cases in order to reduce computational costs. The homogeneous configuration was chosen for mesh convergence study, due to its higher velocities. Even though gravity is not taken into account in the 2D approach, the thermal gradients are limited in this configuration and the 2D results can be compared with the experimental data. Three cases were simulated, using a coarse mesh (9 193 cells), a medium mesh (15 605 cells) and a fine mesh (34 378 cells). Less than 1 • C of difference in the averaged temperature was noted over the integration time. The medium mesh was selected to design the 3D mesh. Time step sensitivity was studies on this medium mesh by dividing the time step by a factor of 2. This leads to less than 1 • C of difference in the averaged temperature.

Full 3D case and 3D case were compared over 15 seconds in the homogeneous case; the differences in temperature are negligible. In order to limit the computation cost, complete simulations were performed on the 3D cases, i.e. half of the tank. This strategy was also used by the JRC [START_REF] Melideo | CFD simulations of filling and emptying of hydrogen tanks[END_REF] and ALAT [START_REF] Zaepffel | CFD analysis of the different flow regimes occurring during the filling of a hydrogen vehicle tank[END_REF].

Results and discussion

Velocity and temperature fields

Figures 5, 6, 7, 8 show snapshots of the velocity magnitude (upper panel) and the temperature field (lower panel) at three different times over a filling scenario, for the 2D and 3D homogeneous cases, and the 2D and 3D heterogeneous cases.

A x-axis symmetry is applied to the 2D cases to ease visualization. For the homogeneous case, the velocity colorbar is scaled between 0 and 10 m/s to highlight the jet penetration and the velocity decay over time. For heterogeneous case, the velocity colorbar is scaled between 0 and 1 m/s, for the same purpose.

Complementary videos of 3D cases are available in the Supplemental Data.

For each case: the maximum magnitude of velocity is at the top left corner, an averaged value of the thermal field is indicated at the bottom left corner and an isoline of the averaged value of the thermal field is superimposed on the thermal field (bottom panel). Supplementary video related to this article can be found at https://doi.org/10.1016/j.ijhydene.2022.03.26

A B C

In the 2D simulations (figures 5 and 7), due to the axisymmetry assumption, gravity is neglected. High temperature spots appear in the regions which are not reached by the cold jet, like the areas above and below the pipe, close to the liner and left boss wall. The jet remains parallel to the x-axis.

For the 3D homogeneous case (figure 6), the jet is flapping up and down in the (x,y) plane. In the video, it can be seen that the maximum temperature is not localized in the upper part but is oscillating around the injection area, with the same frequency as the jet flapping frequency. The thermal field is homogeneous, i.e. not vertically stratified.

For the 3D heterogeneous case (figure 8), before 50 s the jet is moving upward, as in the initial phase of flapping, but then it starts deflecting toward the lower part. A vertical gradient of temperatures appears, which is depicted by the horizontal isoline of temperature (lower panels at t = 310 s and t = 610 s in figure 8).

Global quantities

In the experiments, temperatures are sampled by a limited number of probes. An arithmetic average on probes may not be a representative value for an averaged temperature of the entire thermal field. Using the CFD, it is possible to average directly the temperature field

T av-f ield = i T i v i V tank , (2)
by weighting each cell temperature T i by the cell volume v i over the total volume V tank , or to average the 10 temperature values T probe,i where probes are experimentally located

T av-probes = 10 i T probe,i 10
.

(3)

Due to its higher thermal disparity, the heterogeneous case is more relevant to compare over these two averaging approaches. The difference is lower than 1.5 • C before t = 100 s and lower than 0.5 • C after t = 100 s. Therefore, the averaged temperature from the experimental probes will be regarded as the reference value for further comparisons. The same question arises for the absolute pressure. CFD results show that the pressure field is close to uniform in the tank. Therefore, a single measurement point in the tank is already relevant to characterize the bulk pressure.

Figure 9 presents a comparison of the averaged temperature and pressure between the 2D, 3D and experimental cases, for the homogeneous case A) and for the heterogeneous case B). Figure 10 presents the heat flux q (4) which can be calculated by

q = i λ ef f ∆T s i , (4)
A) B) C) where s i is the cell surface at wall, ∆T is the difference of temperature between the one at wall and the one at the first cell center and λ ef f [W/m/K] is the effective thermal conductivity defined as

λ ef f = λ + λ t , (5)
where λ is the thermal conductivity and λ t is the turbulent thermal conductivity,

λ t = c p µ t P r t , (6)
with the turbulent Prandtl number P r t = 0.85; c p [J/kg/K] is the specify heat capacity and µ t [kg/m/K] the turbulent dynamic viscosity. Figure 10 displays the averaged surface temperature at fluid-solid interfaces using the following expression to weight the average using a local surface:

T wall = i T i s i S | f luid-solid wall . (7)
In the above expression, T i is the cell temperature at wall, s i is the cell surface at wall and S is the total surface of the wall. Concerning the homogeneous case, figure 9A shows that the averaged temperature issued from the 3D field is within ±1 • C from the experimental value.

In the 2D case, the averaged temperature is higher and farther from the experimental value, compared to the 3D case. This difference of averaged temperature is induced by a difference of heat transfer. Figure 10A shows that the global heat flux (dotted line) is slightly larger for the 3D case than for the 2D case.

This small difference is due to less heat flux through bosses and plugs for the 2D case. The larger flux through the liner does not counterbalance this trend. For the 2D case the gravity is neglected and therefore no natural convection is possible. In areas where heat transfer is led by natural convection, like plug cavity, the heat flux is underestimated. In addition, the 2D assumption constrains the jet on the x-axis, while it is free to flap in the (x,y) plane in the 3D case. This may result in a better mixing of the bulk and avoid some stagnant hot or cold spots. As shown in figure 10C, the deviations of heat flux is associated with large deviations in wall temperature predictions between the 2D and 3D simulations.

Concerning the heterogeneous case, the probe averaged temperature issued from the 3D simulations is within ± 3.5 • C from the experimental value before t = 300 s , and ± 1 • C after, as showing in figure 9B. The 2D case shows again a higher averaged temperature. In the heterogeneous case, gravity plays a role even more important than in the homogeneous case. As can be seen in figure 8, due to buoyancy forces, the jet tends to plunge towards the lower part of the tank. In the 2D axisymmetric case, this phenomenon and the vertical thermal stratification cannot be captured. This leads to a non physical thermal disparity and to an underestimation of the velocity close to the bottom wall (figure 7). A consequence is a significant difference of heat transfer on every walls, between the 2D and 3D cases, as highlighted in figure 10B. As shown in figure 10C, this leads to lower temperatures on solid walls, even though the averaged temperature in the fluid region is higher.

Local quantities

Figures 11 and 12 present a comparison of probe temperatures issued from the 3D simulations and the experiments, for the homogeneous and the heterogeneous cases, respectively. In the homogeneous case (figure 11), the simulation globally matches the experimental values. As shown in figure 6 no thermal stratification occurs and the probe values remain close to the averaged temperature value.

In the heterogeneous case [START_REF] Molkov | Physical model of onboard hydrogen storage tank thermal behaviour during fuelling[END_REF], the simulation results are close to the experimental values until t = 150 s, when the thermal gradient starts to develop. Then the simulation tends to underestimate the thermal gradient: in the upper part of the tank (TT760,TT764,TT767), the simulation values are lower than the experimental values and in the lower part of the tank (TT761,TT765,TT768), the simulation values are higher than experimental values. Turbulence is playing a major role in the process of thermal diffusion and its modeling may explain the differences in the heterogeneous case simulation. The kω SST model predicts high levels of turbulent viscosity in regions supposedly laminar. This may result in higher effective thermal diffusion and then in lower gradients of temperature.

Homogeneous cases with lower thermal gradients are less impacted by this phenomenon. Turbulence models capturing the laminar-turbulent transition may be a solution. The kω SST Scale-Adaptive Simulation (SAS) [START_REF] Menter | A scale-adaptive simulation model for turbulent flow predictions[END_REF] has been mentioned by ALAT [START_REF] Zaepffel | CFD analysis of the different flow regimes occurring during the filling of a hydrogen vehicle tank[END_REF][START_REF] Bourgeois | The temperature evolution in compressed gas filling processes: A review[END_REF] or the Γ-Θ model [START_REF] Langtry | Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes[END_REF] has been used by the JRC [START_REF] Melideo | Effects of the injector direction on the temperature distribution during filling of hydrogen tanks[END_REF] in the latter study. The emergence of thermal gradients was captured but not their precise magnitudes. The impact of such advanced models remains to be investigated.

Conclusion

Two filling scenarios of a horizontal gaseous tank, leading to two different thermal behaviors, a homogeneous case and a heterogeneous case, have been simulated and compared to experimental data issued from the HyTransfer project [START_REF] European Commission | Hytransfer[END_REF]. The simulations have been carried out in 2D using an axisymmetry assumption and in 3D by considering half of the tank.

In the homogeneous case, where no thermal stratification occurs, the 2D and 3D simulation averaged temperature are close to the experimental measurements. Due to the restrictions imposed to the fluid under 2D assumption, the inner jet cannot cool areas close to the injector as in the 3D case. This results in a hot spot close to the injector, 10 • C above the averaged temperature at 150 s. A phenomenon of jet flapping is highlighted by the 3D simulations. This phenomenon seems responsible for the periodic motion of the hot spot close to the injector, where the temperature is only 2 • C above the averaged temperature at 150 s. Future work should focus on clarifying jet oscillation mechanism and its link to heat transfer.

In the heterogeneous case, where thermal stratification occurs, the 2D simulation does not capture the vertical heterogeneity and the thermal gradient onset cannot be detected. The 3D simulation predicts an averaged temperature close to the experimental one. It captures the deflection of the jet, which is one of the key processes in the emergence of the thermal gradients. The instant at which the thermal gradients appear is well estimated. The thermal gradients are however underestimated. At the end of the filling process (630 s), the differences of maximum temperature between the simulation results and the experimental measurements at probe TT764 and TT761 are +15 • C and -8 • C, respectively. It is suspected to be a consequence of the kω SST turbulence model which tends to over predict the turbulence level and the associated thermal diffusion.

Advanced models like the kω SST SAS model, should be investigated to improve the prediction of thermal gradients. From a phenomenological perspective, CFD provides insights into the flow physics and thermal behavior inside the tank. Such insights are expected to impact the development of 0D models and help optimizing tank filling strategy.

Figure 1 :

 1 Figure 1: Hexagon 37 L tank with a 3 mm injection cut through the (x,y) plane. Probe locations are indicated by circles: probes located in the (x,y) plane in white and probes located outside the (x,y) plane in black.

Figure 2 :

 2 Figure 2: Experimental values from hyTransfer [24]: A) the homogeneous case, B) the heterogeneous case.

3 ]Figure 3 :

 33 Figure3: Hydrogen density as a function of pressure and temperature, issued from NIST[START_REF] Lemmon | Nist standard reference database 23: reference fluid thermodynamic and transport properties-refprop[END_REF].

Figure 4 :

 4 Figure 4: Computational meshes: A) a global overview of all regions in the (x,y) plane for the homogeneous case, B) slice at the inlet plane for the heterogeneous case, C) zoom around the injector for the heterogenous case.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 5678 Figure 5: Homogeneous 2D case. For each time instant: A) t = 10 s, B) t = 100 s and C) t = 150 s. The upper panel represents the velocity magnitude field, the lower panel represents the temperature field. The black isoline corresponds to averaged temperature. Ax-axis symmetry is applied to ease visualization.

  2D simulation T av-f ield 3D simulation T av-f ield 3D simulation T av-probes Experimental T av-probes 2D simulation p av-f ield 3D simulation p av-f ield Experimental pressure ±1 [ • C] Experimental T av-probes

Figure 9 :

 9 Figure 9: Comparison between the averaged values of the temperature (left axis) and pressure (right axis) issued from the present simulations and the experiment, in the fluid region, for the A) homogeneous and B) heterogeneous cases.

Figure 10 :

 10 Figure 10: A), B) wall heat flux and C), D) temperature at fluid-solid interfaces as functions of time, for the A), C) homogeneous and B), D) heterogeneous cases.

Figure 11 :

 11 Figure 11: Comparison between the local temperature values issued from the 3D simulation and the experiment, in the homogeneous case. Probe locations are indicated in the upper panel.

Figure 12 :

 12 Figure 12: Comparison between the local temperature values issued from the 3D simulation and the experiment, in the heterogeneous case. Probe locations are indicated in the upper panel.

Table 2 :

 2 Thermophysical properties of the solid regions.

	Region	Material	Thermal conductivity [W/m/K]	Specific heat capacity [J/kg/K]	Density [kg/m 3 ]
	Composite Carbon Fiber	0.616	783	1365
	Liner	Polymer	0.414	1799	914
	Boss	Aluminium 6061	167	896	2700
	Plug	Stainless Steel 316L* 16.2	120	7990
	Pipe	Stainless Steel 316L* 16.2	120	7990

*An erroneous value was used for the specific heat capacity of the stainless steel 316L, which should be close to 500 [J/kg/K].

Table 3 :

 3 Number of cells for each case, in the fluid and solid regions.

Acknowledgment

This study was carried on as part of a PhD work funded by Air Liquide Innovation Campus Paris and the National Association of Recherche and Technology (ANRT). It was performed using HPC resources from CALMIP (grants 2019 and 2020, p19025)


Polycarbonate membrane

Gold film

Counter electrode Potentiostat ! Ag

• Classical TEM imaging does not allow distinguishing between Co and Cu layers • EFTEM images at the Co and Cu energy edges reveal the position of the different layers.

• Variations of diameter and thickness along the wire • Simulations on 20 layers show that : • Due to reduced magnetic constants as compared to bulk, electrodeposited Co layers show multiple remanent states as a function of aspect ratios.

•
• Inter layer dipolar coupling favors antiparallel coupling of the layers but also allows transition between uniform and vortex states.

 For spin torque oscillators : possibility to obtain two different working frequencies at remanence by switching the direction of the saturation field.

• Drawback : the control of the electrodeposition conditions must be improved in order to use it in a reproducible way! q k (°)

Conclusion

Thank you !

• No statistical analysis but seem that H s has no effect.

• For diameters ≤ 70 nm, the dipolar coupling avoid vortex states and uniform sates (OP) aligned with the wire axis become the most stable ones for both direction of Hs. • For large diameter (≈ 90 nm) vortex states with the core perpendicular to the wire axis (VIP) become the most stable for both direction of Hs

Geometry used in simulations

Layers like nanocylinders with titlted surfaces.



  Averaged tilt angle ≈ 20°from wire axis and up to 45°C hemical analysis by Energy Filtered Transmission Electron Microscopy (EFTEM) EFTEM @ Cu-L edge (931eV) EFTEM @ Co-L edge (778eV) Superposition of (b) and (c) HRTEM FFT High Resolution TEM (HREM) • random distribution of Co and Cu polycrystals. • grain mean size of 7 ± 2 nm • fcc phases for Co and Cu (as expected) Structural analysis by TEM methods Electrostatic part ϕ elec (x) = C E ∫V(x,y,z) dz Magnetic part 𝜙 𝑠 (𝑥) = 𝜙 𝑒𝑙𝑒𝑐 (𝑥) + 𝜙 𝑀𝑎𝑔 (𝑥) • Integration of 𝑩 along the electron path  Sensitive to in plane components of 𝑩  2D projection of the 3D magnetic state Incident beam After the sample Ψ 0 = exp(ik.r)Ψ(r)= A S (r) exp i((k.r) + ϕ S (r))The information on magnetization is in ϕ s (r) Need for micromagnetic simulations and phase shift reconstruction to recover the full 3D magnetic state𝜙 𝑀𝑎𝑔 = -𝑒 ℏ න 𝐵 ⊥ 𝑥, 𝑧 𝑑𝑥𝑑𝑧• OOMMF simulations of magnetic remanent states after saturation // or ⊥ to the wire axis•Calculation of magnetic induction 𝐵 and reconstruction of electron phase shift 𝜙 𝑀𝑎𝑔 image : 𝐵 = µ 0 (𝑀+ 𝐻 𝐷 ) => 𝜙 𝑀𝑎𝑔 = -𝒆 ℏ  𝑩 ⊥ x, z . dx. dz Difficulty : high number of geometrical and magnetic inputs ! -Geometry extrapolated from EFTEM : layers tilted around Y and X axis (more complex geometry than a simple disk!) -Adjustable magnetic parameters : • µ 0 Ms = 1.5 T (obtained from other measurements) • exchange value A ≈ 22.10 -12 J/m (bulk = 30 . 10 -12 J/m ) • uniaxial crystal anisotropy field H k because of polycrystals (value K u , direction θ Hk ) with K u between 100 and 200 .10 3 J/m 3 (bulk = 400 .10 3 J/m 3 )  Values lower than bulk because of Cu impurities in Co Z Y X Cell size 2 x 2 x 2 nm 3 From 1 to 20 layers of Co Calculation of magnetic states and electron phase shift reconstruction Use of CALMIP parallel computing center of Toulouse, France • No effect of saturation field H S on the remanent states of multilayers • ≈ 60% of layers in antiparallel (AP) states and ≈ 40% in vortex states (VOP).

•

  et al., Nano Lett. 16, 1230, 2016. AP VOP Co 25 nm /Cu 15 nm : effect of dipolar interlayer coupling ϕ 64 +/-2 nm • 70% of the layers show AP configuration whatever H s : favored by stronger dipolar coupling • 20% of layers switch between IP and VOP states when rotating H S (only possible if H k close from the normal of the layers). • 10% of out of plane state (OP) whatever H S : impossible without dipolar coupling ϕ 64 No effect of Hs and 95% of vortex states (VOP state, with access to chirality for tilted layers) • Simulations on multilayers show that VIP states disappear for Co 50 nm /Cu 50 nm because of dipolar coupling ϕ 84 +/-3 nmPhase diagram of single layer MultilayerN. Biziere et al., J. Appl. Phys. 126, 163906 (2019) Co 50 nm /Cu 50 nm : increased aspect ratio• No effect of Hs and 95% of vortex states (VOP state, with access to chirality for tilted layers)• Simulations on multilayers show that VIP states disappear for Co 50 nm /Cu 50 nm because of dipolar coupling • VIP states observable for very thick Co layers with diameter ≥ 80 nm et al.,J. Appl. Phys. 126, 163906 (2019) • Preliminary results : quantitative agreement but calculated induction overestimated (10-20%) • Interest : removing geometrical uncertainties allow focusing on magnetic parameters of individual layers Simulated B field with OOMMF Experimental B field recovered from 3D tomography Implementation in OOMMF of the layer shape Electrostatic potential of Co 25 nm /Cu 15 nm nanowire (spatial resolution 7 nm) Toward a 3D reconstruction • Principle of 3D holographic vector field electron tomography : numerical reconstruction of 3D geometry and B field components of a nano-object via multiple tilted holograms D. Wolf et al., Communications Physics 2, 87, 2019.

-

  Exchange value A = 22.10 -12 J/m (bulk = 30 . 10 -12 J/m ) -K u =130.10 3 J/m 3 Values lower bulk because of Cu impurities in Co Phase diagrams of a single Co layer • 4 different states can be expected depending on aspect ratio • Strong impact of H k and H s => possibility to get information on H k by changing H S

Co 25 nm /Cu 45 nm : system of reference VOP Multilayer Phase diagram of single layer

• Anti parallel states are stable whatever H S if θ k ≥ 20°• vortex states are stable whatever H S only if θ k ≤ 20° significant of a preferential orientation of H k toward the wire axis ? N. Biziere et al., J. Appl. Phys. 126, 163906 (2019) ϕ 64 +/-2 nm ϕ 70 +/-2 nm • 70% of the layers show AP configuration whatever H s : favored by stronger dipolar coupling • 20% of layers switch between AP and VOP states when rotating H S (only possible if H k close from the normal of the layers).


S Mantion

N Biziere

Magnonic crystals are micro/nano-structured magnetic systems that can be used for passive microwave devices such as filters, waveguides, sensors etc.. The Heusler-based Co 2 MnSi (CMS) alloy is considered as a promising material for magnonic applications due to its relatively high saturation magnetization and low Gilbert damping. Here, we demonstrated with micromagnetic simulations that its cubic crystal anisotropy allows microwave reconfigurable operations and filtering properties at zero bias field in a simple model consisting of a squared antidot magnonic crystal. We showed that a low amplitude (few mT) and ultra-short (≈ 1 ns) initialization field can be used to recover different quasi-uniform remanent states that can couple or not to a microwave field. A first experimental demonstration of the reconfigurable operation at zero magnetic field has been observed by ferromagnetic resonance on a squared antidot lattice with square size 200 nm and spacing 600nm.

Spin waves and Magnonic crystals

E-beam lithography

Ion beam etching

Removing of the resist mask

• Relatively high saturation magnetization (𝜇 0 𝑀 𝑠 ≈ 1.3𝑇) [2] • Relatively high magnetocrystalline anisotropy (𝑲 𝒄 ≈ 𝟏𝟎. 𝟏𝟎 𝟑 to 𝟏𝟕. 𝟏𝟎 𝟑 𝑱. 𝒎 -𝟑 ) [2] • Very low Gilbert damping (2.10 -3 down to 7.10 -4 )

• High Curie temprature (𝑇 𝑐 ≈ 860𝐾)

•

  dynamic states in a CMS based magnonic crystal: A numerical study First experimental measurement of the reconfigurable operation (to be submitted)

• 3 Frequency

 3 Spin polarization above 90% Allow to work at high frequencies Possible reconfigurable properties at zero bias field Long propagation distances for spin waves Process with semiconductor facilities Magnon spintronics applications • Micromagnetic simulations performed with Mumax3 • Model system: a squared antidot lattice in a Co 2 MnSi thin film of thickness 20nm Magnetic parameters: 𝑴 𝒔 = 𝟏. 𝟐𝟔 𝒎𝑻, 𝜸 = 𝟐𝟖. 𝟕 𝑮𝑯𝒛. 𝑻 -𝟏 , 𝑲 𝒄 = 𝟏𝟕. 𝟏𝟎 𝟑 𝑱. 𝒎 -𝟑 , 𝑨 = 𝟏𝟖. 𝟏𝟎 -𝟏𝟐 𝑱. 𝒎 -𝟏 [2] Case 1 : CMS with cubic anisotropy On/off switching of SW modes! Strong cubic anisotropy such as with Co 2 MnSi thin film is required for simple reconfiguration operation ([3] S. Mantion et al., Phys. Rev. Applied 17, 044054 (2022)) Heusler-based Co 2 MnSi thin films Magnetization equilibrium position determined by 𝐻 𝑒𝑓𝑓 Magnonic crystal (MC): magnetic material with periodically modulated properties Creation of frequency band gaps in SW dispersion (analog to Bragg diffusion) Perspectives: • Fabrication of a similar device with different s and d in order to compensate the lower 𝐾 𝑐 obtained • Measurements of the filtering applications at remanence with Magneto Surface Spin Waves (MSSW) • Measurements of FMR and propagating SW modes with varying electrical currents through the MC Electrical contacts Micro antenna Spin Wave transduction [4] Ferromagnetic material Growth via ultra high vacuum magnetron sputtering on MgO [110] substrate Extinction of most of the SW modes but attenuations lower than expected  lower 𝑲 𝒄 induces domain wall at remanence Frequency shifts and/or attenuations of particular SW modes Nonlinear effects Reconfigurable dynamic states at remanence Nanofabrication of the magnonic crystal 𝐻 = 0 𝑚𝑇 Objective: Find quasi-uniform magnetic states at remanence ➔ How we proceeded: Calculation of remanent magnetic states vs initialization field 𝐻 0 (saturation at 1T → 0T) Fix s/d ratio Application of a uniform RF pulse 𝒉 𝑹𝑭 along x Calculation SW mode spectrum: 𝑭𝑭𝑻(𝑴 𝒙,𝒚,𝒛 𝒕, 𝒙, 𝒚, 𝒛 ) temporal sinc pulse Result: Optimized reconfigurable states obtained for s/d ≤ 1/band gaps generation at remanence for Magneto Surface Spin Waves (MSSW) in cubic crystal Co 2 MnSi Filtering properties at remanence in a CMS based magnonic crystal [3] SW dispersion in the case of an infinite MC Application of a localized microwave excitation ℎ 𝑅𝐹 to excite propagating SW towards the MC Quasi-Uniform remanent states: • Switching times below 1.2ns • Pulse amplitudes lower than 6 mT • Decrease of s → decrease of amplitude field pulse Remanent state 𝑀 𝑠 = 1.28 𝑇 𝛾 = 28.8 𝐺𝐻𝑧. 𝑇 -1 𝑲 𝒄 = 𝟏𝟎. 𝟏𝟎 𝟑 𝑱. 𝒎 -𝟑 𝐴 = 18. 10 -12 𝐽. 𝑚 -1

Acknowledgments: These works were performed using HPC resources from CALMIP (Grant 2021-p1554) [1] A V Chumak et al, J. Phys. D: Appl. Phys. 50 (2017) 244001