Dorian Dupuy

Adrien Toutant
email: adrien.toutant@univ-perp.fr

Françoise Bataille

Study of the large-eddy simulation subgrid terms of a low Mach number anisothermal channel flow

The subgrid terms of the low Mach number equations are investigated in a strongly anisothermal low Mach number flow. The filtered low Mach number equations are established in three formulations in order to compare the unweighted classical filter and the density-weighted Favre filter on the one hand, and the filtering of the momentum conservation equation and the velocity transport equation on the other hand. In the three formulations, we establish the filtered equations of mass conservation, momentum conservation, energy conservation, of the ideal gas law and of the resolved kinetic energy transport equation. The magnitude of all subgrid terms is assessed a priori in the three formulations using the results of direct numerical simulations of a strongly anisothermal fully developed turbulent channel flow. The classification of the subgrid terms gives the relevance of various effects of the temperature gradient.

Introduction

In solar power towers, the solar flux is concentrated towards the solar receiver, wherethrough its energy is transferred to the heat transfer fluid. The optimisation of the internal geometry of the solar receiver is a key challenge for the efficiency of solar power towers. The thermal exchange towards the fluid should be maximised while the pressure loss should be minimised. This long-term goal would benefit from accurate numerical simulations of the low Mach number strongly anisothermal turbulent flows found in solar receivers [START_REF] Serra | High-temperature gradient effect on a turbulent channel flow using thermal large-eddy simulation in physical and spectral spaces[END_REF]. However, high numerical costs prevent the direct numerical simulation (DNS) of all scales of turbulence in the conditions of a real solar receiver. An effective alternative is the thermal large-eddy simulation (LES). The LES resolves the largest scales of the flow and models the effect of the small scales on the large scales. The method requires subgrid-scale models suited to low Mach number strongly anisothermal turbulent flows [START_REF] Aulery | Energy transfer process of anisothermal wall-bounded flows[END_REF][START_REF] Aulery | Spectral analysis of turbulence in anisothermal channel flows[END_REF][START_REF] Dupuy | Turbulence kinetic energy exchanges in flows with highly variable fluid properties[END_REF].

The subgrid terms can be investigated a priori from high-resolution three-dimensional flow fields. A priori tests were carried out with particle image velocimetry measurements [START_REF] Liu | On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet[END_REF][START_REF] Carper | Subfilter-scale fluxes over a surface roughness transition. part ii: A priori study of large-eddy simulation models[END_REF] and DNS results, in incompressible isotropic homogeneous turbulence [START_REF] Clark | Evaluation of subgrid-scale models using an accurately simulated turbulent flow[END_REF][START_REF] Abba | Analysis of subgrid scale models[END_REF], rotating turbulence [START_REF] Lu | A priori tests of one-equation LES modeling of rotating turbulence[END_REF], channel flows [START_REF] Piomelli | Model consistency in large eddy simulation of turbulent channel flows[END_REF][START_REF] Li | A-priori analysis of LES subgrid scale models applied to wall turbulence with pressure gradients[END_REF], in two phase divergence-free flows [START_REF] Toutant | DNS of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for LES two-phase flow modelling[END_REF], in passive and active scalar decaying homogeneous turbulence [START_REF] Chumakov | Dynamic structure subgrid-scale models for large eddy simulation[END_REF][START_REF] Ghaisas | A priori evaluation of large eddy simulation subgridscale scalar flux models in isotropic passive-scalar and anisotropic buoyancy-driven homogeneous turbulence[END_REF] and in flows with purely compressible effects, in a temporal shear layer [START_REF] Vreman | A priori tests of large eddy simulation of the compressible plane mixing layer[END_REF][START_REF] Vreman | Subgrid-modelling in LES of compressible flow[END_REF][START_REF] Vreman | Direct and large-eddy simulation of the compressible turbulent mixing layer[END_REF] and in freely decaying homogeneous isotropic turbulence [START_REF] Pino Martín | Subgrid-scale models for compressible large-eddy simulations[END_REF]. In flows with variable fluid properties, the analysis can be carried out using two large-eddy simulation filters, the unweighted classical filter and the density-weighted Favre filter [START_REF] Favre | The equations of compressible turbulent gases[END_REF]. The Favre filter is used by most authors in compressible flows [START_REF] Garnier | Large eddy simulation for compressible flows[END_REF]. The classical filter has been employed by Yoshizawa [START_REF] Yoshizawa | Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling[END_REF], Sun and Lu [START_REF] Sun | A large eddy simulation approach of compressible turbulent flow without density weighting[END_REF], Boersma and Lele [START_REF] Boersma | Large eddy simulation of compressible turbulent jets[END_REF], Bodony and Lele [START_REF] Bodony | On using large-eddy simulation for the prediction of noise from cold and heated turbulent jets[END_REF] and Sidharth et al. [START_REF] Sidharth | Filtered velocity based LES of mixing in high speed recirculating shear flow[END_REF]. In addition, the set of subgrid terms to model depends on the formulation of the filtered governing equations, of the energy equation in particular but also of the momentum equation [START_REF] Sidharth | Stretched-vortex based subgrid-scale modeling for variable-density flows[END_REF]. Using a priori tests, Vreman et al. [START_REF] Vreman | A priori tests of large eddy simulation of the compressible plane mixing layer[END_REF] and Martín et al. [START_REF] Pino Martín | Subgrid-scale models for compressible large-eddy simulations[END_REF] assessed the amplitude of all subgrid terms involved in the compressible Navier-Stokes equations for different formulations of the energy conservation equation.

In the literature, the subgrid terms have not been investigated for low Mach number strongly anisothermal turbulent flows. In this paper, we extend the analysis to this configuration. In strongly anisothermal flows, the large temperature differences require taking into account the variations of the fluid properties. If in addition the flow is turbulent, the variations of the fluid properties are related to the velocity fluctuations. The figure 2 is an illustration of the correlation between the velocity and temperature. The isosurface of figure 2 is not only an isotherm but also a surface of isodensity, isoviscosity and isoconductivity. Both the velocity and temperature profiles are turbulent and exhibit a wide range of length scales. The correlations between velocity and temperature lead to additional subgrid terms associated with the nonlinearities of the viscous term and of the heat flux. In the literature, these additional subgrid terms are always neglected (see for example [START_REF] Serra | High-temperature gradient effect on a turbulent channel flow using thermal large-eddy simulation in physical and spectral spaces[END_REF][START_REF] Yahya | Turbulent forced convective flow in an anisothermal channel[END_REF][START_REF] Al-Zurfi | LES of rotational effects on film cooling effectiveness and heat transfer coefficient in a gas turbine blade with one row of air film injection[END_REF][START_REF] Cheung | A fully-coupled simulation of vortical structures in a large-scale buoyant pool fire[END_REF][START_REF] Serra | Thermal large eddy simulation in a very simplified geometry of a solar receiver[END_REF][START_REF] Sanchez | Large eddy simulation of thermal boundary layer spatial development in a turbulent channel flow[END_REF][START_REF] Li | Large eddy simulation of compound angle hole film cooling with hole length-to-diameter ratio and internal crossflow orientation effects[END_REF]). However, to the best of our knowledge, there is no study to justify this assumption in low Mach number flows. We investigate whether these additional subgrid terms may be neglected. This study is essential for the flows found in solar receivers, characterised by a strong coupling between turbulence and temperature, along with high variations of the fluid properties (density, viscosity and thermal conductivity) with temperature [START_REF] Toutant | Turbulence statistics in a fully developed channel flow submitted to a high temperature gradient[END_REF].

We study the subgrid terms using the results of direct numerical simulations of a strongly anisothermal fully developed turbulent channel flow. The investigation is based on a particular form of the Navier-Stokes equations under the low Mach number hypothesis, called low Mach number equations [START_REF] Dupuy | Study of the sub-grid terms of the large-eddy simulation of a low Mach strongly anisothermal channel flow[END_REF]. The filtering of the low Mach number equations gives rise to specific subgrid terms. Three formulations of the filtered low Mach number equations are investigated, which leads to three specific sets of subgrid terms. The analysis is carried out using the classical filter and the Favre filter. With the classical filter, we compare the filtering of the momentum conservation equation to the filtering of the velocity transport equation. We assess the magnitude of each term in the three formulations. This provides insights on the asymmetry of the subgrid and non-subgrid terms between the two sides of the channel and on the relative importance of the various competing physical processes. The study takes into account both the subgrid terms associated with non-commutation of the filter with the derivative and the nonlinearities arising from the large variations of the fluid properties. The objective is to identify the subgrid terms that can be neglected and the subgrid terms that should be modelled in the three formulations.

We describe the low Mach number equations in section 2. In section 3, we study the filtering of the low Mach number equations in three formulations. The subgrid terms derived in the general case are then estimated in the channel flow configuration. The detailed channel flow configuration can be found in section 4 as well as the numerical method used to compute the subgrid terms. The results are discussed in section 5.

Low Mach number equations

We consider a turbulent flow at low Mach number (M a < 0.3). The low Mach number hypothesis let us use Paolucci's method [START_REF] Paolucci | On the filtering of sound from the Navier-Stokes equations[END_REF] to remove acoustic effects from the Navier-Stokes equations. Each variable of the Navier-Stokes equations is written as a power series of the squared Mach number. Keeping only the smaller-order terms leads to the so-called low Mach number equations. The pressure is split in two parts: the thermodynamical pressure P , mean pressure in the domain, and the mechanical pressure P 0 , associated with momentum variations. The thermodynamical pressure is constant in space. The gas is air, an ideal gas for the purpose of this study. The effects of gravity are neglected.

Those considerations lead to the low Mach number equations, given by:

• Mass conservation equation ∂ρ ∂t + ∂ρU j ∂x j = 0, (1)
• Momentum conservation equation

∂ρU i ∂t = - ∂ρU j U i ∂x j - ∂P ∂x i + ∂Σ ij ∂x j , (2)
• Energy conservation equation

∂U j ∂x j = - 1 γP 0 (γ -1) ∂Q j ∂x j - ∂P 0 ∂t , (3)
• Ideal gas law

T = P 0 ρr , (4)
with ρ the density, T the temperature, Σ ij the shear-stress tensor, Q j the conductive heat flux, γ the heat capacity ratio, r the ideal gas specific constant, t the time, P the mechanical pressure, P 0 the thermodynamical pressure, U i the i-th component of velocity and x i the Cartesian coordinate in i-th direction. Einstein summation convention is used and δ ij is the Kronecker delta.

