to accelerate convergence to steady state. Dual time stepping is employed for unsteady simulations. Several pre-conditioning methods for low Mach number flows have been implemented, and for incompressible flows both the Artificial Compressibility method or a fully colocated SIMPLE scheme can be used. Many well validated turbulence models (0 to 5 equations) have been implemented, as well as Hybrid RANS-LES, LES models and a transition model. NSMB is also employed for hypersonic applications, and includes several levels of chemistry modeling for air, CO2 and N2, more details are given in this paper. To simplify the grid generation for complex geometries NSMB uses the so-called fully patch grid approach and it is not needed that grid lines are continous across block interfaces. A futher simplification of the grid generation for complex geometries and in particular for moving components is the chimera approach discussed in more detail in this paper. The next sections briefly summarize recent developments of NSMB.

III. The Chimera method

One of the important elements in CFD simulations is the spatial discretization. This is not a problem for a simple and static geometry, but it becomes a major issue when the geometry is complex and contains moving parts, and this is especially the case when using structured grids. In this case, the grid generation is the limiting factor for the simulations. To overcome this problem, several solutions have been proposed like the Immersed Boundary Method (IBM), Volume Of Fluid (VOF) or Chimera methods. The last one, chimera method (also called overset grid method) appears increasingly for complex research applications and more and more for large industrial applications.

The Chimera method was initially proposed by Benek et al. [START_REF] Benek | A flexible grid embedding technique with application to the euler equations[END_REF] and the method was subsequently improved and employed for many cases [START_REF] Henshaw | An adaptive numerical scheme for high-speed reactive flow on overlapping grids[END_REF][START_REF] Henshaw | Moving overlapping grids with adaptive mesh refinement for high-speed reactive and non-reactive flow[END_REF][START_REF] Ak Kapila | A study of detonation diffraction in the ignition-and-growth model[END_REF][START_REF] Pg Buning | Numerical simulation of the integrated space shuttle vehicle in ascent[END_REF][START_REF] Hinatsu | Numerical computation of unsteady incompressible flow in complex geometry using a composite multigrid technique[END_REF][START_REF] Meakin | Moving body overset grid methods for complete aircraft tiltrotor simulations[END_REF][START_REF] Pearce | Development of a large scale Chimera grid system for the space shuttle launch vehicle[END_REF][START_REF] Maple | A new approach to domain decomposition, the beggar code[END_REF][START_REF] Jespersen | Recent enhancements to OVERFLOW[END_REF][START_REF] Meakin | Chapter 11 : Composite overset structured grids[END_REF][START_REF] Meakin | Unsteady aerodynamic simulation of multiple bodies in relative motion[END_REF][START_REF] Henshaw | Parallel computation of three-dimensional flows using overlapping grids with adaptive mesh refinement[END_REF] ). Its principle is to decompose the complex full domain into simple subdomains independently meshed by curvilinear grids. The only constraint is a superposition of grids. Besides simplifying the mesh generation, this technique offers also a powerful solution to manage moving bodies.

The chimera processus has been summarized for structured grids by Meakin. [START_REF] Meakin | Chapter 11 : Composite overset structured grids[END_REF] The process of chimera method can be divided into four main steps like described by Landmann [START_REF] Landmann | A highly automated parallel Chimera method for overset grids based on the implicit hole cutting technique[END_REF] (fig. 1):

• 1 st step : Detection of overlapped cells. This process consists in finding which cells are overlapped by other cells. The test is based only on the coordinates of the cells.

• 2 nd step : Determination of the overlapped cell status. Three kinds of overlapped cells exists [START_REF] Landmann | A highly automated parallel Chimera method for overset grids based on the implicit hole cutting technique[END_REF] : a) Interpolated cells : the flow data are obtained by interpolation from other overlapping grids.

b) Calculated cells : the values come from the solution of the equations governing the flow. c) Hole cells (or blanked/masked cells) : the values of this kind of cells are never used in the discretisation schemes or overlapped boundary conditions and they are not interpolated nor calculated.

• 3 rd step : Calculation of the interpolations parameters. In this step, the donor cells and the associated weights that compose the interpolation are evaluated for each valid overlapped cell.

• 4 th step : Calculation of the interpolation values.

A three-dimensional structured multi-block chimera method has been implemented in NSMB in order to simulate initially the free falling sphere in a tube. The developped chimera method is fully automatic with few, if not none, of user input and pre-processing and described in the following sub-sections.

III.A. Detection of overlapped cells

The search of overlapped cells is based on a test of inclusion on the coordinates. In order to accelerate the search, the mapping of coordinates in a virtual uniform Cartesian grid is performed. This method is well known as the bucket method . [START_REF] Meakin | A new method for establishing intergrid communication among systems of overset grids[END_REF][START_REF] Samet | The design and analysis of spatial data structures[END_REF][START_REF] Timo | Numerical techniques for complex aeronautical flows[END_REF] An inverse mapping algorithm creates an index array that link virtual grid to the real coordinates. The search begins on the virtual grid and it continues on the physical associated grid points. The creation of link array is time consuming but the gain in the search time is significant (at least a factor 20).

This algorithm is well known but the number of the virtual cells (N v ) is a parameter not clearly defined. A simple test is performed in order to determinate the optimum value of N v in function of the number of real coordinates N . We consider two uniform three-dimensional Cartesian grids (N x = N y = N z = N = 128) Figure 1: Main steps of the chimera method.

overlapping entire the cube 1x1x1 domain. Identical results were obtained on a coarser and a finer grid. The cpu-time is determined for range 1 ≤ N v ≤ 256. In addition to the compromise between cpu-time decrease for virtual grid creation and cpu-time increase for searching the overset, the size required by the link array must be taken into account. The results, represented in figure 2, show that the optimum value for the virtual cell number N v is the half of the number of physical coordinates N .The same test is performed for a grid with tangent hyperbolic distribution and, in this case, the optimum value correspond to N .

Figure 2: Representation of the cpu-time used for the indexation of real coordinates in virtual grid, cpu-time for searching overset and the size of the index between both grids versus the number of virtual box N v in one direction.

Another issue concerns the points that are inside a solid body. In this case, the cells are not detected as overlapped cells. To remove these cells, three methods are implemented. For simple cases (cylinder, sphere, channel) an analytical function is used. This solution is the quickest but the least flexible. For general cases, a second method is implemented using the closest wall. The dot product of the vector from the nearest wall cell center to cell center and the vector of the associated wall normal vector is calculated. If the dot product is negative the cell is in the solid region. This method requires the knowledge of the closest wall point and it is CPU time consumption. For the configuration of a sphere in a tube, having 6 millions of cells and 65 000 wall cells, the second method is 30 times more expensive than the first one. A third method are implemented to combine the rapidity and the flexibility of the both methods. The solution take a simple way : if the issue is the absence of solid volume, the addition of solid discretisation solves the problem. The figure 3 illustrates the concept by representing the added solid mesh associated to the slap of a high lift profil (MDA 30P-30N). The solid mesh fits the fluid mesh to share the same wall distribution but it requires any precaution as skew or size ratio. The detection follows the same process than for a fluid overlap. To simplify the storage and the reading, the solid grids are in the same format than the fluid mesh (structured defiinition). This third approach of solid detection gathers simplicity, rapidity and accuracy.

Figure 3: Scheme of the solid detection by added solid mesh for a high lift configuration, with the slat solid mesh and the fluid around the main profil (left) and the cells labeled as inside solid of the fluid mesh (right).

III.B. Determination of the type of overlapped cells

The second step of the chimera method is to define the type of the overlapped cells (calculated, interpolated, hole). In this step, among several overlapped cell, the best candidate to predict the flow state (calculated cell) and to give the flow solution to the other cells (interpolated cells) is determined. The criterion that justifies the choice is based on the best resolution of the physical problem. Four techniques are distinghuished, based on the criteria associated to the following quantities : user-defined grid ranking, user-defined cell quality, cell size and cell size in the wall-normal direction.

For a simple overset, the definition on an overlapped hierarchy of each block is enough for the cell type selection. All the cells of the highest ranking of the grid are defined as calculated cells. This solution is simple and fast but not suitable for complex overlapped layers where the definition of chimera grid ranking by the user is then needed. Based on the fact that the cell size is a criterion indicating the best discretization, Siikonen et al. [START_REF] Timo | Numerical techniques for complex aeronautical flows[END_REF] or Liao et al. [START_REF] Liao | A multigrid overset grid flow solver with implicit hole cutting method[END_REF] justify this choice by assigning the calculated cells type to the smallest ones. This criterion is locally adapted, automatic and requires no user input. A third criterion is added to the previous one and based on a user-defined cell quality. This way the user can force a specific kind of cell (calculated or interpolated) by associating to this cell a high quality value. This technique can force cell type by protecting or immunizing cells but the user input is binding and does not follow the grid movement. The cell size in the wall-normal direction was introduced by Landmann. [START_REF] Landmann | A highly automated parallel Chimera method for overset grids based on the implicit hole cutting technique[END_REF] The selection use the intersection segment between normal nearest wall and cell boundary. This technique is fully automatic and it is based on an important physical quality : the boundary layer. This test requires to search for the nearest wall and to calculate the intersection for each cell. In order to find a lighter, accurate and automatic criterion, we developed a similar criterion based on the distance of a cell to the nearest local wall. The nearest local wall is the wall in the same block where the cell is defined (figure 4 (left)). With this criterion, the overlapped cell with the smallest local near wall distance is calculated and the others are interpolated. This approach guarantees the resolution of Navier-Stokes equations in the nearest wall region and consequently each boundary layer is accurately computed. Moreover, this technique is based on the local wall distance and this quantity does not change with a moving grid (an update is however needed for deforming grids) whereas the criterion based on the wall normal intersection need an update after each grid movement. In the case of the sphere in a circular tube, the criterion based on the cell size complicates the generation of the mesh. The same problem may be encountered in the configuration of a sphere near a plane wall. The cell aspect ratio of the near wall mesh is not the same as that of the spherical wall mesh and the cell volume is very difficult to control. In our case of a sphere in a tube with a cylindrical grid overlapped by a spherical grid, the criterion based on the volume does not give a proper solution. The new criterion based on the local near wall distance give the best answer. The figure 5 represents the visualization of the calculated cells in the 3 planes crossing the sphere center for the both cases.

III.C. Buffer and hole cells

The next step is the addition of a buffer layer. The buffer cells are in between interpolated and calculated cells. Their role is to move away the interpolation communications and thereby to avoid implicit interpolation. An implicit interpolation is an interpolation that uses an interpolated cell as a donor cell. The more the interpolation needs a large number of donor cells, the more the buffer layer needs to be wide. The importance of wide overset is admitted and described in several papers. [START_REF] Sahu | Numerical simulation of three-dimensional transonic flows[END_REF][START_REF] Desquesnes | On the use of a high order overlapping grid method for coupling in CFD/CAA[END_REF] The set-up of the buffer layer is implemented with an user-input. The user give a fixed number of layer and the algorithm changes the type of the cell from interpolated to calculated over this layer.

The last type of possible cells is the hole cell. This kind of cell is a non-necessary cell and it concerns interpolated cells. This cell does not appear in any discretization schemes. In order to save computation time, the interpolation process is avoided for this cell. After the buffer layers, a user-input imposes a number of interpolated cells in the hole region. These new cells are considered as ghost-cells. The minimum number of ghost-cells is two because the numerical schemes used in NSMB require two neighbor cells (for most simulations a fourth order central scheme with artificial dissipation was used).

III.D. The interpolation

The interpolation is a key element of the chimera method. It allows the passage of flow information between overlapped grids and it creates the link. The influence of interpolation has already been investigated in several studies (Delfs, [START_REF] Delfs | An overlapped grid technique for high resolution CAA schemes for complex geometries[END_REF] Sherer & Scott 31 and Sengupta et al. [START_REF] Sengupta | Solving Navier-Stokes equation for flow past cylinders using single-block structured and overset grids[END_REF] ). The simplest treatment of an interface consists in the interpolation of all state variables. [START_REF] Benek | A flexible grid embedding technique with application to the euler equations[END_REF][START_REF] Hinatsu | Numerical computation of unsteady incompressible flow in complex geometry using a composite multigrid technique[END_REF][START_REF] Chesshire | Composite overlapping meshes for the solution of partial differential equations[END_REF][START_REF] Hs Tang | On nonconservative algorithms for grid interfaces[END_REF][START_REF] Freitas | Simulation of fluid-structure interaction using patched-overset grids[END_REF] However, this approach does not guarantee the global mass conservation as explained by Wang et al. [START_REF] Wang | A fully conservative interface algorithm for overlapped grids[END_REF] Some studies found a solution to ensure the global conservation for compressible flows [START_REF] Mohan | A relaxation approach to patched-grid calculations with the euler equations[END_REF] or for incompressible flows (Wright and Shyy [START_REF] Wright | A pressure-based composite grid method for the Navier-Stokes equations[END_REF] or Tang et al. [START_REF] Hs Tang | An overset-grid method for 3D unsteady incompressible flows[END_REF] ) with a correction of interpolated values. The accuracy of interpolation is important and different schemes of interpolation are Figure 5: Visualization of the calculated cells for the chimera cell selection before cell selection (1 st row), using cell volume criterion (2 nd row) and using the nearest local wall distance criterion (3 rd row) for the cross section of the tube (1 st line), the streamwise plane in the symmetric plane (2 nd line) and streamwise plane perpendicular to the symmetric plane (3 rd line) used : interpolation based on Lagrange polynomials, [START_REF] Desquesnes | On the use of a high order overlapping grid method for coupling in CFD/CAA[END_REF][START_REF] Sherer | High-order compact finite-difference methods on general overset grids[END_REF][START_REF] Sengupta | Solving Navier-Stokes equation for flow past cylinders using single-block structured and overset grids[END_REF][START_REF] Henshaw | A fourth-order accurate method for the incompressible Navier-Stokes equations on overlapping grids[END_REF][START_REF] Zang | A composite multigrid method for calculating unsteady incompressible flows in geometrically complex domains[END_REF] trilinear interpolation, 39, 42-44 quadratic interpolation [START_REF] Burton | Analysis of a fractional-step method on overset grids[END_REF] or tetravolumic interpolation. [START_REF] Timo | Numerical techniques for complex aeronautical flows[END_REF][START_REF] Ch Tai | Parallel computation of unsteady incompressible viscous flows around moving rigid bodies using an immersed object method with overlapping grids[END_REF] The interpolation scheme depends on the physical problem and the overset width. A highly-accurate interpolation needs many donor cells and it requires a wide overset. In our case of the sphere in a circular pipe, the gap between the sphere wall and tube wall does not provide a wide overset nor a large stencil of donor cells. For this reason, we chose some interpolation that use a low number of donor cells. The choice focus itself on the trilinear, tetravolumic and weighted inverse distance interpolations.

The inverse distance weighted interpolation is based on the distance between the points. The implementation is simple and its advantage is the flexibility of the stencil of points. The formula used is given by the following relation :

f M = N i=1 fi d α i N i=1 1 d α i (1)
where N is the number of the stencil points, f i the value at the point i, α is the weighting exponent (in our case equal to 2) and d i is the distance from the stencil point i to the interpolation point M .

The second interpolation is the tetravolumic interpolation based on the tetrahedral volumes composed between the four nearest cells. The interpolation weights are given by the opposite tetrahedral volume :

f M = i=A,B,C W i f i with W A = V BCM V ABC , W B = V ACM V ABC and W C = V ABM V ABC
with W the weight of the interpolation and V ACM the volume of the triangle (ABC). In three dimensions, the triangle volume is replaced by the tetrahedral volume.

The third interpolation is the trilinear one. We use a trilinear interpolation based on a linear interpolation on each direction. In two-dimensional example of the figure 6, the first linear interpolation yields the values f Q and f P then a second linear interpolation provides f M . In three-dimensions, three steps of linear interpolation are required. The developed algorithms were inserted in the École Polytechnique de Montréal 2D RANS code NSCODE. Verification of the interpolation schemes were made in Pigeon [START_REF] Pigeon | Développement d'une méthode d'accélération par grilles virtuelles récursives pour l'assemblage de maillages chimères[END_REF] and the methods robustness was shown for a multi-element airfoil slat/gap optimization problem. [START_REF] Levesque | An Overset Grid 2D/Infinite Swept Wing URANS Solver Using Recursive Cartesian Bucket Method[END_REF]

❡ ❡ ★ ★ ★ ★ ★ ★ ★ ★ ❍ ❍ ❍ ✏ ✏ ✏ ✏ ✏ ✡ ✡ ✡ ✡ ✡ 'receptor' cell center 'donor' cell center B A C M A C B D P M Q e 2

III.E. Application of the chimera method

The method described in the previous paragraph was validated using three different test cases : the simulation of flow around a sphere in an infinty domain, a flow around a sphere translating uniformly along a flat wall, and finally a sphere falling freely in a tube filled with a Newtonian fluid at rest. With each test case the complexity rises, so one part of a single flow, then an interaction between two walls defined by two independent chimera meshes, and finally the case of an interaction fluid / solid of a free falling body represented by on overset grid. These test cases were also selected for the availability of numerical or/and experimental results for comparison.

III.E.1. Flow past a sphere

The flow around a sphere is the topic of large number of experimental studies and numerical studies. [START_REF] Taneda | Experimental investigation of the wake behind a sphere at low Reynolds numbers[END_REF][START_REF] Rh Magarvey | Transition ranges for three-dimensional wakes[END_REF][START_REF] Provansal | Transition to turbulence in the wakes of axisymmetrical objects[END_REF][START_REF] Bouchet | Hydrodynamics forces acting on a rigid fixed sphere in early transitional regimes[END_REF][START_REF] Johnson | Flow past a sphere up to a Reynolds number of 300[END_REF][START_REF] Ghidersa | Breaking of axisymmetry and onset of unsteadiness in the wake of a sphere[END_REF][START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] The flow state is well known and clearly defined by a single parameter, the Reynolds number (Re = U d/ν, with U the freestream velocity, d the diameter of the sphere and ν the kinetic viscosity). For low Reynolds number (Re < 212) the flow is steady and axisymmetric, then for Re > 212 the axisymmetry is broken into a steady non-axisymmetric flow and the Hopf bifurcation appears at Re = 273 when the flow becomes unsteady with horse-shoe vortex shedding.

This case test is the first step of the validation and it demonstrates the ability and the accuracy of the chimera method for a simple three dimensional computation largely referenced. The simulations are performed by two set of grids : an chimera grid composed of a cartesian background grid overlapped by a spherical classical body-fitted mesh.

The representation of the isolines of streamwise velocity in the streamwise plane (figure 7) shows the wake behind the sphere with a recirculation. The shape of the vortex is a torus like described in the literature (Johnson et al. [START_REF] Johnson | Flow past a sphere up to a Reynolds number of 300[END_REF] for example). The dashed line on the figure 7 represents the external border of the spherical grid. The flow information crosses correctly this overset border.

The comparison of the results with results from literature is based on three physical quantities : the drag coefficient (C D ), the recirculation length (L) and the separation angle (θ). For this range of Reynolds numbers, the flow is axisymmetric and the lift force is zero. Comparing the drag coefficients with the results of Bouchet et al. [START_REF] Bouchet | Hydrodynamics forces acting on a rigid fixed sphere in early transitional regimes[END_REF] shows a good agreement, differences do not exceed 0.52% and this demonstrates the good resolution of the overlapped grids. The shape descriptions (L and θ) are in good agreement with the literature results and are shown in Fig. 8. The results are compared with datas obtained by experimental studies (Taneda, [START_REF] Taneda | Experimental investigation of the wake behind a sphere at low Reynolds numbers[END_REF] Nakamura 56 and Roos & Willmarth 57 ) and by numerical studies (Tomboulides et al., [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] Mittal, [START_REF]Planar symmetry in the unsteady wake of a sphere[END_REF] Bagchi & Balachandar 59 and Bouchet et al. [START_REF] Bouchet | Hydrodynamics forces acting on a rigid fixed sphere in early transitional regimes[END_REF] ).

III.E.2. Flow around a sphere translating along a plane wall