In the low Mach number equations, the energy conservation equation has a distinctive form. The local conservation of energy is imposed by a constraint on the divergence of the velocity. Namely, the dilatation of the fluid is proportional to the difference between the conductive heat flux and the global variation of thermodynamical pressure. To obtain equation (3), we use the low Mach number hypothesis to approximate the compressible energy conservation equation in enthalpy form as [START_REF] Nicoud | Conservative high-order finite-difference schemes for low-Mach number flows[END_REF] ∂ρC p T ∂t

+ ∂ρU j C p T ∂x j = ∂P 0 ∂t - ∂Q j ∂x j , (5)
where C p is the heat capacity at constant pressure. We then use the ideal gas law (4) to substitute in this equation ρT with P 0 /r. Using the fact that the thermodynamical pressure is constant in space, we isolate the divergence of the velocity and combine the two temporal terms.

We assume that air is a Newtonian fluid to compute the shear-stress tensor,

Σ ij = µ(T ) ∂U i ∂x j + ∂U j ∂x i - 2 3 µ(T ) ∂U k ∂x k δ ij . (6)
The heat flux is given by

Q j = -λ(T ) ∂T ∂x j . (7)
The variations of viscosity with temperature are accounted for by the Sutherland's law [START_REF] Sutherland | The viscosity of gases and molecular force[END_REF],

µ(T ) = µ 0 T T 0 3 2 T 0 + S T + S , (8)
with µ 0 = 1.716 • 10 -5 , S = 110.4 and T 0 = 273.15 K. The thermal conductivity is deduced from the Prandlt number P r and the heat capacity at constant pressure C p , both assumed constant with P r = 0.76 and C p = 1005 J kg -1 K -1 . The ideal gas specific constant is r = 287 J kg -1 K -1 .

We will study the filtering of the low Mach number equations in three formulations that we shall call the Classical formulation, the Favre formulation and the Velocity formulation.

Filtering of the low Mach number equations

The large-eddy simulation is based on the idea of scale separation. Theoretically, the separation is carried out by the application of a filter, denoted ( • ), on the Navier-Stokes equations. We restrict our discussion in this paper to a spatial filter. The filter is taken to verify the properties of conservation of constants, a = a with a a constant, and of linearity, φ + ψ = φ + ψ for any φ and ψ [START_REF] Sagaut | Large eddy simulation for incompressible flows: an introduction[END_REF]. Note however that the filter may be inhomogeneous and thus not commute with derivation. In the following, this spatial filter will be referred to as the classical filter. The formulation of the filtered low Mach number equations involves two kinds of subgrid terms. The first kind arises from the non-commutation of the filter with the derivative and are related to the variations of filter width, that is from the LES mesh inhomogeneity. They are denoted in the form C β α . The second kind arises from nonlinearities, that is from the non-commutation between the filtering and the multiplication operator. They are denoted in the form F α . Most of them are related to the large variations of the fluid properties, with the notable exception of the subgrid term associated with momentum convection.

Classical formulation

In the Classical formulation, the low Mach number equations are filtered with the classical filter and expressed in terms of classical-filtered variables. The filtered low Mach number equations in the Classical formulation are given by:

• Mass conservation equation ∂ρ ∂t + ∂ ∂x j ρU j + F ρU j + C j ρU j = 0, (9)
• Momentum conservation equation

∂ ∂t ρU i + F ρU i = - ∂ ∂x j ρU j U i + F ρU j U i -C j ρU j U i - ∂P ∂x i -C i P + ∂ ∂x j Σij + F Σ ij + C j Σ ij , (10)
• Energy conservation equation

∂U j ∂x j + C j U j = - γ -1 γP 0 ∂ ∂x j Qj + F Q j + C j Q j - 1 γP 0 ∂P 0 ∂t , (11)
• Ideal gas law

T = P 0 r 1 ρ + F 1/ρ , (12)
with the classical filter counterparts of the shear-stress tensor and of the heat flux given by

Σij = µ(T ) ∂U i ∂x j + ∂U j ∂x i - 2 3 µ(T ) ∂U k ∂x k δ ij , (13)
Qj = -λ(T ) ∂T ∂x j , (14)
with µ the dynamic viscosity, λ the thermal conductivity and the following subgrid terms:

C j ρU j = ∂ρU j ∂x j - ∂ρU j ∂x j (15) C j ρU j U i = ∂ρU j U i ∂x j - ∂ρU j U i ∂x j (16)
C i P = ∂P ∂x i - ∂P ∂x i (17)
C i Σ ij = ∂Σ ij ∂x j - ∂Σ ij ∂x j (18)
C j U j = ∂U j ∂x j - ∂U j ∂x j (19)
C j Q j = ∂Q j ∂x j - ∂Q j ∂x j (20)
F ρU j = ρU j -ρU j (21)

F ρU i = ρU i -ρU i (22)
F ρU j U i = ρU j U i -ρU j U i (23)
F Σ ij = Σ ij -Σij (24)
F Q j = Q j -Qj (25)
F 1/ρ = 1 ρ - 1 ρ (26)
In addition, it is useful to express the transport equation of the resolved kinetic energy ρE = 1 2 ρU i U i . Indeed, the total energy conservation is not explicitly stated in the system of equations ( 9)-( 12) but implied by the momentum conservation equation [START_REF] Carper | Subfilter-scale fluxes over a surface roughness transition. part ii: A priori study of large-eddy simulation models[END_REF] and the energy conservation equation [START_REF] Cheung | A fully-coupled simulation of vortical structures in a large-scale buoyant pool fire[END_REF]. The resolved kinetic energy transport equation is the part of total energy conservation related to the momentum conservation. It is obtained from equation [START_REF] Carper | Subfilter-scale fluxes over a surface roughness transition. part ii: A priori study of large-eddy simulation models[END_REF] multiplied by U i and equation [START_REF] Calvin | An object-oriented approach to the design of fluid mechanics software[END_REF],