To complete the validation of the three-dimensional chimera method, we consider a configuration with an interaction between two walls yielding a border of overlapping.We simulate the flow past a sphere in a uniform translation parallel to a plane wall. This configuration has been previously simulated by Zeng et al. [START_REF] Zeng | Wall-induced forces on a rigid sphere at finite Reynolds number[END_REF] The parameters for this configuration are the distance between the sphere center and the wall (L/d) and the Reynolds number Re = U.d/ν, with U the velocity of the sphere translation and ν the kinematic viscosity of the fluid. The geometry is meshed with a Cartesian grid, which is refined along the plane wall, and a spherical grid which is refined near the sphere wall. The criterion for chimera cell selection is based on the definition of a boundary layer for the background grid where the cells are always calculated. Elsewhere the cells of spherical mesh are calculated. The motion of the sphere is implemented with the wall condition of the plate and using an inlet boundary condition. The mesh for parallel simulations is decomposed into 32 sub-blocks. The range of the Reynolds number is 10 ≤ Re ≤ 250 and the studied distance is L/d = 1 and L/d = 0.75.

The presence of the plane wall breaks the axisymmetric geometry and for all Reynolds numbers a lift force exists. Two mechanisms are the source of this force. The first is linked to the strong interaction between the wall and the wake of the sphere. Figure 9 shows the isovalues of the streamwise velocity in the wall normal plane for the case Re = 200 and L/d = 1.00. The distribution of the flow shows that the up/top asymmetric has been lossed, resulting in a lift force directed away from the wall. The second mechanism is linked to the high velocity in the gap associated to a low pressure. This behavior gives a force with opposite direction (attractive force) to the first mechanism. The resultant force tends to push the sphere away from the wall. The streamlines plotted along the streamwise wall normal plane (figure 10) represents the vortex structure for Re = 200 and L/d = 0.75. The same case performed by Zeng et al. [START_REF] Zeng | Wall-induced forces on a rigid sphere at finite Reynolds number[END_REF] is plotted below our results for comparison (Fig. 10). The overset border is not visible in the results, proving that the chimera implementation is correct. The results of the aerodynamic coefficients are in close agreement with those reported by Zeng et al, [START_REF] Zeng | Wall-induced forces on a rigid sphere at finite Reynolds number[END_REF] differences are below 1%, see Fig. 10. . The axisymmetry breaking arising in the wake generates a steady lift and the sphere is also exposed to the repulsive wakewall interaction. No stable state corresponding to static equilibrium of the lift and repulsion from the wall were found. Instead, the sphere tends to oscillate in a plane before switching to a helical trajectory. The latter presents the counterpart of the steady oblique trajectory of an unconfined sphere. The unconfined sphere undergoes a secondary instability at about G = 190. This instability leads to an unsteady periodic regime and to a chaotic trajectory (figure 4(right)).

The ARV was an ESA funded study for a both a cargo and crew transportation system to service the International Space Station (ISS). The project was abandoned in 2012. The ARV geometry is a blunt axisymmetrical capsule similar to the ARD (Advanced Re-entry Demonstrator) flown by ESA in 1998. The ARV-APEX cover is extracted from the module by two drogue parachutes that are modelled through forces acting on the APEX cover. The unsteady 6DOF CFD simulations were made using the chimera method on three different grids:

• a back ground mesh that is fixed in time;

• a mesh around the ARV capsule that moves in time following a prescribed initial velocity;

• a mesh around the ARV cover that moves in time using the the prescribed forces from the two drogue parachutes together with the aerodynamic forces acting on the cover.

Figure 11 shows a picture of these three chimera grids. The main question to be answered in these simulation was whether the separation process could lead to dynamic instabilities. Figure 12 shows the velocity field around the ARV and the APEX cover at four different instances in time, and one can observe the movement of the APEX cover due to the wake of the capsule.

IV. Turbulence modelling

IV.A. Current modelling in NSMB

The NSMB solver contains a large variety of turbulence modelling closures including LES, RANS, URANS and hybrid RANS-LES approaches. In the context of RANS-URANS, NSMB includes standard one-equation modelling as for example the Spalart-Allmaras 61 and its variants for improved near-wall treatment (Edwards and Chandra 62 ) and compressibility corrections of Secundov (Shur et al 63 , Spalart [START_REF] Spalart | Trends in turbulence treatments[END_REF] ), two-equation modelling including the k -ε and its various low-Reynolds number versions (example: the Chien modelling [START_REF] Chien | Predictions of channel and boundary-layer flows with a low-Reynolds number turbulence model[END_REF] ), the k -ω modelling of Wilcox [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] and the k -ω -SST 67 , as well as non-linear models among which the Abe, Jang and Leschziner 68 and the DRSM (Differential Reynolds Stress Modelling) of Launder, Reece anf Rodi [START_REF] Launder | Progress in the development of a Reynolds-stress turbulence closure[END_REF] and its low-Reynolds number version by Shima et al [START_REF] Shima | Low-Reynoldsnumberr Second-Moment Closure Without Wall-Reflection Redistribution Terms[END_REF] , as well as the SSG (Speziale, Gatski, Sarkar) modelling [START_REF] Speziale | Modeling the Pressure-Strain Correlation of Turbulence: an Invariant Dynamical Systems Approach[END_REF] . A specific attention has been paid in order to better capture non-equilibrium turbulence effects in high Reynolds number flows around bodies, where LES approaches are quite costly. NSMB involves the SAS-Scale Adaptive Simulation [START_REF] Menter | Scale-Adaptive Simulation Method for Unsteady Flow Predictions[END_REF] and the OES -Organized Eddy Simulation modelling [START_REF] Braza | Turbulence properties in the cylinder wake at high Reynolds number[END_REF][START_REF] Bourguet | Anisotropic organised eddy simulation for the prediction of nonequilibrium turbulent flows around bodies[END_REF][START_REF] Haase | DESider A European Effort on Hybrid RANS-LES Modelling: Results of the European-Union Funded Project[END_REF] including a tensorial eddy-viscosity concept derived from DRSM and able to account for near-wall turbulence stress anisotropy. OES modelling is particularly efficient in predicting the dynamics of unsteady coherent structures developing from the near-wall region and convected in the wake and ensures an improved prediction of the aerodynamic/hydrodynamic forces and of their fluctuations. Moreover, NSMB includes various hybrid RANS-LES methods, in order to take benefit from advanced URANS in the near-region and of LES in the outer region characterized by convected vortex dynamics farther from the wall. In this context, the DES (Detached Eddy Simulation 76 ), DDES (Delayed Detached Eddy Simulation 77 ) and IDDES (Improved Delayed Detached Eddy Simulation 78 ) approaches have been implemented. An example of DDES around a supercritical wing using NSMB can be found in the AIAA paper by Grossi et al [START_REF] Grossi | Prediction of transonic buffet by delayed detached-eddy simulation[END_REF] , in the context of the European research program ATAAC, Advanced Turbulence simulations for aerodynamic application challenges, http://cfd.mace.manchester.ac.uk/twiki/bin/view/ATAAC/WebHome. A considerable improvement of the detached eddy simulation has been accomplished within NSMB by implementing the OES modelling in the URANS part of the DES and DDES in case of fixed and moving bodies including fluid-structure interaction [START_REF] Shinde | Numerical simulation of the fluidstructure interaction in a tube array under cross flow at moderate and high Reynolds number[END_REF][START_REF] Moussaed | VMS and OES-based hybrid simulations of bluff body flows[END_REF] . The DDES-OES modelling has been coupled with the ALE -Arbitrary Lagrangian-Eulerian Method 81 accounting for mesh movement and deformation during the fluid-structure interaction process for moving/deforming configurations. The DDES-OES modelling has been also successfully applied in the simulation of the transonic buffet around supercritical wings in collaboration with Dassault Aviation in the context of the TFAST European research program, Transition location effect on shock wave boundary layer interaction, http://tfast.eu/ , (Szubert et al [START_REF] Szubert | Numerical study of the turbulent transonic interaction and transition location effect involving optimisation around a supercritical aerofoil[END_REF], where in addition, the transition location effect upstream of the shock-boundary layer interaction has been studied and an optimisation of its position has been achieved by means of kriging method. These recent developments have been able to capture the fluid-elastic instability leading to flutter in a cylinders bundle at high Reynolds number, (Shinde et al, [START_REF] Shinde | Numerical simulation of the fluidstructure interaction in a tube array under cross flow at moderate and high Reynolds number[END_REF] figures ?? and 13 left) as well as to predict the forces fluctuations around a wing of type Airbus A3xx in the context of the Smart Morphing Center activities, www.smartwing.org, in collaboration with Airbus -Emerging technologies and Concepts, Zilli, 2015. [START_REF] Zilli | Numerical Simulation and modeling of morphing deformation on a supercritical airfoil wing at high Reynolds number[END_REF] Among recent developments in turbulence modelling available in NSMB, a considerable effort has been paid in the development of stochastic forcing in the context of upscale turbulence modelling in order to enhance an eddy-blocking effect between the separated shear layers , by means of re-injection of turbulence from higher-order POD (Proper Orthogonal Decomposition) modes (Szubert et al [START_REF] Szubert | Shock-vortex shear-layer interaction in the transonic flow around a supercritical airfoil at high Reynolds number in buffet conditions[END_REF] ) -the IOES (Improved OES) approach. This method, inspired from inverse turbulence cascade and upscale turbulence modelling (Hunt et al [START_REF] Hunt | Vortical interactions with interfacial shear layers[END_REF] ) produces a thinning of the shear layers and a reduction of the wakes width with considerable improvement of the wall pressure prediction and of the aerodynamic forces (Szubert et al [START_REF] Szubert | Shock-vortex shear-layer interaction in the transonic flow around a supercritical airfoil at high Reynolds number in buffet conditions[END_REF] ). The OES method has been successfully applied in the case of the Dassault Aviation V2C wing (TFAST European program), in comparison with the DDES method (Szubert et al [START_REF] Szubert | Numerical study of the turbulent transonic interaction and transition location effect involving optimisation around a supercritical aerofoil[END_REF] ), where a considerable improvement of the forces and suction effect has been achieved. A view of the POD modes used for the stochastic forcing in the V2C transonic buffet case are presented in the figure 14.