∂ρE ∂t - ∂EF ρU j ∂x j + F ρU j ∂E ∂x j -EC j ρU j + ∂U i F ρU i ∂t -F ρU i ∂U i ∂t = - ∂ ∂x j ρU j E + U i F ρU j U i + F ρU j U i ∂U i ∂x j -U i C j ρU j U i - ∂U i P ∂x i + P ∂U i ∂x i -U i C i P + ∂ ∂x j U i Σij + U i F Σ ij -Σij + F Σ ij ∂U i ∂x j + U i C j Σ ij . (27)
This equation gives the contribution of the subgrid terms of the filtered momentum conservation equation to the balance of the resolved kinetic energy. It will be used to assess the magnitude of their energetic contribution. The contribution of the subgrid terms with regard to total energy is given by the subgrid terms of the resolved kinetic energy transport equation [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re τ = 590[END_REF] and the energy conservation equation [START_REF] Cheung | A fully-coupled simulation of vortical structures in a large-scale buoyant pool fire[END_REF].

Favre formulation

The use of the Favre filter ( • ) is common when working with variable density flows. It is a variable change in which filtered variables are weighted by density. For any variable φ, the Favre-filtered variable φ is defined as φ = ρφ/ρ. With the Favre filter, we avoid the subgrid terms associated with the nonlinearities of the form ρφ in the convective term of the mass conservation equation, the time derivative term of the momentum conservation equation and the ideal gas law. On the other hand, a subgrid term is added when a variable φ appears without the density, as the velocity in the energy conservation equation. The subgrid terms associated with the nonlinearities of the shear-stress tensor and of the conductive heat flux are modified by the use of the Favre filter, but it is not obvious to what extent this affects the behaviour and importance of the subgrid terms. Finally, the subgrid terms associated with the non-commutation of the filter with the derivative are not modified by the use of the Favre filter.

In the Favre formulation, the low Mach number equations are filtered with the classical filter and expressed in terms of Favre-filtered variables. The filtered low Mach number equations in the Favre formulation are given by:

• Mass conservation equation ∂ρ ∂t + ∂ρ U j ∂x j + C j ρU j = 0, (28)
• Momentum conservation equation

∂ρ U i ∂t = - ∂ ∂x j ρ U j U i + ρG U j U i -C j ρU j U i - ∂P ∂x i -C i P + ∂ ∂x j Σ ij + G Σ ij + C j Σ ij , (29)
• Energy conservation equation

∂ ∂x j U j + ρG U j /ρ + C j U j = - γ -1 γP 0 ∂ ∂x j Q j + G Q j + C j Q j - 1 γP 0 ∂P 0 ∂t , (30)
• Ideal gas law

T = P 0 ρr , (31)
with the Favre filter counterparts of the shear-stress tensor and of the heat flux given by

Σ ij = µ( T ) ∂ U i ∂x j + ∂ U j ∂x i - 2 3 µ( T ) ∂ U k ∂x k δ ij , (32)
Q j = -λ( T ) ∂ T ∂x j , (33)
and the following subgrid terms specific to the Favre formulation:

G U j U i = U j U i -U j U i (34)
G Σ ij = Σ ij -Σ ij (35)
G U j /ρ = U j /ρ -U j /ρ (36)
G Q j = Q j -Q j ( 37
)
The subgrid terms F ρU j and G U j /ρ are closely related,

F ρU j ρ = -ρG U j /ρ . (38)
These two subgrid terms express explicitly the correlation between density and velocity.

The resolved kinetic energy transport equation is obtained from equation ( 29) multiplied by U i and equation [START_REF] Nicoud | Conservative high-order finite-difference schemes for low-Mach number flows[END_REF],

∂ρE ∂t -E C j ρU j = - ∂ ∂x j ρ U j E + ρ U i G U j U i + ρG U j U i ∂ U i ∂x j -U i C j ρU j U i - ∂ U i P ∂x i + P ∂ U i ∂x j -U i C i P + ∂ ∂x j U i Σ ij + U i G Σ ij -Σ ij + G Σ ij ∂ U i ∂x j + U i C j Σ ij . (39)
Compared to the Classical formulation, the Favre formulation is a priori more suited to the momentum conservation equation, the mass conservation equation and the ideal gas law but is less appropriate for the energy conservation equation.

Velocity formulation

The Velocity formulation is based on the velocity filtering approach suggested by Sidharth and Candler [START_REF] Sidharth | Stretched-vortex based subgrid-scale modeling for variable-density flows[END_REF], Sidharth et al. [START_REF] Sidharth | Filtered velocity based LES of mixing in high speed recirculating shear flow[END_REF]. The momentum conservation equation in the low Mach number equations is rewritten before filtering as the velocity transport equation,

∂U i ∂t = - ∂U j U i ∂x j + U i ∂U j ∂x j - 1 ρ ∂P ∂x i + 1 ρ ∂Σ ij ∂x j . (40)
The equations are then filtered with the classical filter and expressed in terms of classical-filtered variables. The filtered low Mach number equations in the Velocity formulation are given by the mass conservation equation ( 9), the energy conservation equation [START_REF] Cheung | A fully-coupled simulation of vortical structures in a large-scale buoyant pool fire[END_REF], the ideal gas law [START_REF] Chumakov | Dynamic structure subgrid-scale models for large eddy simulation[END_REF] and the velocity transport equation:

ρ ∂U i ∂t = -ρ ∂ ∂x j U j U i + F U j U i -ρC j U j U i + ρU i ∂U j ∂x j + ρF U i ∂ j U j + ρU i C j U j - ∂P ∂x i -C i P -ρF ∂ i P/ρ + ∂ ∂x j Σij + F Σ ij + C j Σ ij + ρF ∂ j Σ ij /ρ , (41)
with the following subgrid terms:

C j U j U i = ∂U j U i ∂x j - ∂U j U i ∂x j ( 42
)
F U j U i = U j U i -U j U i (43)
F U i ∂ j U j = U i ∂U j ∂x j -U i ∂U j ∂x j (44)
F ∂ i P/ρ = 1 ρ ∂P ∂x i - 1 ρ ∂P ∂x i ( 45
)
F ∂ j Σ ij /ρ = 1 ρ ∂Σ ij ∂x j - 1 ρ ∂Σ ij ∂x j (46)
Upon filtering, the velocity transport equation gives rise to different subgrid terms than the momentum conservation equation. The time derivative of the velocity density product subgrid term F ρU i does not appear in the Velocity formulation. However, there is three additional subgrid terms associated with the divergence of the velocity and the correlation of the density with other terms of the equation.

The resolved kinetic energy transport equation is given by

ρ ∂E ∂t = -ρ ∂ ∂x j U j E + F U j U i U i + ρF U j U i ∂U i ∂x j -ρU i C j U j U i + ρE ∂U j ∂x j + ρU i F U i ∂ j U j + 2ρEC j U j - ∂P U i ∂x i + P ∂U i ∂x i -U i C i P -ρU i F ∂ i P/ρ + ∂ ∂x j Σij U i + F Σ ij U i -Σij + F Σ ij ∂U i ∂x j + U i C j Σ ij + ρU i F ∂ j Σ ij /ρ . ( 47
)
This equation is obtained from equation ( 41) multiplied by U i and equation ( 9).

4 Numerical study configuration

Channel flow configuration

We consider a fully developed three-dimensional anisothermal channel flow, as shown in figure 1. This geometry is one of the simpler that reproduces the distinctive features of low Mach number strongly anisothermal turbulent flows. It is therefore well suited to the study of the subgrid scale specificities of these flows. The channel is periodic in both the streamwise (x) and spanwise (z) directions. The temperatures of the two plane channel walls are T 2 = 586 K at the hot wall and T 1 = 293 K at the cold wall. This creates the temperature gradient in the wall-normal direction (y). Two mean friction Reynolds number are considered, Re τ = 180 and Re τ = 395, where Re τ is defined as the average of the friction Reynolds numbers Re τ,ω calculated at the hot and cold wall,

Re τ,ω = U τ h ν ω , (48)
with U τ = ν ω (∂ y U x ) 0.5 ω the friction velocity and ν ω the wall kinematic viscosity. We show in figure 2 an isotherm and the velocity magnitude on a plane for a given timestep at Re τ = 180.

Numerical settings

To provide the data required to compute the subgrid terms of the low Mach number equations, direct numerical simulations of the fully developed channel flow described in 4.1 are carried for the two selected friction Reynolds number (180 and 395). At Re τ = 180, the domain size is 4πh × 2h × 2πh, and the mesh used contains 384 × 266 × 384 grid points. At Re τ = 395, the domain size is 4πh × 2h × 4/3πh and the mesh used contains 768 × 512 × 512 grid points. In both cases, the mesh is regular in both homogeneous directions and follow a hyperbolic tangent law in the wall-normal coordinate direction. The wall-normal grid coordinates are given by

y k = L y 1 + 1 a tanh k -1 N y -1 -1 tanh -1 (a) , (49)
with a the mesh dilatation parameter and N y the number of grid points in the wall-normal direction. The cell sizes in wall units are ∆ + x = 8.5, ∆ + y = 0.13 -4.2 and ∆ + z = 4.2 at Re τ = 180; ∆ + x = 9.1, ∆ + y = 0.25 -4.1 and ∆ + z = 4.7 at Re τ = 395. The two meshes have the same level of refinement. The small differences are due to the constraints of the numerical method (multigrid solver) and parallelism. A finite volume method is used with a third-order Runge-Kutta time scheme and a fourth-order centred momentum convection scheme. This is performed using the TrioCFD software [START_REF] Calvin | An object-oriented approach to the design of fluid mechanics software[END_REF]. The numerical set-up is validated through a mesh convergence study and by comparison of our results in the incompressible case to the reference data of Moser et al. [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re τ = 590[END_REF] and Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at Re τ = 180[END_REF]. These numerical results have been validated against experimental data [START_REF] Clark | A study of incompressible turbulent boundary layers in channel flow[END_REF][START_REF] Eckelmann | The structure of the viscous sublayer and the adjacent wall region in a turbulent channel flow[END_REF][START_REF] Kreplin | Behavior of the three fluctuating velocity components in the wall region of a turbulent channel flow[END_REF][START_REF] Antonia | Low-reynolds-number effects in a fully developed turbulent channel flow[END_REF]. We provide in figure 3 a comparison of the spatial turbulence kinetic energy terms computed by our numerical procedure to the reference data of Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at Re τ = 180[END_REF] at Re τ = 180. Similarly, the results have been compared to the reference data of Moser et al.

[27] at Re τ = 395 (not shown here). This validates our numerical method at the incompressible limit. In the anisothermal configuration, the same code has been validated against experimental data for a similar friction Reynolds number and temperature gradient by Bellec et al. [START_REF] Bellec | Large eddy simulations of thermal boundary layer developments in a turbulent channel flow under asymmetrical heating[END_REF].

Filtering process

To compute the subgrid terms a priori, we filter explicitly the DNS flow field. At Re τ = 180, three filters of varying width are investigated, from now on called "filter A", "filter B" and "filter C". The three filters are three-dimensional box filter, also known as top-hat filter, of uniform width in both homogeneous directions and nonuniform width in the wall-normal direction. The filter sizes in wall units and the number of grid points corresponding to the filters are,

• filter A: ∆ The filter A corresponds to an average over five DNS cells in the three directions. It is computed using the following discrete approximation of the box filter:

φ(x i , y k , z j , t) = 1 25 (y k +3 -y k -2 ) i+2 i =i-2 k+2 k =k-2 j+2 j =j-2 φ(x i , y k , z j , t) (y k +1 -y k ) .
(50) The variations of the filter width in the wall-normal direction follow those of the DNS mesh.

The filters B and C are constructed to follow a hyperbolic tangent law (49) in the wall-normal direction. We cannot use the same method as for filter A to compute the filters B and C because the filter width is a non-integer multiple of the DNS cell size. In order to carry out the filtering with an arbitrary filter length, the DNS data are first interpolated using a cubic spline. Then, the top-hat filter is computed from the interpolated values without mesh restrictions. The spline interpolation adds an additional filtering to the box filter. However, this additional filter can be neglected given the resolution of the direct numerical simulation. The interpolation and filtering are not computed in three dimensions but sequentially in the three spatial directions. This is possible because the box filter is separable, that is, can be expressed as the product of three one-dimensional filters. Using the filter width of filter A, the method given similar results to equation (50).

The filtered quantities are not computed very close the domain boundary, where not enough points are available to apply the filter. This is justified by the fact that in practice, the large-eddy simulation of the channel would be carried out with DNS precision very close to the wall.

Results and discussion

We assess the importance of the subgrid terms in the Classical formulation, the Favre formulation and the Velocity formulation. This provides valuable data to determine which terms can be neglected and which terms should be modelled. The strength of the subgrid terms are investigated according to their quadratic mean, or root mean square, compared to the non-subgrid terms. Numerically, the quadratic mean is computed by an average on the two homogeneous directions and on 100 non-consecutive time steps that cover all thermodynamic configurations of the flow. The results are converged as the mean computed on 50 time steps is identical to the mean computed on 100 time steps. Since the flow is homogeneous in the streamwise and spanwise directions, the analysis is carried out as a function of the wall-normal coordinate y, scaled by the height of the channel and in the classical wall scaling

y + = Re τ y h = yU τ ν ω . ( 51
)
The subgrid terms are first investigated with the filter A. Then, the effect of the filter width on the results is examined.

Magnitude of the subgrid terms

In this section, the magnitude of subgrid terms is assessed at Re τ = 180 and Re τ = 395 with the filter A. The results at Re τ = 180 and Re τ = 395 are identical with regard to the classification of the subgrid terms. We thus only provide the results with the filter A at Re τ = 180, for the mass conservation equation (figure 4), the streamwise momentum conservation equation or the streamwise velocity transport equation (figure 5), the spanwise momentum conservation equation or the spanwise velocity transport equation (figure 6), the wall-normal momentum conservation equation or the wall-normal velocity transport equation (figure 7), the resolved kinetic energy transport equation (figure 8), the energy conservation equation (figure 9), and the ideal gas law (figure 10). In each case, the results are given in the entire channel with a linear abscissa and at the cold side with a logarithmic abscissa, to emphasise the near-wall region. Only the cold side is shown with the logarithmic abscissa because a large-eddy simulation of the channel would be less resolved at the cold side than at the hot side. Indeed, given the dependence of density, viscosity and conductivity on temperature, the temperature gradient generates an asymmetry between the hot and cold sides of the channel. The dynamic viscosity and the friction velocity are higher at the hot wall and the friction Reynolds number is larger at the cold side. All terms of the mass conservation equation, momentum conservation equation, resolved kinetic energy transport equation and energy conservation equation have a larger amplitude at the cold side than at the hot side. In addition, the local maxima of the profiles of most subgrid terms are closer to the wall at the cold side than at the hot side. This asymmetry can in a large part be explained by the local variations of the mean fluid properties [START_REF] Toutant | Turbulence statistics in a fully developed channel flow submitted to a high temperature gradient[END_REF]. The profiles are also subject to a low Reynolds number effect [START_REF] Dupuy | Turbulence kinetic energy exchanges in flows with highly variable fluid properties[END_REF]. Since the three filters used are symmetric with respect to the centre of the channel, the cold side is less resolved in wall units than the hot side. This contributes to the fact that the subgrid terms have a larger amplitude at the cold side than at the hot side in the mass conservation equation, momentum conservation equation, resolved kinetic energy transport equation and energy conservation equation.

Mass conservation equation

The mass conservation equation (figure 4) gives a relationship between the time derivative of density and the divergence of the mass flux. With the classical filter (Classical or Velocity formulation), the subgrid term associated with the densityvelocity correlation ∂ j F ρU j is very small at the centre of the channel and remains one order of magnitude smaller than the non-subgrid terms near the wall, where it is the most significant. The filter-derivative non-commutation subgrid term C j ρU j has a smaller amplitude in most part of the channel but has the same order of magnitude at the wall. With the Favre filter, it is the only subgrid term in the mass conservation equation.

Since the statistical average of the time derivative of density is zero, the divergence of the resolved mass flux ∂ j ρU j is in balance with the subgrid terms. With the classical filter, the statistical average of the density-velocity correlation subgrid term has the same order of magnitude as the divergence of the resolved mass flux, while the filter-derivative non-commutation subgrid term is negligible. The modelling of the subgrid term ∂ j ρU j is necessary to take this behaviour into account in a large-eddy simulation. Therefore, we consider that with the classical filter, the modelling of the density-velocity correlation subgrid term ∂ j F ρU j is more important than the modelling of the filter-derivative non-commutation subgrid term C j ρU j .

Momentum conservation equation

To study the subgrid terms related to the momentum conservation, we investigate the subgrid terms as they appear in the streamwise, spanwise and wall-normal momentum conservation equations (figures 5, 6 and 7 respectively), or velocity transport equation in the case of the Velocity formulation, and the energetic contribution of the subgrid terms from the resolved kinetic energy transport equation (figure 9). In either case, the time derivative term and the convective term predominate at the centre of the channel. Nevertheless, these two terms in a large part cancel each other out. To be considered negligible, a subgrid term should thus at least be negligible compared to the third largest non-subgrid term, namely the pressure term.

Two non-subgrid terms are negligible throughout the entire channel, namely the velocity-dilatation product in the spanwise and wall-normal velocity transport equations and the pressure-dilatation product in the resolved kinetic energy transport equation. The velocity-dilatation product in the streamwise velocity transport equation is also small but is not negligible. This shows the small influence of dilatation on the flow dynamics in low Mach number strongly anisotermal flows. This is consistent with the negligible effect of dilatation on the turbulence kinetic energy budget [START_REF] Dupuy | Turbulence kinetic energy exchanges in flows with highly variable fluid properties[END_REF].

∂ t ρ -U - x ∂ j ρ -U - j U - x ∂ x P - ∂ j Σ v xj C j ρU j U x C j Σ xj ∂ t F ρU x ∂ j F ρU j U x ∂ j F Σ xj
∂ t ρ -U - x ∂ j ρ -U - j U - x ∂ x P - ∂ j Σ v xj C j ρU j U x C j Σ xj ∂ t F ρU x ∂ j F ρU j U x ∂ j F Σ xj (
∂ t ρ -U ~x ∂ j ρ -U ~j U ~x ∂ x P - ∂ j Σ ^xj C j ρU j U x C j Σ xj ∂ j G ρU j U x ∂ j G Σ xj
∂ j ρ -U ~j U ~x ∂ x P - ∂ j Σ ^xj C j ρU j U x C j Σ xj ∂ j G ρU j U x ∂ j G Σ xj (
ρ -∂ t U - x ρ -∂ j U - j U - x rudu ∂ x P - ∂ j Σ v xj rcuu rucuj C j Σ xj rdfuu rfudu rf1d ∂ j F Σ xj
+ ρ -∂ t U - x ρ -∂ j U - j U - x rudu ∂ x P - ∂ j Σ v xj rcuu rucuj C j Σ xj rdfuu rfudu rf1d ∂ j F Σ xj rf1v (
∂ t ρ -U - z ∂ j ρ -U - j U - z ∂ z P - ∂ j Σ v zj C j ρU j U z C j Σ zj ∂ t F ρU z ∂ j F ρU j U z ∂ j F Σ zj
∂ t ρ -U - z ∂ j ρ -U - j U - z ∂ z P - ∂ j Σ v zj C j ρU U z C j Σ zj ∂ t F ρU z ∂ j F ρU j U z ∂ j F Σ zj (
∂ t ρ -U ~z ∂ j ρ -U ~j U ~z ∂ z P - ∂ j Σ ^zj C j ρU j U z C j Σ zj ∂ j G ρU j U z ∂ j G Σ zj
∂ j ρ -U ~j U ~z ∂ z P - ∂ j Σ ^zj C j ρU j U z C j Σ zj ∂ j G ρU j U z ∂ j G Σ zj (
+ ρ -∂ t U - z ρ -∂ j U - j U - z rudu ∂ z P - ∂ j Σ v zj rcuu rucuj C j Σ zj rdfuu rfudu rf1p ∂ j F Σ zj rf1v (
∂ t ρ -U - y ∂ j ρ -U - j U - y ∂ y P - ∂ j Σ v yj C j ρU j U y C i P C j Σ yj ∂ t F ρU y ∂ j F ρU j U y ∂ j F Σ yj
∂ t ρ -U - y ∂ j ρ -U - j U - y ∂ y P - ∂ j Σ v yj C j ρU j U y C i P C j Σ yj ∂ t F ρU y ∂ j F ρU j U y ∂ j F Σ yj (
∂ t ρ -U ~y ∂ j ρ -U U ~y ∂ y P - ∂ j Σ ^yj C j ρU j U y C i P C j Σ yj ∂ j G ρU j U y ∂ j G Σ yj
∂ j ρ -U ~j U ~y ∂ y P - ∂ j Σ ^yj C j ρU j U y C i P C j Σ yj ∂ j G ρU j U y ∂ j G Σ yj (
In the momentum conservation equation (figures 5, 6 and 7), the most significant subgrid term in the three formulations is the subgrid term associated with momentum or velocity convection, ∂ j F ρU j U i , ∂ j ρG U j U i and ∂ j F U j U i . It is larger than viscous term in the bulk and cannot be neglected at the wall. This subgrid term has a larger amplitude at the centre of the channel in the Classical than in the Favre and Velocity formulations but a slightly lower amplitude near the wall. It is harder to model in the Classical formulation since it is a triple correlation whereas it is a double correlation in the Favre and Velocity formulations.

Three subgrid terms have a medium-sized amplitude. In the Classical formulation, the time derivative subgrid term ∂ t F ρU i is the second most significant subgrid term. In particular, it has the same order of magnitude as the viscous term at the centre of the channel. Present in the three formulations, the filter-derivative noncommutation subgrid term associated with the convective term, C j ρU j U i or C j U j U i , has a lower amplitude than any non-negligible non-subgrid term but by less than one order of magnitude. It appears very small in the spanwise and wall-normal directions, but is more significant in the streamwise momentum equation. In the Velocity formulation, this subgrid term is followed by the velocity divergence filterderivative non-commutation subgrid term C j U j , which has a similar behaviour. Since the subgrid terms C j U j U i and C j U j appear with opposite sign in the velocity transport equation, the subgrid terms associated with the filter-derivative non-commutation are less significant in the Velocity formulation.

The subgrid term associated with the nonlinearities of the viscous diffusion, ∂ j F Σ ij and G Σ ij , can be considered insignificant in the three formulations. Its amplitude is increased by around one order of magnitude in the Favre formulation but remains one order of magnitude smaller than any non-subgrid term in the entire channel. Assuming that the small-scale variations of velocity are significant, this suggests that the influence of small-scale variations of viscosity on the behaviour of the flow is negligible. The filter-derivative non-commutation subgrid term associated with the pressure and viscous terms C i P and C i Σ ij are also negligible. All three of the additional subgrid terms in the Velocity formulation, F U i ∂ j U j , F ∂ i P/ρ and F ∂ j Σ ij /ρ , are also found to have a very small amplitude. The most significant of the three is the subgrid term F ∂ i P/ρ . The Velocity formulation thus appears more interesting than the classical formulation with regard to the subgrid-scale modelling as the time derivative subgrid term is replaced with three negligible subgrid terms.

The resolved kinetic energy transport equation (figure 9) confirms the classification of the subgrid terms in the momentum conservation equation. The most significant subgrid term in the momentum conservation equation have a significant energetic contribution, reaffirming the importance of the subgrid terms associated with momentum or velocity convection, ∂ j F ρU j U i , ∂ j ρG U j U i and ∂ j F U j U i , and, in the Classical formulation, of the subgrid term associated with the time derivative term, ∂ t F ρU i .

Energy conservation equation

Omitting the time derivative of the thermodynamical pressure, constant in space, the energy conservation equation (figure 9) expresses the equality, up to a constant scalar factor, of the divergence of the velocity and of the heat flux. The difference between the divergence of the resolved velocity and heat flux represents the effect of the subgrid terms. With the classical filter (Classical or Velocity formulation), the only significant subgrid term is the filter-derivative non-commutation subgrid term C j U j associated with the divergence of the velocity. The subgrid terms ∂ j F Q j and C j Q j associated with the heat flux are at least 50 times smaller than the non-subgrid terms throughout the entire channel. This is consistent with the assumption that the variations of conductivity over a small control volume can be neglected. With the Favre filter, the additional subgrid term associated with the density-velocity correlation G U j /ρ has a large amplitude, of the same magnitude as the non-subgrid terms. In other words, the divergence of the Favre-filtered velocity is a poor approximation of the divergence of the velocity because of the small-scale variations of density and velocity. The density-velocity correlation is more significant in the energy conservation equation in the Favre formulation than in mass or momentum conservation equation in the Classical or Velocity formulations. The nonlinearity error ∂ j G Q j associated with the heat flux is significantly larger with the Favre filter than with the classical filter but remains rather small.

Ideal gas law

The ideal gas law (figure 10) is used to compute the filtered temperature from the filtered density. With the Favre filter, there is theoretically no subgrid term in

Equation Formulation

Large Medium Small or very small Mass Classical, Velocity

∂ j F ρU j , C j ρU j Favre C j ρU j Momentum Classical ∂ j F ρU j U i ∂ t F ρU i , C j ρU j U i C i P , ∂ j F Σ ij , C i Σ ij Favre ∂ j ρG U j U i C j ρU j U i C i P , ∂ j G Σ ij , C i Σ ij Velocity ρ∂ j F U j U i ρC j U j U i , ρU i C j U j ρF U i ∂ j U j , ρF ∂ i P/ρ , C i P , ρF ∂ j Σ ij /ρ , ∂ j F Σ ij , C i Σ ij Energy Classical, Velocity C j U j ∂ j F Q j , C j Q j Favre G U j /ρ C j U j ∂ j G Q j , C j Q j
Ideal gas law Classical, Velocity F 1/ρ Favre Table 1 -Classification of the subgrid terms in the three formulations. the ideal gas law. With the classical filter (Classical or Velocity formulation), the subgrid term F 1/ρ it is found negligible. The ideal gas law can thus be used without model with both the classical and Favre filter.

Summary

The subgrid terms are classified according to their quadratic average in table 1. The large subgrid terms are the most significant subgrid terms, and should be modelled first. The medium subgrid terms have a smaller amplitude. The modelling of the small or very small subgrid terms should only be considered after all large and medium subgrid terms are modelled, and may not be recommended as their amplitude may not be larger than even the modelling error of the larger subgrid terms.

Effect of the filter width

The effect of the filter width on the subgrid term is investigated at Re τ = 180 by comparing the amplitude of the subgrid terms with the filters A, B and C. As the filter size is increased, the magnitude of all subgrid term is increased compared to the non-subgrid terms. In the mass conservation equation, the momentum conservation equation and the ideal gas law, the increase is not sufficient to affect the conclusions of the classification of the subgrid terms, since the order of magnitude of the subgrid terms remains similar for the three filters. This supports the assumption that, while the variations across the channel of viscosity and thermal conductivity have a significant impact on turbulence [START_REF] Serra | Thermal large eddy simulation in a very simplified geometry of a solar receiver[END_REF], the effects of their small-scale variations on the flow are negligible.

In the energy conservation equation, the filter width has a strong influence on the magnitude of the filter-derivative non-commutation subgrid term C j U j associated with the divergence of the velocity. While the magnitude of this subgrid term is moderate with the filter A (figure 9), it is very large with the filter C (figure 11), as its amplitude is larger than that of the resolved heat flux from the logarithmic region to near the centre of the channel. The results with the filter B (not presented here) are in-between. The differences are primarily due to the increased anisotropy of the filter when the filter width is increased. The three filters follow a hyperbolic tangent law in the wall-normal coordinate direction, with a larger dilatation parameter the larger the filter width. With the filter C, the subgrid term C j U j is rather small in the viscous sublayer and becomes very large farther from the wall. This suggests that a filter more uniform in the logarithmic region is preferable to limit the influence of this subgrid term.

The conclusions of the classification of the subgrid terms (1) are valid for the filter A and B. For the filter C, C j U j is a large subgrid term. We consider its modelling necessary.

Conclusion

The low Mach number equations are suited to turbulent flows with a low Mach number but subjected to large variations of the fluid properties. They are characterised by a distinctive form of the energy conservation equation, that does not let us categorically choose between the classical filter and the Favre filter. In this study, we filter the low Mach number equations and identify the specific subgrid terms. Then, we investigate a priori the magnitude of all subgrid terms using the flow field from direct numerical simulations of a strongly anisothermal turbulent channel flow. The temperature gradient creates an asymmetry between the hot and cold sides regarding the amplitude and position of the maxima of the subgrid terms, explained by the local variations of the mean fluid properties and a low Reynolds number effect. Regardless of the formulation, more than half of the subgrid terms are found negligible. In particular, the effect of small-scale variations of viscosity and thermal conductivity may be neglected. The two most significant subgrid terms are the subgrid terms associated with momentum convection and with the density-velocity correlation. Due to the mesh inhomogeneity, the modelling of some filter-derivative non-commutation may also need to be considered, depending on the width of the selected filter. The classical filter is found more appropriate if the momentum equation is expressed as the velocity transport equation. The Favre filter removes the need for the modelling of the density-velocity correlation from the mass conservation equation but requires the modelling of an additional subgrid term in the energy equation, which has a very significant amplitude. The density-velocity correlation thus needs to be modelled in both cases.

Figure 1 -

 1 Figure 1 -Biperiodic anisothermal channel flow.

Figure 2 -

 2 Figure 2 -Surface of the 400 K isotherm and velocity magnitude in a plane normal to the flow direction (Re τ = 180).

Figure 3 -

 3 Figure 3 -Validation of the spatial turbulence kinetic energy terms.

+ x = 43

 43 , ∆ + y = 0.8 -21, ∆ + z = 21 and 77 × 53 × 77 grid points; • filter B: ∆ + x = 68, ∆ + y = 0.5 -25, ∆ + z = 34 and 48 × 50 × 48 grid points; • filter C: ∆ + x = 91, ∆ + y = 0.13 -43, ∆ + z = 45 and 36 × 40 × 36 grid points. At Re τ = 395, only filter A is examined. The filter sizes in wall units are ∆ + x = 47, ∆ + y = 1.5 -21 and ∆ + z = 24. This corresponds to a mesh with 154 × 102 × 102 grid points.

Figure 4 -

 4 Figure 4 -Root mean square of the terms of the mass conservation equation as a function of the wall-normal coordinate in the Classical and Velocity formulations (a) and in the Favre formulation (b) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  a) In the Classical formulation.

  b) In the Favre formulation.

  c) In the Velocity formulation.

Figure 5 -

 5 Figure 5 -Root mean square of the terms of the streamwise momentum conservation equation as a function of the wall-normal coordinate in the Classical formulation (a), in the Favre formulation (b) and of the terms of the streamwise velocity transport equation in the Velocity formulation (c) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  a) In the Classical formulation.

  b) In the Favre formulation.

  c) In the Velocity formulation.

Figure 6 -

 6 Figure 6 -Root mean square of the terms of the spanwise momentum conservation equation as a function of the wall-normal coordinate in the Classical formulation (a), in the Favre formulation (b) and of the terms of the spanwise velocity transport equation in the Velocity formulation (c) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  a) In the Classical formulation.

  b) In the Favre formulation.

  In the Velocity formulation.

Figure 7 -

 7 Figure 7 -Root mean square of the terms of the wall-normal momentum conservation equation as a function of the wall-normal coordinate in the Classical formulation (a), in the Favre formulation (b) and of the terms of the wall-normal velocity transport equation and in the Velocity formulation (c) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  In the Favre formulation.

  In the Velocity formulation.

Figure 8 -

 8 Figure 8 -Root mean square of the terms of the resolved kinetic energy transport equation as a function of the wall-normal coordinate in the Classical formulation (a), in the Favre formulation (b) and in the Velocity formulation (c) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  With the Favre filter.

Figure 9 -

 9 Figure 9 -Root mean square of the terms of the energy conservation equation as a function of the wall-normal coordinate in the Classical and Velocity formulations (a) and in the Favre formulation (b) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  With the Favre filter.

Figure 10 -

 10 Figure 10 -Root mean square of the terms of the ideal gas law as a function of the wall-normal coordinate in the Classical and Velocity formulations (a) and in the Favre formulation (b) at Re τ = 180 with the filter A. The amplitude is scaled by the maximum value in the domain in the three formulations.

  With the classical filter.

  With the Favre filter.

Figure 11 -

 11 Figure 11 -Root mean square of the terms of the energy conservation equation as a function of the wall-normal coordinate in the Classical and Velocity formulations (a) and in the Favre formulation (b) at Re τ = 180 with the filter C. The amplitude is scaled by the maximum value in the domain in the three formulations.

Acknowledgment

The authors gratefully acknowledge the CEA for the development of the TRUST platform. This work was granted access to the HPC resources of CINES under the allocations 2016-c20162a5099 and 2017-A0022A05099 made by GENCI.


Victor Munoz Cuartas

Marianne Perrin

Marie-Laetitia Pastor

Hélène Welemane

Arthur Cantarel

Moussa Karama

Victor Munoz

Determination of the elastic properties in CFRP composites: comparison of different approaches based on tensile tests and ultrasonic characterization

Keywords: carbon fibre composite, elastic constants, ultrasonic characterization, tensile test, stiffness tensor

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

Introduction

Fibre-reinforced composites are widely used for many structural applications. The primary benefits of the composite components are the reduction of weight and the simplification of assembly [START_REF] Soutis | Fibre reinforced composites in aircraft construction[END_REF]. The determination of the mechanical properties is essential for ensuring performance to the composite structures. The determination of their stiffness is critical to provide performance to the structure. In addition, the knowledge of complete elastic stiffness matrix is important for modeling and evaluating the mechanical behavior of composite materials under loading conditions [START_REF] El Bouazzaoui | Development of damage in a 2D woven C/SiC composite under mechanical loading: II. Ultrasonic characterization[END_REF].

For now, tensile test with strain gauge measurements is the technique normalized by ISO standards to identify composites elastic properties. Such conventional method is destructive in nature and provides only a part of elastic constants when thin plates like laminated composite structures are considered. Accordingly, ultrasonic techniques based on the measurement of Corresponding author, Ph.D., E-mail: marianne.perrin@iut-tarbes.fr ultrasonic wave velocities provide an interesting and non destructive way to address such issue.

First experimental studies based on ultrasonic technique were done by [START_REF] Zimmer | Determination of elastic constants of a unidirectional fiber composite using ultrasonic velocity measurements[END_REF]. They have used transmission contact technique to obtain dynamic elastic stiffness necessary to derive the unidirectional composite elastic behavior. Several ultrasonic bulk wave methods have also been developed to obtain phase velocities in anisotropic plates, especially when only one or two sides of the sample are accessible [START_REF] Chu | On determination of orthotropic material moduli from ultrasonic velocity data in nonsymmetry planes[END_REF], Vishnuvardhan et al. 2007).

Immersion ultrasonic technique represents a particularly suitable method in the case of composite materials with small thickness [START_REF] Munoz | Indetification of the elastic properties in CFRP composites[END_REF][START_REF] Reddy | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF]. For a transversely isotropic composite, measurements in symmetry planes are sufficient to entirely determine the stiffness tensor (that is the five independents constants) [START_REF] He | Characterization of nonlinear shear properties for composite materials using digital image correlation and finite element analysis[END_REF]. Various numerical methods and/or technical devices have then been implemented to interpret experimental results (Balasubramaniam andWhitney 1996, Kawashima et al. 1998).

Generally speaking, two experimental protocols can be used for such immersion technique to get velocities measurements, namely through-transmission and back-reflection techniques [START_REF] Reddy | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF]. Through-transmission method requires two transducers: one to send the wave through the sample and the second one to receive the transmitted wave. As the ultrasonic wave travels through the test sample, the wave is reflected in part as it encounters a medium of different acoustic impedance. Then the transmitted wave is received by the transducer and displayed or stored for analyses.

In the back-reflection technique, a transducer working in pulse/echo mode is associated to a large flat reflector which is positioned parallel to the transmitter. The back-reflected wave travels exactly along the same path as the incident wave in the opposite direction [START_REF] Rokhlin | Double through-transmission bulk wave method for ultrasonic phase velocity measurement and determination of elastic constants of composite materials[END_REF]. When the sample angle is changed, the position of the incident wave on the back reflector changes. Contrary to through-transmission, it is yet not necessary in this case to move the reflector or the transmitter/receiver transducer due to the large dimensions of the reflector.

The aim of this paper is to compare the back-reflected immersion ultrasonic method with classical mechanical characterizations based on tensile tests. Strain gauges are used to measure the axial and lateral strains and stiffness constants are then estimated from the elastic parts of the stress-strain response. At the same time, a digital image correlation system is also implemented to corroborate the gauge data. Regarding immersion technique, a specific device including a rotation system has been set up to study the response of the material under various incident waves. These three results of stiffness tensor measurement (strain gauges, digital image correlation during tensile tests and ultrasonic characterization) are finally compared and discussed.