Moreover, a specific attention has been paid in implementing transition modelling in NSMB, in the context of the PhD thesis of Damien Szubert (Szubert, 2015 86 ). The γ -R θ transition modelling of Langtry and Menter [START_REF] Langtry | Correlation-Based Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes[END_REF] has been implemented and led to simulation of transonic buffet around the V2C -3D configuration within the TFAST project (figure 13 right).

IV.B. Simulations of turbulent flows

The following results concern the recent developments in turbulence modelling within NSMB. The DDES-OES modelling has been employed in a tandem of cylinders configuration for the prediction of the aerodynamic forces, predominant frequencies associated with landing gear noise, as well as the MIV -Movement Induced Vibration developing by the downstream cylinder due to the shearing mechanism created by the upstream cylinder at Reynolds numbers in the range 100 000 -170 000. An analysis of the flow dynamics for the static configuration can be found in an experimental study of Elhimer et al [START_REF] Elhimer | Coherent and turbulent processes in the bistable regime around a tandem of cylinders including reattached flow dynamics by means of high-speed PIV[END_REF] by means of Time-Resolved PIV carried out in the S4 wind tunnel of IMFT.

The mesh used is of 12 M and was provided by NTS -St Petersbourg in the context of the ATAAC EUprogram. Figure 15 shows instantaneous vorticity fields colored by pressure, at organized large wavelength. The spectra in figure 16 show the development of the von Kármán mode and of the Kelvin-Helmholtz vortices mode in comparison with experiments by Jenkins et al [START_REF] Jenkins | Characterization of Unsteady Flow Structures Around Tandem Cylinders for Component Interaction Studies in Airframe Noise[END_REF] . Figure 16 (left) shows the turbulence structure developed after free vertical vibration (in the direction of lift force), spontaneously developed by the downstream cylinder allowed for one-degree-of-freedom vibration, when the reduced velocity (u * = U inf /(f 0 D)) becomes higher than 3, U inf being the upstream velocity and f0 the structural eigenfrequency. The simulations have been carried out for a Scruton number (mass-damping) of 1. Figure 16 (right) shows the progressive amplification towards negative damping (i.e. flutter) of the downstream cylinder's motion, as the reduced velocity increases. This evolution also shows the progressive increase of the phase-lag as the fluid-structure interaction regime leaves out the VIV (Vortex Induced Vibration and enters to MIV (Movement Induced Vibration) behaviour for reduced velocities higher than 3.

V. Icing modelling

In-flight aircraft icing is a main research area within the aeronautic community. Ice accretion mostly results in severe aerodynamic degradation and is the major cause of flight incidents. Because of the rigorous certification process and since all configurations cannot be experimentally tested, aircraft in-flight icing must be simulated via computational fluid dynamics.

By impacting on aircraft structures, the supercooled water droplets carried by air can immediately freeze upon impact or flow as a water film depending on the surface temperature. Two main types of ice known as rime ice and glaze ice are formed depending on the environmental conditions.

The dry rime ice usually forms at low flow speeds and low flow temperatures when the entire droplet freezes at its impact location. Higher flow temperatures and flow speeds induce the formation of wet glaze ice. In that situation a fraction of the droplet freezes at impact point and the other portion flows as a thin water film runback and freezes downstream. The modification of the geometry due to icing results in large aerodynamic degradation.

Ice accretion codes used by industries are commonly based on a 2D inviscid panel method for the flow simulation, a boundary layer method based on empirical relations for the heat transfer coefficient evaluation, a Lagrangian formulation to simulate the droplet trajectories and the so-called Messinger model to predict the ice thickness. A better accuracy can be obtained by repeating this one-way coupling process until the desired exposure time (multi-step approach). [START_REF] Hasanzadeh | Quasi-steady con-vergence of Multistep NavierStokes Icing Simulations[END_REF] In NSMB, the supercooled water droplets field and impingement locations are obtained by means of an Eulerian approach similarly to FENSAP-ICE [START_REF] Beaugendre | Development of a sec-ond generation in-flight icing simulation code[END_REF] and ONICE3D. [START_REF] Montreuil | ECLIPPS: 1. three-dimensional CFD prediction of the ice accretion[END_REF] The ice thickness is calculated by resolving the Shallow Water Icing Model (SWIM) for three-dimensional and multi-block ice accretion simulations. [START_REF] Beaugendre | Development of a sec-ond generation in-flight icing simulation code[END_REF] A Mesh deformation is used to track the ice/air interface through time for multi-steps calculations. These approaches allows a seamless integration to Navier-Stokes solver as the entire process is computed on the same mesh used for the fluid computation and parallelized with Message Passing Interface for complex multi-block configurations.

V.A. Modelling of Droplet Field

An Eulerian approach is developed for modelling droplet trajectories and impingement on a solid structure. The Eulerian model solves a droplet velocity field and a volume fraction distribution in the whole computational domain. Droplets are considered spherical and not to be subject to deformation or breaking. The droplet phase is sufficiently diluted to neglect interactions between droplets. Temperature of the droplet phase is set to be constant and equal to free-stream temperature, heat transfer with the surrounding air is neglected. The only forces acting on droplets are drag, gravity and buoyancy.

The governing equations (which correspond to an Eulerian formulation of a Lagrangian transport of particles) for the conservation of mass and momentum of the droplets are written as follows 94 :

∂α ∂t + ∇.(αu) = 0 ∂αu ∂t + ∇. (αu ⊗ u) = α C D Re d 24K (u a -u) + α 1 -ρa ρ 1 F r 2 g (2)
α is defined as the non-dimensionalized volume fraction of water and u the non-dimensionalized velocity field of droplets. u a is the non-dimensionalized velocity of air, ρ the density of water, ρ a the density of air, g the gravity vector. F r = U ∞ (Lg) is the Froude number, U ∞ the freestream velocity of air, L the characteristic length (typically the airfoil chord length), K = ρdU ∞ /18Lµ is an inertia parameter and µ the dynamic viscosity of air. d is the median diameter of the droplets. The first term on the right-hand-side of the momentum equation accounts for the drag acting on the droplet or particle based on low-Reynolds number behaviour for spheres. The droplets Reynolds number ( Re d ) is defined based on the slip velocity between the air and droplet and the droplet diameter.

Re d = ρd|u a -u|U ∞ µ (3)
The drag coefficient of the droplets C d is given by :

     C d = 0.44 Re d , if Re d > 1000 C d = 24 Re d , if Re d < 0.1 C d = 24×(1+0.15×Re 0.657 d ) Re d (4)
A third order upwind scheme is used and the implicit droplet equations are solved with a BiCGSTAB algorithm. The freestream values of droplet velocity and volume fraction are imposed as boundary conditions at the far field. A switching boundary condition is applied at the wall : Neumann when the incoming droplet flux is positive and Dirichlet otherwise. Finally the collection efficiency coefficient that characterizes the configuration's ability to capture incoming water droplets on the solid surfaces is defined as :

β = αu.n (5)

V.B. Modelling of Ice Thickness

To calculate the properties of the ice and the liquid water film on the solid wall the mass and energy conservations are solved. The input to the equations are the wall shear stress τ wall and the heat transfer coefficient h c obtained by the flow solver and the collection efficiency β and the droplet velocity field u computed by the droplet module. The Shallow Water Icing Model icing model (SWIM) implemented is based on the Messinger model ( 95). The Mass and energy conservation equations are solved for each wall cell. Radiation is neglected in the energy conservation equation as it mostly accounts for a low contribution.

The system of equation writes :

         ∂ρ w h ∂t + ∇.(ρ w ūf h) = ṁimp -ṁevap -ṁice ∂ρ w C w hT ∂t + ∇.(ρ w C w ūf hT ) = qimp -qevp -qice -qcnv (6)
h is the water film height, ρ w the water density, C w the specific heat of water, T the surface temperature in Celsius, ṁevp the mass rate of evaporation, ṁimp the impinging mass rate of water droplets and ṁice the resultant mass rate of ice accretion. qimp is the heat rate of impacting droplets, qevp the evaporating heat rate and qice the heat rate of freezing water. The main contribution comes from the heat convection term qcnv .

The impinging mass rate contribution writes ṁimp = LW C. V ∞ . β and the evaporating mass rates ṁevp are determined by a parametric model developed by Kreith et al. [START_REF] Kreith | Principles of heat transfer[END_REF] The heat rate contributions are :

             qimp = ṁimp C w T ∞ + |u| 2 2 qevp = ṁevp Le+Ls 2 qice = ṁice (C i T -L f ) qcnv = h c (T -T ∞ ). (7)
C i is the specific heat of ice, L e is the latent heat of evaporation, L s is the latent heat of sublimation and L f is the latent heat of fusion.

ūf is defined as the mean velocity of the water film and writes [START_REF] Beaugendre | Development of a sec-ond generation in-flight icing simulation code[END_REF] :

ūf (x) = h 2µ w τ wall (x), (8)
with τ wall the wall shear stress, x the wall surfaces coordinates and µ w is the viscosity of water.

To insure a physical solution and close the system, the following compatibility relations must be verified 92 :

             h ≥ 0, ṁice ≥ 0, hT ≥ 0, ṁice T ≤ 0. (9)
The heat transfer coefficient is a key parameter in icing simulations. In NSMB, this coefficient is evaluated from the heat fluxes q wall by imposing a constant wall temperature T wall . The analytical definition of the heat transfer coefficient writes :

h c = q wall T ∞ -T wall . (10)
The heat transfer coefficient is highly dependent to the chosen roughness parameter. As a first approximation, the constant equivalent sand grain roughness parameter used in this paper is determined by the correlation from Shin and Bond. [START_REF] Shin | Experimental and computational ice shapes and resulting drag increase for a NACA 0012 airfoil[END_REF] The convective terms of the Eq. 6 are discretized via finite volume method with a first order Roe scheme. The system of two equations is solved explicitly in time via a Runge-Kutta scheme and is parallelized with MPI.

V.C. Results