Specimen and experimental procedure

Specimen

In view of its increasing use in aeronautics, a carbon fibre reinforced laminate is considered for this study. The M10R/38%/UD150/CHS composite is made of 14 unidirectional plies of prepreg leading to a thickness composite ℎ = 2 mm. Sample fabrication is carried out using the manual lay up technique. Then the sample is cured at 125°C during 90 minutes at a pressure of 2 bars. The dimensions of the unique sample used for the ultrasonic characterization are 150 x 100 x 2 mm 3 . The density of this composite is measured experimentally, namely 𝜌 = 1449 kg/m 3 . This material is representative of a transversely isotropic composite.

Experimental procedure for the tensile test method

The elastic behavior determination based on mechanical tests requires to apply axial tensile load with different orientation regarding the symmetry axis of the material. In this way, laminated samples with three axis fibre directions (0°, 90° and 45°), (see Fig. 1) have been cut according to the standard ISO 527-5 for unidirectional composites (Standard NF EN ISO 527-5 2009). With those tests, only 4 elastics constants can then be determined since loading in the direction of axis 3 (orthogonal to the laminate plane) cannot be done for thin composites plates. According to the standard, the displacement rate is 2 mm/min for 0° samples, 1,5 mm/min for 45° samples and 1 mm/min for 90° samples. To check the reproducibility of the response, three specimens are considered for each loading direction.