Rime icing is simulated on a NACA23012 airfoil and compared with experimental result from Broeren et al. [START_REF] Broeren | Effect of High-Fidelity Ice-Accretion Simulations on Full-Scale Airfoil performance[END_REF] The O-grid mesh is generated by NS-GRID software presently developed by Ecole Polytechnique de Montreal. NS-GRID is an elliptic smoothing mesh generation code that control orthogonality, skewness and spacings via control functions. [START_REF] Hasanzadeh | Quasi-steady con-vergence of Multistep NavierStokes Icing Simulations[END_REF] The NACA23012 mesh has a chord of c = 1.828m and contains 512 × 256 cells divided into two blocks with a farfield extended to 50 chords.

The parameters of the simulations are an angle of attack of 2.0 • , the Mach number is M a = 0.24, the Reynolds number is Re = 12.60 × 10 6 , the freestream temperature is T = 247.85K, the median droplet diameter is d = 30µm, the liquid water content is LW C ∞ = 0.55g/m 3 and the total icing time is 600s. Visualization of Liquid Water Content distribution is shown in Fig. 19. Ice shape obtained after a single icing step are compared with experimental results in Fig. 20. The simulated ice shape is in excellent agreement with experimental data. One should not that taking into account a droplet size distribution to simulate the water droplet spray should provide even more precise results on the upper impingement limit where a slight underestimation of the ice thickness is observed.

Glaze icing is simulated on a two-dimensional NACA0012 airfoil and compared with NASA LEWICE2D 99 and the experimental data available. 100 The NACA0012 O-grid mesh was generated by NS-GRID and has a chord of c = 0.5334m. The mesh contains 512 × 256 cells divided in two blocks. The constant sand grain roughness used in the computation is set to k s = 0.55mm. Calculation is performed on two cores. Impingement efficiency distribution curve along the airfoil surface is shown on Fig. 21. Compared to LEWICE, present results showcase a slight higher value of the lower impingement limits. This observation is in accordance with results of several studies of the Eulerian droplet model. 94,101 The calculated glaze-type ice shape is generated with three icing steps and compared with available data in Fig. 22. This test was repeated on a NACA00012 extruded mesh. Identical results were obtained along wing span partially validating the three-dimensional implementation. The last test case is a three-dimensional simulation performed on an ONERA-M6 swept wing. The 4 blocks mesh contains 884736 cells. The mean aerodynamic chord is c = 0.53m, the semi-span is b = 1m and the computational domain extends to 20 × c. The flow is simulated at an angle of attack of 6 • , a Mach number of M = 0.15, a Reynolds number of Re = 2.17 × 10 6 , a static temperature of T = 265K and with a constant equivalent sand grain roughness number of k s = 0.55mm. The parameters of the water spray are a mono-dispersed droplet diameter of d = 20µm and a liquid water content of LW C = 1.0g/m 3 . Icing is simulated at a total time of t = 500s. Impingement curves are plotted on three sections A, B & C respectively at 20%, 60% and 90% of span length and compared with Lagrangian results from Sang et al. 102 in Fig. 23. The present results show a good agreement. The three-dimensional ice shape at exposure time t = 500s is represented in Fig. 24. According to our knowledge there is unfortunately no experimental or numerical data regarding ice shapes to compare with on this configuration.

VI. 2nd High Lift Prediction Workshop

NSMB was used to perform studies in the 2nd AIAA High Lift Prediction Workshop (HiLiftPW2) to better understand the influence of the numerical schemes on numerical predictions. 103 The model is the DLR-F11 model (Figure 25), a semi-span three element high lift wing representative for a commercial wide-body twin-jet. 104 While all participants ran the mandatory test cases in RANS mode, NSMB was one of few that ran URANS simulations. This allowed characterizing the flow in steady or unsteady states. In addition, emphasis was placed on the study of effects of initial conditions on final solution.

The grid is the original 1-to-1 structured grid provided by the HiLiftPW-2 committee. Three grid sizes are proposed: coarse, medium and fine with respectively 9.5, 32 and 101 million cells. Each grid level is obtained by a uniform growth of around 1.5 time in each coordinate direction. The non-dimensional normal wall spacing of the first cell is 1.88 10 6 , 1.29 10 6 and 0.83 10 6 respectively for the coarse, medium and fine grid, and the farfield is located at around 100 reference chords.

NSMB settings were : second order central scheme for diusive terms, fourth order central scheme with artficial dissipation of Jameson for convective terms. For unsteady calculations, a dual time stepping scheme with backward second order interpolation is used. The system is solved with an implicit LU-SGS (Lower-Upper Symmetric Gauss Seidel) scheme with 100 inner iterations with 2 level multigrid. The nondimensionalized timestep was 0.06. Three variant of the fully turbulent Spalart-Allmaras turbulent model are tested : the original Spalart-Allmaras and the Edwards and SALSA modifications.

Here, only results concerning the effect of initial conditions are discussed. The dataset is as listed in table 1, all referring to the medium mesh and using the steady solver with the SA-Edwards turbulence model. The trend of the lift coefficient is plotted in figure 26 for the three initial conditions. The difference between the initial conditions occurs only after α = 21 o . After this angle-of-attack, the lift is overpredicted with the results initialized from freestream, whereas the opposite occurs with the initial condition set at 16 degrees. The global trend of the lift differs from that observed in the experiment. The stall curve presents a drought around 20 degrees before the maximum lift for all initial conditions. An examination of the wing surface pressure distributions shows that this drought is related to the fuselage/wing junction. A linear eddy viscosity model using the Boussinesq relation can be replaced by a nonlinear stress term like the Quadratic Constitutive Relation (QCR), a scheme described by Spalart

Angle of attacks

VII. Recent developments for hypersonic flow simulations

NSMB code was from the start of its development designed for use for hypersonic applications (in particular re-entry space vehicles). NSMB includes different levels of chemistry modeling for Air, CO2 and N2. Available models for air consider 5, 7, 11 or 13 chemical species, and a large number of non-equilibrium chemistry models are available, among them the models of Park 1989, 106 Park 1993, 107 Park 2001, 108 Gar-diner 109 and Dunn and Kang 110 are available. For all these models the equilibrium constants for the chemical reactions are computed using the polynomials provided in the references. Chemical equilibrium can be computed using polynomials (TGAS 111 /VGAS 112 ), or by calculating the chemical composition for equilibrium. Thermochemical non-equilibrium has been implemented using the Landau-Teller equation for the translationalvibrational energy exchange with the relaxation time calculated using the semi empirical Millikan and White formula. All these models are available for inviscid and viscous (both laminar and turbulent) flows.

Transport coefficients are computed using the Blottner model 113 or the Gupta-Yos model. 114 Different methods of catalytic wall boundary conditions have been implemented, including super catalytic walls and partially catalytic walls.

NSMB includes an option to adapt the far field boundary to the bow shock. This option removes the so-called carbuncle phenomenon introduced by the numerical scheme due to the non-alignment of the grid with the bow shock. The bow shock adaption procedure requires a C-type of grid, and the shock should be located in the blocks with a far-field boundary condition. During the solution process the position of the shock on a grid line is determined using the Mach number as variable (other variables were used but were less accurate or less robust). Once the shock position is known the far field boundary is moved along the original grid line to a position near the shock location. Then grid lines are clustered around the shock position, followed by a redistribution of the grid points on the remainder of the original grid line. In general the grid needs to be adapted two to three times to obtain a stable shock position. For capsule type geometries the bow shock grid adaption improves considerably the solution on the heat shield, and in particular the heat flux. Figure 27 shows the Mach number contours in the symmetry plane for the ARV capsule, and one clearly observes the improvement due to the grid adaptation. Figure 28 shows the heat flux contours on the heat shield for a chemical non-equilibrium calculation at M ach = 25.1 for the original and adapted grid. One clearly observes the improvements in the heat flux due to the adaption of the grid to the bow shock.

Roughness and blowing modeling was implemented in NSMB during the European Union funded project RASTAS-SPEAR (Radiation-Shapes-Thermal Protection Investigations for High Speed Earth Re-entry),. [START_REF] Vos | Earth re-entry capsule CFD simulations taking into account surface roughness and mass injection at the wall[END_REF] For roughness modeling the so-called equivalent sand-grain approach is used, which consists of two steps. In the first step the rough surface is reduced to an equivalent reference surface, in the second step the turbulence model is modified to reproduce the effect of rough surfaces on the skin friction and heat transfer. Reviews of surface roughness modeling can be found in. [115][116][117] Both the Spalart-Allmaras turbulence model 61 and the k -ω family of turbulence models [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF][START_REF] Menter | Zonal Two-equation k-w Turbulence Models for Aerodynamics Flows[END_REF] were modified to account for rough walls. Calculations were made for the Earth Re-entry Vehicle shown in Fig. 28 using the Spalart-Allmaras turbulence model for a smooth wall for a rough wall using 2 roughness heights (250 and 500 µm). As can be seen in Fig. 29, for a laminar flow the computed heat flux is in excellent agreement with the heat flux obtained from empirical correlations. For a turbulent flow using the Spalart Allmaras turbulence model the computed heat flux follows the laminar result until y = 0.23 m, indicating that the flow at the nose is laminar and that transition to turbulence takes place further downstream. Once transition takes place the heat flux increases rapidly. This increase in heat flux is faster and higher when using roughness modeling. The turbulent correlations also assume a laminar flow at the nose (until y = 0.08 m), followed by a turbulent flow further downstream. In a recent project funded by the European Space Agency on Ablation modeling NSMB was coupled to the Mutation++ library developed at the Von Karman Institute. 121 This library is written in C++ and designed to provide efficient algorithms for the computation of various properties important in the fields of hypersonics and combustion Computational Fluid Dynamics (CFD) including:

• Thermodynamic Properties using -The statistical mechanics approach using partition functions for the different energy modes of a species. Rigid-Rotor and Harmonic Oscillator models have been implemented

-Multi-temperature lookup tables for high fidelity models based on quantum mechanics -High fidelity curve fits, and in particular the NASA 7 and 9 coefficient polynomial data bases

• Multicomponent Transport Properties -Use of collision integrals for each collision pair of species in the mixture -Wilke, Gupta-Yos, Armaly-Sutton and Blottner mixture rules -State of the art methods to solve the system, LDLT for direct solution and Conjugate-Gradient for iterative solution

• Finite Rate Chemistry in Thermal Non-equilibrium

• A highly robust Multiphase Equilibrium solver (not used by NSMB).