Strain gauges of 350 Ohm (HBM K-LY41-6/350-3-2M) are connected to a data acquisition system HBM Spider 8.30 and the data processing is made with Catman Easy software. Digital image correlation system is also performed during tensile tests. Such technique aims at matching two digital images of a surface observed at two different states of load, generally in the reference state (unloaded) and in a deformed state [START_REF] Sutton | Advances in two-dimensional and threedimensional computer vision[END_REF].

The markers used for matching the images with the DIC system are done by painting white points on the sample surface (see Fig. 2). The processed image with four markers allows then to obtain the strain during the load for each direction (axial and lateral). The DIC uses an Aramis sensor made of two CCD cameras 1392 x 1040 pixels. However for plane structures, only one camera is needed for measuring the 2D strain. These cameras were related to an Aramis software for digital processing. Gauges and monitored markers are positioned on the same face of the specimen.

𝐺 12 = 𝐸 45 -𝐸 1 -𝐸 2 + 𝜈 12 𝐸 1 (1)
As said before, this procedure does not allow to obtain the last elastic property, namely the Poisson ratio , that would entirely determine the elastic behavior of the laminated composite.

Experimental procedure for ultrasonic characterization

General principle

The elastic stiffness matrix for a transversely isotropic material is given by the following tensor written in the coordinate axes described in Fig. 1:

[𝐶 𝑖𝑗𝑘𝑙 ] = [ 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 𝐶 ] (2)
Ultrasonic characterization aims at providing five components representative of the elastic behavior of such material, namely 𝐶 , 𝐶 , 𝐶 , 𝐶 and 𝐶 (𝐶 can be derived from 𝐶 and 𝐶

). This determination relies on the resolution of the Christoffel equation [START_REF] Rose | Ultrasonic waves in solid media[END_REF]:

| 𝐶 -𝜌 ( ) | = (3)
where ( ) is the wave velocity related to the propagation direction of unit vector , ρ is the material density and is identity second-order tensor. Accordingly, this requires the measurement of ultrasonic velocities for different orientations of wave propagation in the sample. For a transversely isotropic composite, measurements in planes of symmetry are sufficient to determine all five independent elastic constants [START_REF] Aristégui | Optimal recovery of the elasticity tensor of general anisotropic materials from ultrasonic velocity data[END_REF]Baste 1997, Chu et al. 1994).

For the immersion ultrasonic technique, water acts as a couplant that transfers the wave from the transducer to the sample under inspection (Castagnède et al. 1990, Rokhlin and[START_REF] Rokhlin | Double through-transmission bulk wave method for ultrasonic phase velocity measurement and determination of elastic constants of composite materials[END_REF]). The transducer is not directly connected with the sample and hence consistent coupling is ensured. In this way, it is possible to measure the wave velocities at different angles of propagation either by adjusting transducer orientation or by rotating the sample [START_REF] Reddy | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF]. Here the sample is held and moved using a rotation system (see Fig. 3). As ultrasonic generator, an Omniscan 32: 128 PR is used with a mono-element transducer connected to it. The transducer acts at the same time as ultrasonic source and receiver (reflection mode). Classical value of frequency is used, namely 5 MHz, which allows mainly to avoid inside reflection; back-reflector and specimen were checked for alignment at normal incidence. The measurements are performed for different angles between plane 2-3 and incident plane (see Fig. 4). For a given , the sample is then rotated in the incident plane for different incident angles 𝜃 𝑖 between the transducer and axis 3 of the sample. The wave passes through the sample in the direction of with angle 𝜃 , so 𝜃 is function of and 𝜃 𝑖 . For each set of angles ( and 𝜃 𝑖 ), the time of flight is measured.

Procedure for the elastic constants estimation

The determination of the velocity of propagation in a medium requires the exact time of travel in the specimen for the given angle of propagation. Such time is determined using the crosscorrelation technique [START_REF] Rao | New methodology for precise ultrasonic velocity measurement and its applications[END_REF]. It requires a reference signal, which corresponds here to the signal recorded without the sample (only water). The time to cross the material is obtained by subtracting the time of flight along the reference path in the coupling medium (water) to the overall time of flight along the path in the sample at a given incident angle. This time difference gives then the ultrasonic velocity (𝜃 ) in the specimen (m/s) for the refracted angle 𝜃 by the following relation (Vishnuvardhan et al. 2007):

( ) = ( 𝜃 𝑖 ) = ( 2 - + ( ) 2 2 ) (4)
where is the velocity in the immersion fluid (water) and 𝜃 𝑖 the incidence angle (radians). The transducer senses both longitudinal and quasi-shear waves. Yet, the longitudinal velocity of the first one is always more important than the quasi-shear velocity of the latter; this allows to distinguish their own contribution in the velocity signal [START_REF] Reddy | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF]. These both waves are generated when the incident wave encounters the sample

In the present case, only two configurations of the incident plane are necessary to derive the elastic constants, namely = and = . Fig. 5 shows the ultrasonic velocities obtained in plane 2-3 (that is for = ). The longitudinal velocity and quasi-shear velocity are almost constant for all angles of propagation, which stands in agreement with the transverse isotropy of the material around axis 1. The slight deviation of velocity (only 8% for a range of 28° of incident angle) can be explained by some misalignment of plies during the manufacturing process. In contrast, regarding the ultrasonic velocities in plane 1-3 (that is for = , Fig. 6), both longitudinal and shear velocities are clearly affected by the incident angle. This confirms the anisotropic character of the material in this plane (Vishnuvardhan et al. 2007). Note finally that, for incident angles bigger than 60°, no signal can be detected. Ultrasonic measurements in plane 2-3 are directly related to the constants 𝐶 and 𝐶 . Component 𝐶 is determined from normal incidence longitudinal velocity, that is for 𝜃 𝑖 = (point A in Fig. 5):

𝐶 = 𝜌 ( = 𝜃 𝑖 = ) (5)
On the other hand, 𝐶 is determined from the average shear velocity 𝑚𝑒𝑎𝑛 ( = ) data measured in 2-3 plane (part B in Fig. 5):

𝐶 = 𝜌 ( 𝑚𝑒𝑎𝑛 ( = )) (6)
𝐶 , 𝐶 and 𝐶 are obtained from both longitudinal and shear velocity data measured in 1-3 plane. Experimental velocity data is fitted in the following two relations obtained from the Christoffel equation's solution (Eq. ( 3)). This allows then to determine the three unknown parameters.

( = 𝜃 𝑖 ) = √ 𝐴+√𝐴 2 𝐵 𝜌 (7)
( = 𝜃 𝑖 ) = √ 𝐴 √𝐴 2 𝐵 𝜌 (8)
where

𝐴 = (𝐶 cos 𝜃 𝑖 + 𝐶 sin 𝜃 𝑖 + 𝐶 ) (9) 𝐵 = 𝐶 𝐶 𝑐𝑜𝑠 𝜃 𝑖 + 𝐶 𝐶 𝑠𝑖 𝜃 𝑖 + 𝑖𝑛 2 [𝐶 𝐶 + 𝐶 -(𝐶 + 𝐶 ) ] (10)
To solve undetermined Eqs. ( 7)-( 8), one can use the nonlinear least-square optimization technique which minimizes the deviations between the experimental and theoretical velocities for the considered angles of propagation [START_REF] Reddy | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF]:

min 𝐶 𝑗𝑘𝑙 ∊ 𝑛 ∑ ( 𝑖 𝑒 -𝑖 ) 𝑚 𝑖= (11)
In Eq. 11, is the number of independent parameters to be extracted (here 5 elastic constants) and is the number of measurements of velocities in different directions (here 137 experimental data). 𝑒 and are the experimental and theoretical phase velocities, respectively. In a last step, component 𝐶 is calculated using this relation between elastic stiffness components:

𝐶 = 𝐶 -𝐶 (12)
Engineering moduli in the coordinate system ( , = , , and ) can be deduced by inverting the stiffness tensor.

Results

Table 1 and Fig. 7 show the elastic constants estimation obtained by means of the three different methods used to characterize the carbon fibre/epoxy composite. Since tensile test with strain gauge measurements is the technique normalized by ISO standards, results are compared to this method. Indeed, we can see a good accuracy between the 0° tensile test with strain gauges and the only data provided by the manufacturer, namely the axial Young modulus . Regarding DIC technique, we can note a very accurate determination of Young moduli (error of 5% on the mean) and

(error of 2% on the mean), whereas results are less precise for Poisson ratio (error of 15 % on the mean) and shear modulus (error of 8% on the mean). Latter moduli are determined from lateral strains measurement. For carbon fibre composite subjected to axial tension, these strains exhibit a low amplitude and are then more affected by measurement inaccuracies. However, elastic constants obtained by DIC clearly stand in the range of standard measurements if we consider the standard deviations of each technique. DIC appears then as an interesting alternative to strain gauge, especially as DIC avoids contact and is more practical to set up.

The value of in US characterization is the most accurate result in regards to the normalized tensile test (error of 11% on the mean) and corresponds to a good estimation if the standard deviation is taken into account. Results on mean values of (error of 22% on the mean) and Poisson ratio (error of 23 % on the mean) seem less accurate but, again, they stand in the same range as the reference according to the standard deviation of both techniques. The higher error is found for (error of 76 % on the mean). All these deviations can be due to the estimation of the time of flight of the wave inside the material. Basically, the time of flight is very sensitive to the time reference measured between the transducer and the back-reflector as explained in section 2.3.2. It should be finally underlined that US technique is the single one that provides Poisson ratio .

Conclusions

This paper intended to compare different methodologies to establish the elastic behavior of composite materials. Transversely isotropic composites were considered here through the study of unidirectional carbon fibre/epoxy laminates. Two strain measurement techniques (gauge as standard, Digital Image Correlation) during tensile tests and the immersion ultrasonic method were analyzed.

Results highlight first the interest of DIC to determine the strain response and its good accuracy on the final determination of elastic constants. Yet, tensile tests do not allow to entirely provide the elastic stiffness and some indetermination remains for structures calculation.

This study has shown also the ability of the immersion ultrasonic technique to derive all elastic constants of a UD laminated composite through the measurement of time delay of waves inside a sample. The main advantage of this technique is that only one specimen is needed instead of cutting many samples with desired size and shape for a classical characterization tensile test. If elastic properties estimations are quite encouraging, some limitations on the precision have been noted that would require further investigations on the influence of the time of flight. This could be improved by adjusting the experimental setup in such a way that the distance between transducer and reflector is fixed during all the acquisition or by using two transducers (through-transmission).

Given that, this immersion ultrasonic method could then be extended with the same principle to more complex anisotropic materials including orthotropic symmetry. This would need another degree of freedom to rotate the sample that will provide the possibility to propagate ultrasonic waves along more azimuthal angles/planes with respect to the fibre axis.

Fig. 1

 1 Fig. 1 Setup for tensile characterization method. Plane 1-2 is the laminate plane

Fig. 2

 2 Fig. 2 Setup for tensile characterization method

  (a) 𝜃 𝑖 = (b) 𝜃 𝑖 ≠ Fig. 3 Back-reflection setup

  path description. : angle between plane 2-3 and incident plan direction. 𝜃 𝑖 : incident angle. 𝜃 : refraction angle

Fig. 5

 5 Fig. 5 Ultrasonic velocity vs. incident angle in plane 2-3 ( = )

Table 1

 1 Calculated elastic constants from the three different methods

	Test		[GPa]	[GPa]			[GPa]
	Tensile test (Gauges)	Mean SD	138 10	8.2 1.3	0.34 0.03		5.1 0.3
	Tensile test (DIC)	Mean SD	130 4	8.4 0.8	0.29 0.04		4.7 0.2
	Ultrasonic (US) characterization	Mean SD	123 5	6.4 1.1	0.42 0.07	0.13 0.02	9.0 0.8
	Manufacturer (Hexcel)	Mean	140


Victor Munoz
email: victor.munozcuartas@enit.fr

Marianne Perrin
email: marianne.perrin@iut-tarbes.fr

Marie-Laetitia Pastor
email: marie-laetitia.pastor@iut-tarbes.fr

Hélène Welemane
email: helene.welemane@enit.fr

Arthur Cantarel
email: arthur.cantarel@iut-tarbes.fr

Moussa Karama
email: moussa.karama@enit.fr

APPLICATION OF NON DESTRUCTIVE TESTING TO THE DEETCTIONS OF AERONAUTICAL DEFECTS IN COMPOSITE STRUCTURES

Keywords: Non-destructive testing, Ultrasonics, Infrared thermography, Composites, Simulated defects

A study of two Non-destructive Testing methods (NDT) was carried out in specimens with different kinds of simulated defects. Ultrasonics test (US) and Infrared Thermography (IRT) were applied with the aim of evaluate the detectability and the accuracy of each method. These techniques acquired great importance in the aeronautics industry because they control the aerostructure without intervening in their physical and mechanical integrity. In the second part of the study, a comparison of both techniques was done with the purpose of analysing their limits and advantages. It appeared that detectability of defects was much better in a sample with flat-bottomed holes defects in the case of Ultrasonic Test. However it was found that Infrared Thermography is far more limited to the thickness of the specimen than the ultrasonic waves. On the other hand, defects were all revealed with IRT in a sandwich composite including Teflon inserts, which was not the case for US.

INTRODUCTION

The use of composites in the aerospace industry has increased dramatically since the 1970s. The primary benefits that composite components can offer are reduced weight and assembly simplification. It happens that different types of mechanical and thermal loads are applied to these structures when they are in service, generating internal stresses. It is in this way that delamination or disbonding may result in the aerostructure. Defects such as random inclusions or undesirable material have also been found during the manufacture process of composite materials [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF]. These flaws may create stress concentration with serious consequences. It is therefore important to test the composite structure to ensure their integrity. These defects have to be revealed in order to make the correct maintenance or replacement. The faster the damage or defect is detected, the safer the aerostructure is.

Many techniques have been used with the purpose of detect such defects and in the operational configuration. Liquid penetrant inspection, magnetic control inspection, eddy currents control, radiographic testing and shearography are some of the methods used for evaluating the health of the material, without altering its properties [START_REF] Cherfaoui | Essais non destructifs[END_REF]. These kinds of controls are called the Non-Destructive Testing (NDT). The most widespread technique is the ultrasonic testing. It is the only technique certified by the aeronautics industry [START_REF] Boro-Djordjevic | Ultrasonic characterization of advanced composite materials[END_REF]. Ultrasonic Test (US) uses high frequency waves to conduct examinations and make measurements. Ultrasonic inspection is used for flaw detection/evaluation, dimensional measurements and material characterization [START_REF] Budnik | Essais non destructifs -Contrôle par ultrasons[END_REF]. Some studies have already shown the detectability and the accuracy of a flaw in size, shape and deepness for different kinds of composite structures used in the aeronautic industry [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF][START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF].

A full-field measurement technique based on Infrared Thermography (IRT) is also used in the NDT field. Infrared thermography is a nondestructive, non-intrusive, non-contact mapping of thermal patterns or thermograms, on the surface of objects. IRT is more widely used in recent years for structural investigation [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF]. The mean difference about IRT with regard to US is that IRT does not give the information about the depth of defects. Some works aimed to compare different infrared thermography configurations, have detailed the accuracy of defects detection [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF][START_REF] Péronnet | Non destructive investigation of defects in composite structures by three infrared thermographic techniques[END_REF].

The aim of this work is to evaluate two different specimens by ultrasonics test and infrared thermography. A carbon-epoxy laminated composite and a carbon-epoxy-glass sandwich composite were tested to verify the accuracy and detectability of their simulated defects by using both NDT methods. Specimens have different kinds of defects and were made in different materials used in the aeronautic industry. Defects have various diameters and locations on the material so that the limit of size and detection in deep can be determined. The second objective is to compare the information obtained by each technique. Finally, a conclusion about the intention of capitalize the information obtained by NDT methods is presented.

SPECIMENS AND EXPERIMENTAL PROCEDURE 2.1 Specimens

The specimen used to compare the detectability and accuracy between the IRT and the US is a carbon/epoxy laminated composite with flat bottomed holes defects and variable thickness. It is called specimen A (figure 1). This sample is divided in 2 sections. The first one is 4 mm width in 50 mm length. The second section is 130 mm length with 8 mm width. This sample contains 15 holes with diameters ranging from 2 mm to 10 mm. These holes have different deepness (from 0.52 mm until 7.5 mm). The aim of these holes is to simulate discontinuity on the specimen that is considered as a continuous medium. A second sample, called specimen B was also used. It was tested only in infrared thermography because of its materials properties. This specimen is a carbon-glass-epoxy skin with foam-core sandwich composite including Teflon inserts. The thickness of its skin is 0.66 mm. It has been demonstrated that the foam-core is a sound insulation, whereby, the ultrasonic waves can not penetrate this specimen to be detected by the receiving transducer [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF]. This specimen contains 18 defects. These simulated flaws are Teflon inserts ranging from 3 mm to 6 mm. They are located at the skin section. These inserts were placed during the stacking of the tissues so that they simulate as best as possible a delamination in the sample. The manufacturing drawing of this specimen is showed in figure 2.

Experimental protocol

Ultrasonics

In ultrasonics, the sound energy is introduced and propagates through the material in the form of waves. When there is a discontinuity (such as inclusion, delamination or disbonding) in the wave path, part of the energy will be reflected back from the flaw surface. The reflected wave signal is transformed into an electrical signal by the transducer and is displayed on a screen. There are two configurations that can be used in an US testing: reflection and transmission. In this study, only reflection configuration is used in specimen A.

This technique allows to get information about the modulus of elasticity by measuring the wave velocity [START_REF] Baste | Evaluation de la matrice d'élasticité des composites orthotropes par propagation ultrasonore en dehors des plans principaux de symétrie[END_REF]. Deepness, size and shape of a detected defect, can also be obtained by this technique.