In the next phase of the project an engineering model will be developed to account for both carbonbased TPS ablation and catalysis. These models will be updated and refined using experiments carried out at the Von Karman Institute in Belgium. Once this process is complete these models will be included in Mutation++, and can then be used by the NSMB CFD code to perform CFD simulations of re-entry vehicles using ablative TPS materials.

VIII. Conclusion

Since the last AIAA paper describing the NSMB solver 1 , the NSMB developpers consortium has been reduced but the development of the this code is still on-going and many major improvements have been performed in the last few years : an fully automatic chimera method for moving bodies, an icing model, many turbulence models for unsteady flows and new models for hypersonic flows. All these developments have been successfully applied to various flow configuration allowing the NSMB solver to be one of the most powerfull research CFD solver in Europe.

Figure 4 :

 4 Figure 4: Left : Definition of local nearest wall distance and global nearest wall distance for the cell C 1 of the mesh 1 in the configuration of two overlapped polar meshes; Right : Flow structure of the free falling of a sphere in a tube.

1 Figure 6 :

 16 Figure 6: Tetra-volumic (left) and tri-linear interpolation (right) interpolation schemes.

Figure 7 :

 7 Figure 7: Isovalues of the streamwise velocity component (left) and the streamlines (right).

Figure 8 :

 8 Figure 8: Drag coefficient C D , recirculation length (L/d) and separated angle versus Reynolds number for the flow past a sphere and 50 ≤ Re ≤ 200.

Figure 9 :

 9 Figure 9: Isovalues of the streamwise velocity in the plane normal to the wall (-0.3 to 1.1 with a increment of 0.2)

Figure 10 :

 10 Figure 10: The drag and lift coefficient for L/d = 1.00 and versus Reynolds numbers and b,c) streamlines in the plane normal to the wall for the case Re = 200, L/d = 0.75 for the present study (left) and extract to Zeng et al. 60 (right)

Figure 11 :

 11 Figure 11: Chimera grid used in the ARV APEX cover separation simulations.

Figure 12 :

 12 Figure 12: Velocity field around the ARV and APEX cover at t = 0.2, 0.4, 0.6 and 0.8 secs.

Figure 13 :

 13 Figure 13: Left : Development of vibratory instability and flutter as the reduced velocity increases; Right : λ 2 criteria colored by µ t /µ of the flow around the Dassault Aviation V2C 3D wing, M = 0.70, Reynolds around 4.10 6 and α = 7 o .

Figure 18

 18 Figure 18 shows the prediction of oblique choc reflexion over a boundary layer, according to the TFAST test-case measured by Technical University of Delft. The mesh is of 30 M points and a first overview of the computations can be found in the CTR Biannual report of 2014 by Szubert et al 90 in collaboration with the group of Prof. P. Moin.

Figure 14 :

 14 Figure 14: POD reconstructed instantaneous turbulent velocity fields of the flow around the Dassault Aviation V2C aifoil, M = 0.70, Re = 3.245 10 6 and α = 6 o , DDES-kω SST.

Figure 15 :

 15 Figure 15: Simulation of the turbulent flow around a tandem of cylinders at Reynolds number 160, 000 and centers distance 3.75 D. Illustration of turbulent motion development and of streamwise vorticity filaments formation (blue structures) arranged spanwisely periodic, according to a large-scale wavelength, (Test-case of the ATAAC European program).

Figure 16 :

 16 Figure 16: Power Spectral Density of the wall pressure at four points located at θ = ±45 o on the cylinders surface and comparison of normalised spectral density with experimens by Jenkins et al 89 .

Figure 17 :

 17 Figure 17: Evolution of the downstream cylinders position as a function of time for increasing reduced velocity and Scruton number of 1. Progressive appearance of phase-lag characterizing the MIV behaviour and negative damping.

Figure 18 :

 18 Figure 18: Oblique shock interaction, upstream Mach number 1.7 and Reynolds number of 36 M. IDDES and DDES computations using the Spalart-Allmaras model. Comparison of the boundary layer shape factor with PIV experiments of the Technical University of Delft (R. Giepman, F. Schrijer and B. van Oudheusden).

Figure 19 :

 19 Figure 19: NACA23012 -Visualization of Liquid Water Content.5

Figure 20 :

 20 Figure 20: NACA23012 -Comparison of ice shape with experimental data from, 985

Figure 21 :

 21 Figure 21: NACA00012 -Comparison of collection efficiency with LEWICE on clean geometry.5

Figure 22 :

 22 Figure 22: NACA00012 -Comparison of ice shape with experimental data from. 100

Figure 23 :

 23 Figure 23: Collection efficiency on an ONERA M6 swept wing 5

Figure 24 :

 24 Figure 24: Icing on a ONERA M6 swept wing

Figure 25 :

 25 Figure 25: DLR-F11 geometry.

  and applied to the HiLiftPW-2 by Pulliam & Sclafani 105 . a) b) Figure 26: Lift coefficient (a) and its amplitude (b) versus the angle of attack for several initial conditions.

Figure 27 :

 27 Figure 27: Mach number contours symmetry plane ARV capsule at M ach = 1.47, left unadapted grid, right adapted grid.

Figure 28 :

 28 Figure 28: Mach numbers in the symmetry plane (top row) and heat flux on the heat shield (bottom row), unadapted grid (left) and twice adapted grid (right). Earth Re-entry Vehicle at M ach = 25.1, h = 48 km, chemical non-equilibrium simulation using a 13 species air model.

  Figure 29 shows the computed heat flux along the wall in the symmetry plane. No blowing was used for these calculations. The computed results are compared with empirical correlations for a laminar flow 118 and a turbulent flow.119 Roughness effect are taken into account in the empirical correlations using the formulation of Powars.120

Figure 29 :

 29 Figure 29: Influence of surface roughness on the heat flux, Earth Re-entry Vehicle at M ach = 25.1, h = 48 km, chemical non-equilibrium calculation.

Table 1 :

 1 Dataset for initial condition study.

	Initial conditions

Table 2 :

 2 Table 2 summarizes the increase in heat flux compared to a smooth wall at 4 different points. Influence roughness on computed heat fluxes.

		Mid Sphere	Sphere-Cone	Mid-cone	Tore
			Junction
		Y = 0.1 m	Y = 0.2 m	Y = 0.37 m
	Spalart Rough 1	5%	2%	25%	27%
	Spalart Rough 2	5%	2%	71%	54%
	Correlation Rough 1	21%	21%	20%	22%
	Correlation Rough 2	39%	42%	41%	43%

of 28 American Institute of Aeronautics and Astronautics

of 28 American Institute of Aeronautics and Astronautics


Meriam Attia

Ouidad Yousfi

Abdelwahed Omri

Does Board Structure Matter for Innovation?


A Paolone

E Placidi

E Stellino

M G Betti

E Majorana

C Mariani

Rago

A Nucara

O Palumbo

P Postorino

I Rago

F Trequattrini

M Granata

J Teillon

D Hofman

C Michel

A Lemaitre

N Shcheblanov

G Cagnoli

F Ricci

Effects of the annealing of amorphous Ta 2 O 5 coatings produced by ion beam sputtering concerning the effusion of argon and the chemical composition

Keywords:

come

Effects of the annealing of amorphous T a_2O_5 coatings produced by ion beam sputtering concerning the effusion of argon and the chemical composition

Introduction

Amorphous Ta 2 O 5 is a material with a large variety of applications in electronic circuits [START_REF] Li | Ultralow power switching of Ta 2 O 5 /AlO X bilayer synergistic resistive random access memory[END_REF], quantum technologies [START_REF] Splitthoff | Tantalum pentoxide nanophotonic circuits for integrated quantum technology[END_REF], gas [START_REF] Liu | New room temperature ammonia gas sensor synthesized by a tantalum pentoxide (Ta 2 O 5 ) dielectric and catalytic platinum (Pt) metals[END_REF] or molecular sensing [START_REF] Magesa | Electrochemical sensing fabricated with Ta 2 O 5 nanoparticle-electrochemically reduced graphene oxide nanocomposite for the detection of oxytetracycline[END_REF], enhancement of Raman scattering [START_REF] Yang | Charge transfer enhancement in the surface-enhanced Raman scattering of Ta 2 O 5 superstructures[END_REF], hydrogen generation activity [START_REF] Liu | Formation and characterization of gray Ta 2 O 5 and its enhanced photocatalytic hydrogen generation activity[END_REF] or hydrophobic activity [START_REF] Sharma | Fabrication and estimation of optical constants of hydrophobic Ta 2 O 5 thin films[END_REF]. Among the physical properties which allow this broad range of applications, the refractive index n~2.1 at 1064 nm is of the highest importance [START_REF] Amato | Highreflection coatings for gravitational-wave detectors: state of the art and future developments[END_REF]. Actually, multilayers reflecting coatings and photonic crystals are conceived and realized combining films of different refractive index in stacked sequences. This possibility has been largely explored for the construction of the coatings of the most demanding optical components ever realized: the mirrors of the gravitational-wave (GW) detectors [START_REF] Amato | Highreflection coatings for gravitational-wave detectors: state of the art and future developments[END_REF][START_REF] Granata | Internal friction and Young's modulus measurements on SiO 2 and Ta 2 O 5 films done with an ultra-high Q silicon-wafer suspension[END_REF][START_REF] Steinlechner | Thermal noise reduction and absorption optimization via multimaterial coatings[END_REF][START_REF] Granata | Mechanical loss in state-of-the-art amorphous optical coatings[END_REF][START_REF] Abernathy | An overview of research into low internal friction optical coatings by the gravitational wave detection community[END_REF][START_REF] Amato | Optical properties of high-quality oxide coating materials used ingravitational-wave advanced detectors[END_REF]. At the present state, the mirrors of the Advanced Virgo and Advanced LIGO detectors are large masses of silica (35 cm diameter, 40 kg mass) with a high reflecting coating composed of alternate layers of amorphous SiO 2 and Ta 2 O 5 [START_REF] Amato | Highreflection coatings for gravitational-wave detectors: state of the art and future developments[END_REF][START_REF] Acernese | status of advanced virgo[END_REF][START_REF] Acernese | advanced virgo status[END_REF].

In this particular application, besides the optical properties, the most important parameter is the mechanical energy loss, because it is the key limiting factor of the interferometer performances due to the associated thermal noise [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF]. Despite the small thickness of the coating (≤6 μm) compared to that of the bulk mirror (about 20 cm), the main source of mechanical energy loss is indeed the thin film, in particular the Ta 2 O 5 layers [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF]. Two principal strategies to reduce the coating loss have been proposed: doping with amorphous TiO 2 or post deposition thermal treatments [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF][START_REF] Harry | Titania-doped tantala/silica coatings for gravitational-wave detection[END_REF]. In the following we will investigate in more detail some effects of the second procedure. Previous reports evidenced that Ta 2 O 5 has a mechanical energy dissipation about 10 times higher than SiO 2 [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF]. The exact value for Ta 2 O 5 depends on the deposition technique; however, after a thermal treatment for 10 hours in air, the elastic dissipation reduces and reaches the same value, independently of the synthesis method [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF]. This procedure has been largely used for the production of the mirrors of the gravitational-wave detectors. Also the refractive index of Ta 2 O 5 slightly depends on the synthesis method; however, the small differences seem not to be erased by thermal treatments [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF]. Notably, longer thermal treatments do not have an impact either on the elastic loss or on the refractive index [START_REF] Amato | Highreflection coatings for gravitational-wave detectors: state of the art and future developments[END_REF].

While the macroscopic effects produced on the optical and mechanical properties by the thermal treatments have been largely investigated, much less is known about the induced microscopic changes. An investigation of Ta 2 O 5 by Extended X-ray Absorption Fine Structure (EXAFS) measurements evidenced that samples annealed in air at 300 and 600 • C displayed only marginal changes of the amorphous structure, when it is investigated on a length scale of 50 nm [START_REF] Bassiri | Order within disorder: the atomic structure of ion-beam sputtered amorphous tantala (a-Ta 2 O 5 )[END_REF]. These results are in agreement with Raman spectroscopy measurements that could not reveal important changes of the Ta 2 O 5 structure on a macroscopic length scale induced by thermal treatments [START_REF] Amato | Highreflection coatings for gravitational-wave detectors: state of the art and future developments[END_REF]. However, Fluctuation Electron Microscopy (FEM) measurements on Ion Beam Sputtered amorphous tantala coatings revealed a medium range order on a typical length scale of 2 nm; this spatial scale increases up to 5 nm in samples thermally treated at 600 • C [START_REF] Hart | Medium range structural order in amorphous tantala spatially resolved with changes to atomic structure by thermal annealing[END_REF], even though a macroscopic crystallization is not revealed at this temperature.

In the present work, we further investigate the effects of thermal treatments on amorphous Ta 2 O 5 produced by Ion Beam Sputtering (IBS), as this production technique is the method of choice to obtain the highly uniform coatings needed for gravitational-waves detectors [START_REF] Amato | Optical properties of high-quality oxide coating materials used ingravitational-wave advanced detectors[END_REF]. Some aspects are deeply investigated: the first is the possible retain of argon atoms in the coating during the synthesis procedure. In fact, in the IBS deposition, an Ar ion beam is used to produce the sputtered matter to be deposited on the substrate and some atoms could be retained either at the coating/substrate interface or in the whole thickness of the film. Moreover, possible subtle changes of the stoichiometry of the coating surface after thermal treatments are investigated. Finally, the variation of refractive index or coating thickness after heating under different conditions is critically studied.

Experimental

The investigated samples consisted in amorphous Ta 2 O 5 films deposited on crystalline Si wafers. They were produced at Laboratoire de Matériaux Avancés in Lyon by means of IBS in a VEECO Spector system. All coatings have been produced using accelerated, neutralized argon ions as sputtering particles. Argon was fed into the ion-beam sources while oxygen was fed into the chamber, for a total residual pressure inside the chamber of the order of 10 -4 mbar. During deposition, the sputtered coating particles impinged on substrates heated up to about 100 • C; the ion energy and current were of the order of 1 keV and 0.5 A, respectively. For the present study, four different thicknesses of the coatings were produced under the same conditions (4.5, 14, 103 and 953 nm) and a choice among them was performed for various spectroscopic studies, according to the desired information to be obtained.

Concomitant mass spectrometry and thermogravimetric (TGA) measurements were performed in a vacuum of the order of 10 -5 mbar in a Setaram Setsys Evolution 1200 system coupled to a Pfeiffer Vacuum quadrupole mass spectrometer, used to detect the possible evolution of argon gas. [START_REF] Palumbo | Decomposition of NH 3 BH 3 at sub-ambient pressures: a combined thermogravimetry-differential thermal analysis-mass spectrometry study[END_REF] The choice to perform measurements in a vacuum was motivated by the much higher sensitivity of the mass spectrometer when working in a pressure lower than 10 -2 mbar. Measurements were conducted between room temperature and 600 • C with a temperature rate of ≈ 3 • C/min on the "as deposited" samples. A specific sample holder was realized by means of a platinum wire able to sustain the large initial mass of samples, which was much higher than the usual one for TGA experiments (5 g in the present case vs. 10 mg).

Surface morphology was investigated ex situ by a VEECO multimode (Nanoscope IIIa) Atomic Force Microscopy (AFM) in the tapping mode by using nonconductive Si tips with a nominal radius of about 5-7 nm.

The X Ray Photoelectron Spectroscopy (XPS) experiments were carried out at the Nanostructures at Surfaces laboratory at Sapienza University in Rome, in an Ultra-High-Vacuum (UHV) chamber with a base pressure in the low 10 -10 mbar range. An X-ray photon source with hν=1253.6 eV was used (PSP TA10 Mg Kα); photoelectrons were measured with a hemispherical VG Microtech Clam-2 electron analyzer in constant pass energy mode set at 50 eV and with overall energy resolution better than 1 eV. The binding energy (BE) scale was calibrated on a freshly sputtered gold foil in electrical contact with the sample, by setting the Au 4f 7/2 core-level at 84.0 eV BE. XPS experiments were performed on the "as-deposited" samples and on small parts of the same samples heated either in vacuum at 300 • C in situ or in air at 500 • C for 10 hours (ex situ).

Mid-infrared reflectivity measurements were performed at ambient conditions using a commercial Bruker Interferometer coupled with an infrared microscope (Bruker Hyperion 2000) and a gold mirror as reference. For each sample, 25 different positions were considered, collecting spectra on an area of 5 × 5 mm 2 with steps of 1 mm in the two perpendicular directions. The incident light spot size was fixed at 200 μm. Spectra were collected with 4 cm -1 resolution and each spectrum was the average of 256 scans for each point. A nitrogen cooled mercurycadmium-telluride (MCT) and a potassium bromide (KBr) crystal were used as detector and beamplitter, respectively.

Raman measurements have been performed at ambient conditions using a Horiba spectrometer equipped with a 532 nm laser, a 100x objective, a grating of 1800 grooves/mm and a spot size of 2 μm 2 . Each spectrum was the average of 12 acquisitions, with an exposition time of 300 s for each of them. Three spectra have been collected for each A,B,C). Since the Raman intensity can be significantly affected by small changes in the laser focus and in the sample orientation, all the spectra were rescaled and the analysis was limited to the phonon-band profiles.

Infrared and Raman spectroscopy measurements were performed on the samples with a thickness of the coating of 953 nm. As reference, the "as deposited" sample was considered and initial spectra in these conditions were acquired. Parts of the same specimen were afterwards heated up, either in air up to 500 • C for 10 hours or in a vacuum of the order of 10 -5 mbar up to 500 • C for 15 hours. The spectra of samples after annealing in both conditions were further acquired, with a special care to measure in the same positions initially investigated.

Results

A preliminary investigation of the roughness of the surface of the samples was conducted by atomic force microscopy (AFM) topographies on the samples with thickness 4.5 and 103 nm. As reported in Fig. 1, the two Ta 2 O 5 coatings exhibit an extremely flat surface with an r.m.s. roughness in the best case as low as 0.15 nm. For the 4.5 nm thick layer we observed the occurrence of some linear ditches on the surface, not observed instead in the 103 nm thick layer. The presence of these grooves is confirmed by the detection of a low Si oxide signal in XPS (see subsequent section).

In order to investigate whether some residual Argon atoms from the deposition are left in the whole coating layer or only at the interface between the substrate and the coating, we performed mass spectrometry measurements. In particular, the samples were heated in a vacuum of the order of 10 of 14 nm (Fig. 2). This fact can be interpreted as the evidence that argon is present in the samples and it can be released in vacuum with a maximum of effusion around 300 • C; at 500 • C all the argon gas present in the specimen is discharged by means of a thermal treatment in vacuum. The negligible intensity observed for the 14 nm thick sample suggests that argon is present in the whole body of the coating, not only at the interface between Ta 2 O 5 and Si. In fact, if argon came only from the interface, the intensity of the argon peak for the two thicknesses would be the same (as the surface areas of the measured samples with different thickness practically coincide). On the contrary, under the hypothesis that Ar resides in the whole coating, the intensity of the peak should scale with the ratio of the thicknesses (≈70, in the present case) and, therefore, with the sensitivity of the used mass spectrometer it should not be detectable when t=14 nm. In conclusions the mass spectrometry experiments provide evidence that argon is present in the whole amorphous Ta 2 O 5 coating. further study the presence of Ar in the samples and its evolution with thermal treatments, X-ray photoelectron spectroscopy analysis was carried out on the samples with Ta 2 O 5 coverage of 4.5 and 103 nm, respectively, both on "as-deposited" sample and after annealing; the latter either in air at 500 • C (ex-situ) or in a vacuum at 300 • C (in situ). The thin-coated sample was chosen in order to investigate also the Ta 2 O 5 /substrate interface, thanks to the surface sensitivity of the XPS technique [START_REF] Paolone | Thermal stability and reduction of iron oxide nanowires at moderate temperatures[END_REF], to be compared with the thicker Ta 2 O 5 sample with the deeply buried interface. The as-received sample with t=4.5 nm exhibits the characteristic XPS peaks of tantalum oxide with the expected stoichiometry (Fig. 3) and minor contributions from Ar and from Si coming from the near-surface interface (Fig. 3 and Table 1).