In this study, a Omniscan 32: 128 PR is used with a 5L64-NW1 multi-elements transducer. The Omniscan is both, ultrasonic source and receiver at the same time. The frequency of the transducer is 5 MHz, which is the common frequency in aeronautics NDT [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF]. For specimen A, the wave velocity is 2983 m s -1 . This experiment was carried out by contact; therefore it is necessary to use a gel coating between the sample and the transducer. Ultrasonics testing enables two-dimensional mapping, accordingly, a freescan was also used as a coding system to detect the position of the transducer in regard to the specimen. Specimen A was placed in a support to hold it while the experiment was being carrying out. Figure 3 shows the ultrasonics device.

Infrared Thermography

Infrared thermography (IRT) provides colourful images and movies of specimens where local changes in surface temperature indicate subsurface defects. Thermal waves are propagated through the material. When they pass by a different medium, a thermal gradient is generated in the specimen. Indeed, the two mediums have different emissivity coefficients. These differences are captured by an infrared camera enabling the emissivity coefficient to be converted to temperature. It is measured on the front of the specimen. Thermal two-dimensional mapping is created and inhomogeneities can be detected [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF][START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF][START_REF] Péronnet | Non destructive investigation of defects in composite structures by three infrared thermographic techniques[END_REF].

This technique allows to get information about the position, size and shape of the flaw. In this line of thought, delamination, disbonding or crack networks can be detected [START_REF] Garnier | The detection of aeronautical defects in situ on composite structures using Non Destructive Testing[END_REF]. Contrary from the US, information about the deepness can not be obtained.

To heat the sample, a 1000 W halogen lamp is used in two configurations: transmission and reflection. Specimen A was tested only in transmission mode (figure 4a) and specimen was tested in both configurations (figure 4a and4b). The infrared sensor was an IR camera whose thermal resolution was 20 mK. A 30 s and 50 Hz movie was recorded in which 15 s corresponds to the heating of the sample and the 15 last seconds to its cooling. The treatment of the movie was done with the ALTAIR software. The analysis of the movie is done in relative mode with the purpose of removing the ambient temperature. Lock-in IRT was also used on both specimens. The idea behind lock-in IRT is that temperature modulation induced at the surface of the inspected component from the outside propagates as a harmonic thermal wave. Lock-in IRT allows better defect inspection than common IRT and it is less sensitive to environmental conditions [START_REF] Choi | Quantitative determination of a subsurface defect of reference specimen by lock-in infrared thermography[END_REF]. The connection of the lock-in IRT is detailed on figure 5. A signal generator is also needed in addition to the classical IRT.

RESULTS

Ultrasonics

As it was said on section 2, only specimen A was tested with ultrasonics, because material properties of specimen B do not let waves to propagate through the foam-core.

Analyses of results are made taking into account two parameters: detectability and accuracy. 14 of the 15 holes were detected by this method. This is the 93% of them. The only defect that was not detected, is the one located at 1 mm deep and 2 mm of diameter. However, it is noticed that defects located near the contact surface (less than 1 mm) are difficultly distinguish from the signal input peak. On the other hand, the detected defects were measured on average, 1 mm larger or smaller than they really are. Defects located at 7.5 mm deep (3 and 6 mm diameter) were measured 1 mm smaller than their real size, and those located at 5,2 mm deep (2, 3 and 6 mm diameter) were measured at the correct dimension. 10 mm diameter defects were all measured 1 mm larger than their real size. Table 1 shows the measured dimension for each range of defects.

Table 1. Defect measurement on specimen A using US Figure 6 shows the C-Scan mapping for this specimen in which measurements were done. Other works [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF] have obtained similar detectability with better accuracy for the same specimen. The gap between the real size and the measured size is 0.5 mm. Those measurements were done by using the -6 dB method. It consists of measuring the defect where the backwall echo, decreases 50%.

Infrared thermography

Specimen A

As in US testing, two parameters were taken into account to analyse the results: detectability and accuracy. A single picture must be selected from the recorded film. That selection depends on the quality of the image for each defect. In this study, for specimen A, 2 different images (figure 7) were selected to measure 2 ranges of defects located in two different sections of the specimen. These ranges of defects are located at 0.52 mm and 2.8 mm deep respectively and they are in the 4 mm section and 8 mm section of the sample. Figure 7 shows that only 6 from the 15 defects, were detected using classical IRT in transmission mode. It is noticed that heat saturation is generated after a short period of time. The 4 mm width section is rapidly saturated and in the 8 mm width section, no other defect is revealed with the increasing temperature.

The size of the measures defects has the same gap than in US testing. They were measured 1 mm larger than they really are. Table 2 shows the detailed size of the defects.

Table 2. Defect measurement on specimen A using IRT 6 mm diameter defect at 2.8 mm deep was detected but no measurements could be done because the temperature difference between defective and healthy parts was not enough to be measured. It is also remarkable that the number of detected defect decreases with deepness. While the defect is deeper, more difficult is the detection. In lock-in IRT, the number of detected defects in specimen A increased. However, no measurement of their size was done because of the temperature difference. In addition to classical IRT, 3 mm diameter defect at 2.8 mm deep and 5.2 mm deep defects (6 and 10 mm diameter) were visible as shown in figure 8.

Specimen B

In the case of specimen B, lock-in IRT has similar results as in classical IRT. All defects were detected and they could be measured in diameter.

Best results were obtained at 0.1 Hz for both specimens tested in lock-in IRT. It was founded that the textured surface of specimen B was confounded with the defects in IRT mapping. Measurements were therefore less accurate than classical IRT in specimen A.

Defects located at the camera side, were easily distinguished; on the other hand, defects located on the opposite side of the camera, were difficultly detected. In addition, the dimension of the measurements was 2 or 3 mm larger that the real size of the defect. Table 3 shows the founded size for each range of defects. Figure 9 shows the result on specimen B tested by lock-in thermography. Nowadays there are not specific methods for measuring the real size of a defect from an IRT mapping. The criterion depends on who does the measure. In some works, the size of defects has been calculated using a system of grey image measurement [START_REF] Péronnet | Identification expérimentale du comportement d'un fuselage composite : detection de défauts pae mesures de champs[END_REF].

ANALYSIS

A comparative analysis was done aimed to explain the reach of each technique. Table 4 summarizes the reach for each technique regarding different parameters. The C-Scan mapping in US shows that detectability is more affected by the diameter of defects rather than thickness. This detection becomes more difficult when the defect is close to the detecting surface (less than 1 mm). In this case the defect echo is confused with the echo of the detection surface.

By the other hand, it is clearly noticed that IRT is very limited to large thicknesses. Defects are not visible beyond 2.8 mm in a carbon-epoxy laminated composite. However, by the use of lock-in IRT, the reach on thickness became to 5.2 mm. At any rate, diameter of defects also affects the detectability in IRT. Figure 7b and 8 shows that flaws located in the same range of deepness with different diameter are not all detected. The heat diffusivity inside the material does not let to see temperature difference in the surface of the sample.

In the case of lock-in thermography, low frequencies are more convenient for detecting the flaws in both specimens. The use of a sinusoidal signal with low frequency (0.1 Hz) allows to optimize the wave spread in the specimen thickness. Other works [START_REF] Péronnet | Non destructive investigation of defects in composite structures by three infrared thermographic techniques[END_REF] with the same sample found that the best frequency for detectability of defects is 0.05 Hz.

For the carbon-glass-epoxy skin with foam-core sandwich composite, the use of infrared thermography is quite difficult since the foam is a thermal insulator. It is recommended to control the 2 skins of the sandwich separately and in reflection mode (figure 4b). Otherwise, the heat flux may not arrive to the detection surface.

CONCLUSION

The subject of this study was to compare the ultrasonics test and the infrared thermography as NDT techniques. By using 2 different samples it was shown that both methods have the same accuracy regarding the diameter of the flaws, which is not case for detectability parameter. 93% of defects in a carbon/epoxy laminated composite were detected by the US while only 60 % of them were detected by the IRT.

In the case of a carbon-glass-epoxy skin with foamcore sandwich composite, all defects were detected by IRT. Results were similar using both, classical IRT and lock-in IRT for this specimen, but the measured size was less accurate than the laminated composite. The textured surface was confused with the defects.

In future works, mechanical characterisation on carbon-epoxy laminated composites will be study by calculating the elastic constants from the velocity of the ultrasonic wave. Then NDT techniques will be used for monitoring the evolution of damage. Static and dynamic test (tensile, fatigue) will be carried out and at the same time, NDT techniques (Ultrasonics, Infrared Thermography and Acoustic Emissions) will monitor the evolution of the damage.

Figure 1 .

 1 Figure 1. Manufacturing drawing of the specimen A

Figure 2 .

 2 Figure 2. Manufacturing drawing of the specimen B

Figure 3 .

 3 Figure 3. Ultrasonic measurement

Figure 4 .

 4 Figure 4. Infrared measurement

Figure 5 .

 5 Figure 5. Lock-in thermography connection

Figure 6 .

 6 Figure 6. Lock-in thermography connection

Figure 7 .

 7 Figure 7. IRT mapping of specimen A Defects at 0.52 mm deep were measured using image 7a and those at 2.8 mm deep using image 7b.It is noticed that heat saturation is generated after a short period of time. The 4 mm width section is rapidly saturated and in the 8 mm width section, no other defect is revealed with the increasing temperature.

Figure 8 .

 8 Figure 8. Lock-in IRT in specimen A

Figure 9 .

 9 Figure 9. Lock-in IRT in specimen B Péronnet E. [5] tested the same sample and found a gap between real size and measured size ranging from 1.5 mm to 3.5 mm. This information is in accordance regarding the experiments on the present study.

Table 3 .

 3 Defect measurement on specimen B using lock-in IRT

	Real diameter (mm) Measured diameter (mm)
	6	8
	3	5

Table 4 .

 4 Characteristics of US and IRT

	Parameter	US	IRT
	Accuracy in size	1 mm (larger)	1 to 2 mm (larger)
	Minimal diameter	2 mm from 1 mm deep	2 mm from 0 mm deep
	Detection in deepness	More than 8 mm in laminated	2,8 mm for a 6 mm diameter defect in
		composites	laminated composites
	Detectability in the laminated	14/15. Non detected a 2 mm	6/15. Non detected defects beyond 5,2
	composite	defect at 1 mm deep	mm deep and 3 mm diameter
	Detectability in foam-core sandwich	NaN	18/18 by analysing both sides of the
	composites		sample
	Time for getting results	About 10 min	Instantly
	General limitation	Detection of defects close to the	Important specimen thickness
		inspection surface