The Argon detection on Ta 2 O 5 is not trivial since the Ar 2p peak overlaps with the Ta 4d doublet and the Ar 3s with the Ta 4f (see Table 1). The best Ar detectable peak is the 2s that exhibits a low sensitivity [START_REF] Yeh | Atomic subshell photoionization cross sections and asymmetry parameters: 1 ⩽ Z ⩽ 103[END_REF]. Nevertheless, by means of Ar 2s detection, and very long integration time, we were able to measure the presence of Ar in the Ta 2 O 5 layer, as observed in Fig. 4.

The XPS spectra of the Ar 2s peak clearly point out the presence of a low percentage of Ar in the layers, estimated to less than 0.3 at.%. This small quantity remarkably reduces (by 50% for the 103 nm thick layer) or even disappears (for the 4.5 nm thick layer) after an annealing of 1 hour at 300 • C in UHV conditions. In the samples annealed in air at 500 • C the Ar 2s component (Fig. 5) appears reduced in intensity by only 50% with respect to the untreated samples, suggesting that the treatment in air is less effective in removing the trapped argon.

The occurrence of Ar in both samples rules out the possibility of an argon accumulation occurring only at the Ta 2 O 5 /Si interface and it is in agreement with indications provided by mass spectrometry. This finding is further strengthened by the behavior of Ar signal reduction under sample annealing: while argon is completed desorbed from the thin layer, we have only a partial desorption in the thick layer due to the volumetric permeation.

The thermal treatments also influence the surface stoichiometry of the thin films. We report in Fig. 6 the comparison of the Ta 4f and O 1s regions of Ta-oxides-based film in the as-received state and after the annealing in vacuum at 500 • C. After the thermal treatment, different spectral components emerge in the low binding energy side of the main Ta 2 O 5 -associated peaks, thus suggesting a different oxidation state for tantalum, induced by heating. Specifically, on the basis of previous results [START_REF] Atanassova | XPS study of N 2 annealing effect on thermal Ta 2 O 5 layers on Si[END_REF][START_REF] Perez | Influence of post-deposition annealing on the chemical states of crystalline tantalum pentoxide films[END_REF], one can attribute these smaller components to the formation of different sub-stoichiometric Ta-oxides phases, namely Ta 2 O 3 , TaO, TaO 2 , in agreement with a previous investigation performed by XPS on Ta 2 O 5 thermally treated at 600 • C [START_REF] Prato | Multitechnique investigation of Ta 2 O 5 films on SiO 2 substrates: comparison of optical, chemical andmorphological properties[END_REF]. It must be noted, however, that the change of stoichiometry is likely to be confined to the first few nm of the coating. Extensive changes in bulk stoichiometry would result in measurable changes in the Raman signal. We will show below that this is not observed.

Table 2

Maxima of the interference fringes and their period in the two wafers, before and after annealing either in air (W1) or in vacuum (W2). The effects of thermal treatments on the coatings were investigated also by means of infrared and Raman spectroscopy on the samples with a coating thickness of 952 nm, in order to gain sensitivity to the coating film properties. Figs. 7 and8 report the comparison of the infrared spectra acquired on the "as-deposited" sample and on the same sample after annealing in air (Fig. 7) or in vacuum (Fig. 8) at 500 • C. In all spectra, the feature at 2300 cm -1 is an artefact due to atmospheric absorption.

In the infrared reflectance spectra, one can observe well defined interferences fringes, due to the superposition of the electric fields reflected at the air-Ta 2 O 5 and at the Ta 2 O 5 -Si interfaces, respectively. The period of these fringes is given by (nd) -1 , being n the real refractive index of Ta 2 O 5 and d the optical path of the radiation within the Ta 2 O 5 film. Table 2 displays the period of the fringes calculated for each sample before and after annealing either in air (W1) or in vacuum (W2).

As shown in Table 2, before the annealing, the spacing of the fringes, Δ_max, is almost the same for the two samples, W1 and W2 (2518 and 2511 cm -1 , respectively); after the annealing, the fringe spacing Δ_max reduces for both samples, although the decrease is significantly larger in the sample annealed in vacuum (63 cm -1 ) than in the one thermally treated in air (42 cm -1 ). We can conclude that in both wafers the annealing affects the dn product, but in the sample annealed in vacuum this effect is more evident. Infrared spectroscopy cannot separate the contribution from the refractive index or the coating thickness to the modulation of the fringes. However, previous measurements on amorphous Ta 2 O 5 deposited on SiO 2 substrates indicated that n lightly decreases after thermal treatments [START_REF] Granata | Amorphous optical coatings of present gravitational-wave interferometers[END_REF][START_REF] Anghinolfi | Optical properties of uniform, porous, amorphous Ta 2 O 5 coatings on silica: temperature effects[END_REF]. Therefore, in the present case, the thickness of the coating films should increase after annealing by ~15-20 nm in the two samples (~1-2%). Such increase of the coating thickness was also suggested by previous spectroscopic ellipsometry measurements on Ta 2 O 5 films deposited on silica substrates [START_REF] Anghinolfi | Optical properties of uniform, porous, amorphous Ta 2 O 5 coatings on silica: temperature effects[END_REF] and interpreted as due to an increase of the volume of the pores present in the structure upon annealing [START_REF] Anghinolfi | Optical properties of uniform, porous, amorphous Ta 2 O 5 coatings on silica: temperature effects[END_REF].

It must be noted, however, that no major changes of the macroscopic structure of Ta 2 O 5 occur during the annealing, either in vacuum or in air. In fact, Figs. 9 and 10 display the comparison of the Raman spectra of samples W1 and W2 (the same measured by infrared spectroscopy) before and after the annealing in air or in vacuum. The two peaks at ~620 and ~640 cm -1 can be attributed to silicon oxide, while the underlying band between 600 and 750 cm -1 has been assigned to amorphous Ta 2 O 5 by comparison with the literature [27]. It can be noted (Figs. 9 and 10) that not only the oxide band-profile does not depend on the position on the sample (Figs. 9a,9d,10a,10d), but it remains also unchanged before and after the annealing procedures (Figs. 9c,10c), thus ensuring a good homogeneity of the sample on the micrometer scale.

Conclusions

The present study investigated the effects of thermal treatments in vacuum or in air on amorphous Ta 2 O 5 deposited on crystalline Si, mimicking the mirrors of gravitational-wave detectors. Evidence is provided that the synthesis by Ion Beam Sputtering leaves some argon atoms (less than 0.3 at%) inside the whole thickness of the Ta 2 O 5 film. Argon is released when the sample is heated in vacuum above 200 • C, however it remains entrapped in the sample if the annealing is conducted in air, even at temperatures as high as 500 • C. The macroscopic structure of the samples is not strongly affected by heating; however, some substoichiometric Ta oxides develop on the surface of the coating after heating in vacuum. Finally, the thickness of the Ta 2 O 5 seems to increase by ~1-2% after annealing, due to the increase of the size of the pores present in the coating.

Fig. 1 .

 1 Fig. 1. AFM topographies of Ta 2 O 5 layers with 4.5 nm (left panel) and 103 nm (right panel) coatings.

Fig. 3 .

 3 Fig. 3. XPS survey spectrum of Ta 2 O 5 coating with a thickness of 4.5 nm deposited on crystalline Si.

Fig. 2 .

 2 Fig. 2. Comparison of the signal of the mass spectrometer for argon gas as a function of temperature for the samples with a thickness of 14 and 953 nm.

Fig. 4 .

 4 Fig. 4. XPS spectra of Ar 2s for Ta 2 O 5 layers with thickness 4.5 nm (a) and 103 nm (b). Top spectra are related to the "as deposited" samples, while bottom spectra to samples after one hour annealing at 300 • C in ultra-high vacuum.

Table 1 5 TaFig. 5 .

 155 Fig. 5. XPS spectra of Ar 2s for Ta 2 O 5 layers with thickness 103 nm, after annealing in air at 500 • C.

Fig. 6 .

 6 Fig. 6. XPS spectra of (a) Ta 4f, and (b) O 1s regions of Ta-oxides-based film in the as-received state and after the annealing treatment in vacuum.

Fig. 7 .

 7 Fig. 7. (a) MIR map of "as deposited" sample before the annealing; the 25 spectra of the map are reported with different shades of blue and are shifted along the z axis. (b) Comparison between the average spectra of the sample before (blue) and after (pink) the annealing in air at 500 • C; the vertical lines indicate the maxima of the curves. (c) MIR map after the annealing; the 25 spectra of the map are reported with different shades of pink and are shifted along the z axis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Fig. 8 .

 8 Fig. 8. (a) MIR map in the "as-deposited" state; the 25 spectra of the map are reported with different shades of green and are shifted along the z axis. (b) Comparison between the average spectra before (green) and after (purple) the annealing in vacuum at 500 • C; the vertical lines indicate the maxima of the curves. (c) MIR map after the annealing; the 25 spectra of the map are reported with different shades of purple and are shifted along the z axis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Fig. 9 .

 9 Fig. 9. (a) Raman spectra collected in three different point on w1 (A, B, C) before the annealing; the 3 spectra are reported with different shades of blue and are shifted along the z axis. (b) Comparison between the average spectra of W1 before (blue) and after (pink) the annealing in air at 500 • C. (c) Raman spectra collected in three different point on W1 (A, B, C) after the annealing; the 3 spectra are reported with different shades of pink and are shifted along the z axis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Fig. 10 .

 10 Fig. 10. (a) Raman spectra collected in three different point on W2 (A, B, C) before the annealing; the 3 spectra are reported with different shades of green and are shifted along the z axis. (b) Comparison between the average spectra of W2 before (green) and after (purple) the annealing in vacuum. (c) Raman spectra collected in three different point on W2 (A, B, C) afterthe annealing; the 3 spectra are reported with different shades of purple and are shifted along the z axis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Acknowledgments

The present work has been performed in the framework of the Virgo Coating Research and Development collaboration. The support of Sapienza Grande Progetto di Ateneo 2017, Amaldi Research Center

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

(MIUR program "Dipartimento di Eccellenza" CUP:B81I18001170001), of European Gravitational Observatory and of Istituto Nazionale di Fisica Nucleare is gratefully acknowledged. We wish to thank S. Brutti for the use of one of his furnaces for some thermal treatments.