A Le Bot

High frequency vibroacoustics: A radiative transfer equation and radiosity based approach

Keywords: High-frequency, geometrical acoustics, statistical energy analysis, radiosity, specular reflection, multiple diffraction

This paper is a review of the theoretical framework for application of radiative transfer equations to structural dynamics and acoustics. It is shown that under the assumption of geometrical acoustics and when the phase of rays is neglected, a representation of sound field in terms of incoherent fictitious sources provides a sufficiently large framework to embody diffuse as well as specular reflection, steady-state or transient phenomena and even diffraction. Two special cases are mentioned. The so-called radiosity method in acoustics corresponds to a purely diffusing boundary and statistical energy analysis when the sound field is diffuse.

Introduction

The interest of the scientific community in high frequency methods applied to vibro-acoustics is relatively recent compared to other fields of physics such as electrodynamics or acoustics. It seems that ray methods have not been so popular in structural dynamics where most engineers prefer to use softwares based on the governing equations solved by the finite element method.

However, this 'brute force' approach of vibro-acoustics rapidly encounters a practical limitation. When the frequency increases, the number of natural modes increases. Even for a numerical model of order of tens of millions of degrees of freedom, a reasonable task on current computers, the maximum analysis frequency is still too low and does not permit to cover the entire audio frequency band. This limitation is particularly unfavorable in the automotive industry or for larger structures such as aircraft or launchers which may contain millions of natural modes within the audio frequency band [START_REF] Shorter | Modeling noise and vibration transmission in complex systems[END_REF]. Furthermore, it might be emphasized that the sensitivity of these computations is poor which raises the question of how useful it is to know the 'exact' sequence of natural modes for systems having a large number of modes.

For these reasons, statistical energy analysis [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF] (SEA) has been developed for fifty years. Introduced early in the sixties by R.H. Lyon [START_REF] Lyon | Power flow between linearly coupled oscillators[END_REF] and many other authors [START_REF] Smith | Response and radiation of structural modes excited by sound[END_REF][START_REF] Newland | Calculation of power flow between coupled oscillators[END_REF][START_REF] Mace | The statistical energy analysis of coupled sets of oscillators[END_REF], this theory is based on statistical concepts and applies when sources are random, uncorrelated and when the number of modes is large enough to be considered as a statistical population. Then under certain circumstances [START_REF] Lafont | Review of statistical energy analysis hypothesis in vibroacoustics[END_REF], sound/vibration reaches a state of diffuse field and equipartition of modal energy (an interesting method to dispense with equipartition is proposed in Refs. [START_REF] Maxit | Extension of SEA model to subsystems with non-uniform modal energy distribution[END_REF][START_REF] Totaro | Extension of the statistical modal energy distribution analysis for estimating energy density in coupled subsystems[END_REF]). Then, a simple law emerges: the energy exchange between two coupled subsystems is proportional to the difference in their modal energies. This has served as a foundation to interpret the modal energy as a 'vibrational temperature' and the coupling power proportionality as a thermal law. SEA follows on from room acoustics and Sabine's formula and as such could be considered as a theory of statistical vibroacoustics [START_REF] Bot | Statistical vibroacoustics and entropy concept[END_REF] in the same way as exist statistical mechanics, statistical electromagnetism and so on.

Among numerous attempts to relax the more restrictive assumptions of SEA, several methods aim to predict the repartition of vibrational energy within subsystems where the field is not diffuse. But the behavior of local energy quantities is more subtle than it appears at first sight [START_REF] Carcaterra | Energy density equations and power flow in structures[END_REF][START_REF] Li | Energy fields of partially coherent sources[END_REF][START_REF] Bot | Geometric diffusion of vibrational energy and comparison with the vibrational conductivity approach[END_REF] and it is generally not possible to obtain a closed set of equations on energy variables without introducing strong assumptions.

To circumvent this difficulty, it has b een prop osed by several authors [START_REF] Bot | A vibroacoustic model for high frequency analysis[END_REF][START_REF] Franzoni | An acoustic boundary element method based on energy and intensity variables for prediction of high frequency broadband sound fields[END_REF][START_REF] Cotoni | Radiation of plane structures at high frequency using an energy method[END_REF][START_REF] Bot | Energy transfer for high frequencies in built-up structures[END_REF][START_REF] Kim | Prediction of sound level at high-frequency bands by means of a simplified boundary element method[END_REF][START_REF] Sarradj | High frequency boundary integral method as an alternative to statistical energy analysis[END_REF][START_REF] Guyader | Integral equation for frequency averaged quadratic pressure[END_REF][START_REF] Herrin | The energy source simulation method[END_REF][START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF][START_REF] Bot | Comparison of vibrational conductivity and radiative energy transfer methods[END_REF][START_REF] Sadoulet-Reboul | A hybrid method for vibroacoustics based on the radiative energy transfer method[END_REF][START_REF] Pratellesi | A hybrid formulation for mid-frequency analysis of assembled structures[END_REF] to start from an integral representation of the vibrational field and to calculate a quadratic quantity by neglecting to various degrees the coherence between the boundary sources. When all sources are incoherent, it results in a boundary integral equation on energy that may be solved by a collocation or Galerkin method. Although they have been developed in various contexts, all these methods suggest more or less the idea that vibrational energy propagates in structures like radiant thermal energy. This thermal analogy is clearly not equivalent to the thermal conductivity proposed in other respects (see Refs. [START_REF] Langley | On the vibrational conductivity approach to high frequency dynamics for two-dimensional structural components[END_REF][START_REF] Sestieri | Space average and wave interference in vibration conductivity[END_REF][START_REF] Bot | Comparison of vibrational conductivity and radiative energy transfer methods[END_REF] for a discussion on this point). The radiative transfer analogy is more apparent in room acoustics with the so-called radiosity method [START_REF] Kuttruff | Simulierte Nachhallkurven in Rechteckräumen mit diffusem Schallfeld[END_REF][START_REF] Kuttruff | Room acoustics[END_REF][START_REF] Kuttruff | Energetic sound propagation in rooms[END_REF][START_REF] Nosal | Investigation of the validity of radiosity for sound-field prediction in cubic rooms[END_REF][START_REF] Zhang | Relaxation of sound fields in rooms of diffusely reflecting boundaries and its application in acoustical radiosity simulation[END_REF]. Furthermore, when the vibrational field tends to be diffuse, it may be checked that the coupling power proportionality of statistical energy analysis is obtained as a limit case [START_REF] Bot | Derivation of statistical energy analysis from radiative exchanges[END_REF]. In this concern, let us also mention the recently proposed method called dynamical energy analysis which investigates the gap between full ray-tracing and statistical energy analysis [START_REF] Tanner | Dynamical energy analysisDetermining wave energy distributions in vibro-acoustical structures in the high-frequency regime[END_REF][START_REF] Chappell | Dynamical energy analysis for builtup acoustic systems at high frequencies[END_REF][START_REF] Chappell | Boundary element dynamical energy analysis: A versatile method for solving two or three dimensional wave problems in the high frequency limit[END_REF].

The purpose of this paper is to review the theoretical framework of an extension of the radiosity method to specular reflection and diffraction which may serve as a basis for a generalization of statistical energy analysis. The outline of the paper is as follows. In Section 2 we introduce some generalities about representation of a vibration field in terms of radiant energy supplied by actual and fictitious sources. In Section 3 are derived the balance equations governing the emissive power of fictitious sources. The special case of diffuse reflection and the radiosity equation is introduced in Section 4. Two problems which enlarge the scope of classical radiosity are then discussed, specular reflection in Section 5 and diffraction in Section 6. The paper finishes with some concluding remarks.

Energy representation of a wave field

The laws of radiative energy transfer is based on two principles; the linear superposition principle and Huygens' principle.

If the vibrational field is sufficiently 'disordered' then we may neglect all correlations between rays and energy variables turn out to be well-suited to describe the field state. Two rays emanating from a single source with a slightly different direction take different paths and are subjected to different events such as reflection, diffraction, diffusion etc. before eventually reaches the same receiver point. They rapidly 'forget' their phase so that after travelling a long distance, their energies may be added signifying that interference effects can be neglected.

Huygens' principle states that each point of reflection, refraction, diffraction etc. behaves like a secondary source from which emanate spherical waves. A more modern form is the so-called Helmoltz-Kirchhoff formula or any representation formula by Green's function. The principle of all these formulas is that in the presence of an obstacle, the field is unchanged if one removes the obstacle and replaces it by appropriate fictitious sources. The field is therefore the sum of a direct field (the field that would be without the obstacle) and the field produced by fictitious sources.

Let Ω be a domain whose regular boundary is Γ and whose set of diffracting points is ∆. We admit that Ω may contain volume sources ⇢.B yv i r t u e of Huygens' principle, we also put surface sources σ on Γ and diffraction sources λ i on ∆ (Fig. 1). The energy density W at any point x and time t is given by (assuming incoherence),

W (x,t)= Z Ω ⇢ (y,t 0 ) 4⇡cR 2 dΩ y + Z Γ σ (y, v,t 0 ) 4⇡cR 2 dΓ y + X i2∆ λ i (v,t 0 ) 4⇡cR 2 (1)
where R = |x -y| is the source-receiver distance, v =(xy)/R the sourcereceiver direction, t 0 = t -R/c the retarded time and c the sound speed (Fig. 1). The domain sources ⇢ have been assumed omnidirectional but σ and λ i depend on the emission direction v.

Inside the domain, we may also introduce the specific radiative intensity I(x, v,t)definedastheenergyflo wperunitsolidangleandunitareanormal to the rays. While the energy density W only depends on the position x and the time t, the specific intensity depends in addition on the direction vector v.C o n s i d e r i n gap o s i t i o nx and an infinitesimal solid angle d! about ad i r e c t i o nv,t h es pe c i fi ci n t e n s i t yi s ,

I(x, v,t)d! = Z d! ⇢ (y,t 0 ) 4⇡R 2 dΩ y + σ (y, v,t 0 ) 4⇡R 2 dΓ+ X i λ i (v,t 0 ) 4⇡R 2 (2)
where the integration is performed over the cone of vertex x and angle d!, dΓ is the part of boundary seen from x into the solid angle d! and the last sum is restricted to diffraction points lying in the cone (Fig. 1). An equivalent representation is obtained for two-dimensional domains by substituting the kernel 1/4⇡R 2 with 1/2⇡R.T h i si su s e f u lf o ri n s t a n c ef o r wave propagation in structural components such as plates, membranes, shells and so on. But when the medium is dispersive, c must be understood as the group velocity since this is the speed of energy flow.

The symbols ⇢, σ, λ may readily be interpreted in terms of energy flux. Consider a point source located at x 0 or equivalently a source density ⇢(y,t)= δ(yx 0 )⇢(t). The intensity vector is radial and has magnitude ⇢/4⇡R 2 at ad i s t a n c eR. By calculating the flux of intensity over an infinitely small sphere of radius ✏,i ty i e l d s , lim

✏!0 Z 4⇡ ⇢(t -✏/c) 4⇡✏ 2 ✏ 2 d! = ⇢(t)( 3 )
So, ⇢ is the power density of domain sources and therefore has units W/m 3 . We may assume that the p ower ⇢ of domain sources is always known. These sources represent actions of external loads and therefore must be specified for a well-posed problem. External loadings are usually expressed in terms of force or moment applied to a structure or volume velocity in a fluid. They may be specified either by their root-mean-square value or power spectral density. In all cases, asymptotic relationships based on point-mobility of the infinite system are available in the literature to convert these values into a power [START_REF] Skudrzyk | Vibrations of a system with a finite or an infinite number of resonances[END_REF][START_REF] Skudrzyk | Sound radiation of a system with a finite or an infinite number of resonances[END_REF][START_REF] Goyder | Vibrational power flow machines into built-up structures,part I: Introduction and approximate analyses of of beam and plate-like foundations[END_REF]. The specific intensity I(x, u,t)l e a v i n gt h es u r f a c ei nd i r e c t i o nu is calculated by dividing the infinitesimal power into the solid angle d! by the apparent normal surface cos ✓dΓw h e r e✓ is the polar angle with the normal to the boundary (Fig. 2). So for an infinitesimal surface dΓa n da ne m i s s i o n angle ✓,

I(x, u,t)d! cos ✓dΓ=lim ✏!0 σ(x, u,t-✏/c) 4⇡✏ 2 ✏ 2 d!dΓ( 4 )
The specific intensity leaving in direction u is therefore σ(x, u,t)/4⇡ cos ✓.A perfectly diffuse source has a constant specific intensity. In two-dimensional systems, the factor 4⇡ is replaced with 2⇡. Alternatively σ/4⇡ is the power per unit area of surface and solid angle and is sometimes called the directional emitted flux or directional emissive power [START_REF] Modest | Radiative heat transfer[END_REF]. It has units of W/m 2 . The radiosity B is the total energy flow leaving per unit area of the surface. However some authors prefer to define the radiosity as the incident energy rate. The radiosity is obtained by integrating the specific intensity over all emission directions, B = R

I cos ✓d!,g i v i n g B(x,t)= Z 2⇡ σ(x, u,t) d! 4⇡ (5)
where the integration is performed over the half-sphere of inward directions.

Radiative transfer equation

Ar e p r e s e n t a t i o no fa ne n e r g yfi e l dl i k ei nE q .( 1 )i m m e d i a t e l yr a i s e st h e question of the determination of fictitious source powers σ and λ.

Let the bidirectional reflectivity [START_REF] Modest | Radiative heat transfer[END_REF] R(v, u) of the boundary be the ratio of radiative intensity I(x, u,t)r e fl e c t e di nd i r e c t i o nu with emission angle ✓ and the incident flux I(x, v,t)cos' from direction v with incidence angle '.T h eb i d i r e c t i o n a lr e fl e c t i v i t yi sa l s or e f e r r e dt oa sreflection matrix [START_REF] Joyce | Exact effect of surface roughness on the reverberation time of a uniformly absorbing spherical enclosure[END_REF]. At a fixed point x of the boundary, the specific intensity I(x, u,t)l e a v i n gi n direction u is the sum of contributions from all incident directions v times R(v, u). It yields,

I(x, u,t)= Z R(v, u)I(x, v,t)cos'd! (6)
This equation expresses conversion in phase space of energy flux at any point of the boundary. But at a boundary point x,t h es p e c i fi cr a d i a t i v ei n t e n s i t y is I = σ/4⇡ cos ✓ while the incident flux is given in Eq. ( 2). Substitution gives,

σ(x, u,t) cos ✓ = Z Ω R(v, u)⇢ (y,t 0 ) cos ' R 2 dΩ y + Z Γ R(v, u)σ (y, v,t 0 ) cos ' R 2 dΓ y + X i R(v, u)λ i (v,t 0 ) cos ' R 2 (7
) This is an integral equation which specifies the unknown σ at any position x and in any direction u as dependent on σ at all other positions y and directions v.

As i m i l a rm e t h o da p p l i e sf o rλ.L e td e fi n et h ebidirectional diffractivity D(v, u) by analogy with the bidirectional reflectivity as the ratio of the emitted power per unit solid angle about u and the incident radiative intensity from direction v.T h et o t a lp o w e re m i t t e di st h es u mo fc o n t r i b u t i o n sf r o m all incident directions,

λ i (u,t) 4⇡ = Z D(v, u)I(x, v,t)d! (8)
By substituting Eq. ( 2),

λ i (u,t)= Z Ω D(v, u)⇢ (y,t 0 ) dΩ y R 2 + Z Γ D(v, u)σ (y, v,t 0 ) dΓ y R 2 + X j6 =i D(v, u)λ j (v,t 0 ) 1 R 2 (9)
Again, the unknowns λ i are determined by an integral equation where both λ j and σ are involved.

Radiosity method

The radiosity method has been introduced in acoustics by Kuttruff in Ref. [START_REF] Kuttruff | Simulierte Nachhallkurven in Rechteckräumen mit diffusem Schallfeld[END_REF][START_REF] Kuttruff | Room acoustics[END_REF]. Originally developed in thermics to compute the temperature of facing walls in thermal equilibrium by the so-called standard method or view factor method, this is a popular method in optics for illumination problems. The radiosity method may be readily transposed in acoustics provided that a finite sound speed is introduced. The advantage is that it provides a fruitful piece of information for transient regimes.

The radiosity method applies to rooms with perfect diffusely reflecting walls. Diffraction is usually not taken into account so that we set λ i =0i n Eqs. [START_REF] Shorter | Modeling noise and vibration transmission in complex systems[END_REF][START_REF] Lafont | Review of statistical energy analysis hypothesis in vibroacoustics[END_REF]. The reflectivity of a perfectly diffusing surface is the constant [START_REF] Joyce | Exact effect of surface roughness on the reverberation time of a uniformly absorbing spherical enclosure[END_REF]

R(v, u)= 1 -↵ ⇡ ( 10
)
where ↵ is the absorptivity of the wall. Most often, we may consider that the absorptivity does not depend on the incidence angle. Furthermore, a diffuse reflection follows Lambert's law which states that the specific intensity of a diffuse source does not depend on the emission direction. The brightness of ad i ff u s es u r f a c ei st h es a m ew h a t e v e rt h eo b s e r v a t i o na n g l e . T h eL a m b e r t law or cosine law is, σ(x, u,t)

4⇡ cos ✓ = K (11)
By a direct integration over all emission directions and a comparison with Eq. ( 14), it immediately follows that the constant K equals B/⇡.T h u s ,b y substituting Lambert's law [START_REF] Carcaterra | Energy density equations and power flow in structures[END_REF] and the diffuse reflectivity (10) into Eq. ( 7), we obtain the radiosity equation,

B(x,t)=(1-↵(x))  E(x,t)+ Z Γ B ✓ y,t- R c ◆ cos ' cos ✓ ⇡R 2 dΓ y (12)
where E(x,t)i st h ec o n t r i b u t i o no ft h ed i r e c tfi e l d ,

E(x,t)= Z Ω ⇢ ✓ y,t- R c ◆ cos ' 4⇡R 2 dΩ y (13)
In the two above equations, ' is the incidence angle at y, ✓ the emission angle at y and R = |x -y| the source receiver distance.

The most important consequence of the radiosity equation is that it allows to determine the sound decay after sources have been switched off even when the sound field is not diffuse. The radiosity equation turns out to be therefore an important generalization of Sabine's theory to predict the reverberation time in rooms of arbitrary shape and absorption factors. Miles [START_REF] Miles | Sound field in a rectangular enclosure with diffusely reflecting boundaries[END_REF] observed that after a rapid transient state, the sound ultimately decays exponentially at the same rate at each location in the enclosure. By introducing an exponential time dependence B(y,t)=B(y)exp(-λt)i nE q .( 1 2 )w i t h E =0,ityields, 14) by two algorithms (Gilbert [45], Kuttruff [START_REF] Kuttruff | Energetic sound propagation in rooms[END_REF])

B(x)=(1-↵(x)) Z Γ B (y)exp ✓ λ R c ◆ cos ' cos ✓ ⇡R 2 dΓ y (
The calculation of reverberation time then reduces to the solving of the above equation in λ.S e v e r a la l g o r i t h m sa n dc l o s e d -f o r ms o l u t i o n sh a v eb e e np r oposed in the literature [START_REF] Gerlach | Der nachhallvorgang als markoffsche kette theorie und erste experimentelle überprüfung[END_REF][START_REF] Carroll | Decay of reverberent sound in a spherical enclosure[END_REF][START_REF] Gilbert | An interative calculation of auditorium reverberation[END_REF][START_REF] Kuttruff | A simple iteration scheme for the computation of decay constants in enclosures with diffusely reflecting boundaries[END_REF][START_REF] Kuttruff | Energetic sound propagation in rooms[END_REF][START_REF] Nosal | Investigation of the validity of radiosity for sound-field prediction in cubic rooms[END_REF]]. An example is proposed in Table 1. The room is rectangular 20×30×5m 3 and has volume V =3 0 0 0m 3 ,s u r f a c ea r e aS =1 7 0 0m 2 and a uniform absorption ↵ =0.01 excepted in one side where a very high absorption ↵ =0.8 is prescribed. The reverberation time T r is estimated by the three standard formula Sabine, Eyring and Millington and by Eq. ( 14) with T r =1 3 .8/λ. In the uniform case (↵ =0 .01 on all surfaces) the sound field is diffuse. All methods give more or less the same result. When the floor 20 × 30 is highly absorbing, a strong difference appears between Sabine, Eyring and Millington formula since the difference between ↵ and -ln(1 -↵)i sl a r g ef o r↵ =0 .8. But, the radiosity equation gives a result of 0.8sw h i c hi ss i g n i fi c a n t l yl e s s than Sabine's result and close to Eyring's result. Finally, for a highly absorbing wall (10 × 5) the trend is reversed and all standard formula predict ar e v e r be r a t i o nt i m el o w e rt h a nt h eo n eo fr a d i o s i t y .

Another important application of the radiosity method is that it provides an alternative method to ray-tracing for the computation of sound pressure level maps in rooms with stationary sound sources. The formal proof of equivalence with ray-tracing and some numerical examples are given in Ref. [START_REF] Bot | Comparison of an integral equation on energy and the ray-tracing technique for room acoustics[END_REF][START_REF] Muehleisen | Steady state acoustic radiosity for the prediction of sound pressure levels in enclosures with diffuse surfaces[END_REF].

Specular reflection

The most intriguing property of the representation (1) is that it accounts for specular reflection [START_REF] Bot | A functional equation for the specular reflection of rays[END_REF]. To illustrate it, let us consider the simple situation of a steady-state point source in front of a perfect infinite mirror (Fig. 3). In Figure 3: Point source facing a perfect mirror. The source is at (0, 1), the receiver at (x, y) and the reflecting point at (ν, 0). The incidence angle is ϕ and the emission angle θ. at w o -d i m e n s i o n a ls p a c e ,t h ed i r e c tfi e l di s ,

W d (x, y)= 1 2⇡c p x 2 +(y -1) 2 (15)
if the source is at x =0andy =1andhasaunitpo w er. Iffurthermorethe reflecting line is y =0t h e nw es e tafi c t i t i o u sl a y e rσ(⌫, ✓)w h e r e⌫ is the abscissa and ✓ the emission angle.

In the right hand-side of Eq. ( 7), the source power ⇢ = δ(x)δ(y -1) is a delta Dirac function and the first integral reduces to R(', ✓)cos'/R.B u t for a perfect mirror, the specular bidirectional reflectivity is [START_REF] Joyce | Exact effect of surface roughness on the reverberation time of a uniformly absorbing spherical enclosure[END_REF],

R(', ✓)= δ(✓ -') cos ' (16)
Furthermore, the distance from the point source to the boundary point x = ⌫, y =0isR = √ ⌫ 2 +1 and the incidence angle is ' =arctan⌫.S u b s t i t u t i n g in Eq. ( 7) gives,

σ(⌫, ✓)=δ(✓ -arctan ⌫) cos ✓ √ ⌫ 2 +1 (17)
The reflected field is given by the second integral of Eq. ( 1). At a boundary point ⌫ the emission angle in view of the receiver point x, y is ✓ =a r c t a n ( x -⌫)/y and satisfies cos ✓ = y/ p (x -⌫) 2 + y 2 .T h er e fl e c t e d field is therefore,

W r (x, y)= Z 1 -1 δ(arctan x -⌫ y -arctan ⌫) cos ✓ √ ⌫ 2 +1 d⌫ 2⇡c p (x -⌫) 2 + y 2 (18)
After simplification,

W r (x, y)= 1 2⇡c p x 2 +(y +1) 2 (19)
So, the reflected field shrinks to a single image-source located at x =0 , y = -1. Of course, this is a well-known property of rays. But the fact that the formation of image-sources is representable by an integral of type R σ/2⇡cRd⌫ points out that Eq. ( 1) accounts for a wide class of phenomena in spite of the fact that interferences have been neglected.

The extension of the radiosity method to specular reflection has been studied for a long time in thermics [START_REF] Eckert | Radiative heat exchange between surfaces with specular reflection[END_REF][START_REF] Sparrow | An enclosure theory for radiative exchange between specularly and diffusely reflecting surfaces[END_REF][START_REF] Modest | Radiative heat transfer[END_REF] and more recently in acoustics [START_REF] Kuttruff | Stationary propagation of sound energy in flat enclosures with partially diffuse surface reflection[END_REF]. The classical solution is based on the image-source technique and the introduction of specular view-factors. But some non-standard solutions have also been proposed recently [START_REF] Franzoni | An acoustic boundary element method based on energy and intensity variables for prediction of high frequency broadband sound fields[END_REF][START_REF] Cotoni | Specular and diffuse reflections of rays in coupled thin plates at high frequencies[END_REF]. In statistical energy analysis, the image-source technique also provides a powerful extension of the method when the sound field is not diffuse in the enclosure [START_REF] Sgard | Prediction of the acoustical performance of enclosures using a hybrid statistical energy analysis: Image source model[END_REF].

Multiple diffraction

Another interesting example of problems to which the representation equation (1) may give an elegant solution is that of multiple diffraction [START_REF] Reboul | Introduction of acoustical diffraction in the radiative transfer method[END_REF][START_REF] Reboul | Radiative transfer equation for multiple diffraction[END_REF]. To illustrate the problem, let us consider a two-dimensional acoustic medium with an incident stationary plane wave on two small obstacles localized respectively at the origin and (0,L) as shown in Fig. 4.

To solve this problem, two fictitious sources λ 1 and λ 2 are set to the diffraction points. Their directional emissive powers λ i /2⇡ are given by Eq. ( 9) where the source term is 2⇡I 0 and no source σ has been placed. Then

λ 1 (✓)=D(', ✓)2⇡I 0 + D(0,✓) λ 2 (⇡) L ( 20
)
λ 2 (✓)=D(', ✓)2⇡I 0 + D(⇡, ✓) λ 1 (0) L (21)
Eqs. [START_REF] Guyader | Integral equation for frequency averaged quadratic pressure[END_REF], [START_REF] Herrin | The energy source simulation method[END_REF] give the emissive powers λ 1 , λ 2 in any direction provided that they are known in grazing directions. We must first determine λ 1 (0) and λ 2 (⇡). Applying the two above equations respectively to ✓ =0a n d✓ = ⇡ The source is at infinite with incidence ϕ, the receiver at (x, y) and the diffracting points at (0, 0) and (0,L). The emission angle is θ.

gives,

λ 1 (0) = D(', 0)2⇡I 0 + D(0, 0) λ 2 (⇡) L (22)
λ 2 (⇡)=D(', ⇡)2⇡I 0 + D(⇡, ⇡) λ 1 (0) L (23)
This set of linear equations is readily solved, giving

Finally, the energy density at any point is found by applying Eq. (1),

W (x, y)= I 0 c + λ 1 (✓ 1 ) 2⇡c p x 2 + y 2 + λ 2 (✓ 2 ) 2⇡c p (x -L) 2 + y 2 (26)
This is the solution to the problem of double diffraction.

Transition to statistical energy analysis

When the field becomes diffuse in a domain, a simple law emerges: The energy leaving the system through the boundary is proportional to the total energy,

P 1!2 = η 12 ωE 1 (27)
where E 1 is the total energy, P 1!2 the outgoing power, ω the frequency (rad/s) and η 12 the so-called coupling loss factor. When two subsytems are connected by a common boundary, the net exchanged energy may be decomposed as a difference of two terms

P 12 = P 1!2 -P 2!1 = ω(η 12 E 1 -η 21 ωE 2 )
. This is the coupling power proportionality law well-known in statistical energy analysis. This phenomenon is correctly predicted by the radiative transfer equations. Assuming for convenience a unique stationary point source located at x 0 and no diffracting point, the total energy is obtained by integrating the energy density W (x)o fE q .( 1 )o v e rt h ed o m a i n ,

E 1 = Z Ω dΩ x 4πcR 2 0 + Z Ω Z Γ σ (y, v) 4πcR 2 dΓ y dΩ x (28)
where R 0 = |x 0 -x| is the source-receiver distance and R = |x -y|.T h e outgoing power is the flux of the normal component of the incident radiative intensity I(x, v)cosϕ,

P 1!2 = α ✓Z Γ cos ϕ 4πR 2 0 dΓ x + Z Γ Z Γ σ (y, v) cos ϕ 4πR 2 dΓ y dΓ x ◆ (29)
Let us consider the canonical problem of a system of volume V lightly damped with reflecting boundaries excepted a small aperture of area S of totally absorbing boundary. We may wonder if the ratio P 1!2 /E 1 tends to the asymptotic limit,

η 12 ω = Sc 4πV α (30)
predicted by statistical energy analysis. In Table 2 is shown the evolution of P 1!2 /E 1 and η 12 ω versus S for a room of 20 × 30 × 5m 3 with an aperture centred on a wall 30 × 5. The aperture has α =0 .99 while α =0f o ro t h e rs u r f a c e s . As i n g l epo i n ts o u r c e of 1 W is located at the centre of the room. To perform the calculation, the radiosity equation ( 12) has been solved in B(y)i ns t e a d y -s t a t ec o n d ition by a collocation method with a mesh of about 1200 triangles. E 1 and P 1!2 have been calculated by Eqs. [START_REF] Kuttruff | Simulierte Nachhallkurven in Rechteckräumen mit diffusem Schallfeld[END_REF][START_REF] Kuttruff | Room acoustics[END_REF] where σ (y, v)=4 B(y)cosθ have been substituted. The number of receiver points to assess E 1 is 36. The ratio of standard deviation and mean value of energies also appears in Table 2. The results clearly indicates that a large discrepancy between η 12 and P 1!2 /E 1 is observed for large apertures where the direct field is dominant aperture (m 2 ) 28 × 4 10 × 4 4 × 4 2 × 2 1 × 1 η 12 ω (s -1 ) 3.1416 1.1220 0.4488 0.1122 0.0280 P 1!2 /E 1 (s -1 ) 1.8220 0.8624 0.3942 0.1094 0.0282 std/mean (%) 126 58.8 26.9 7.45 1.91

Table 2: Comparison of the empirical coupling loss factor P 1→2 /E 1 with the asymptotic limit η 12 ω of statistical energy analysis for a room of size 20 × 30 × 5m 3 with perfectly reflecting (α = 0) walls, floor and ceiling and an aperture (α =0 .99) centred on the wall 30 × 5. The point source is centred in the room.

(large std/mean). But the convergence is good for small apertures when the field becomes diffuse inside the room (small std/mean).

Conclusion

In this study, the theoretical framework for a generalization of the radiosity method in acoustics and statistical energy analysis in vibro-acoustics has been presented in detail. The main features are the followings. The vibrational field, described in terms of energy density and specific intensity, is constructed as a linear superposition of actual and fictitious sources positioned in the domain, on the boundary and at the diffraction points. The emissive power of these sources is driven by the balance of radiant exchanges as in optics or thermics.

The method whose principle is well-known in thermics seems to be also useful in vibro-acoustics. The development of the radiosity method in acoustics and its application to sound illumination problem in room acoustics or reverberation time determination beyond the validity of Sabine's law underline its importance in engineering. Furthermore, its extension to specular reflection and diffraction highlights that it may be also useful in a larger class of problems, including sound radiation, structural response, sound transmission [START_REF] Bot | Energy exchange in uncorrelated ray fields of vibroacoustics[END_REF]. From this point of view, radiative exchange equations look like an a t u r a le x t e n s i o no fs t a t i s t i c a le n e r g ya n a l y s i sw h e ne n e r g yfi e l d sa r en o t diffuse.

Figure 1 :

 1 Figure 1: Actual and fictitious sources ρ, σ and λ.

Figure 2 :

 2 Figure 2: Specific radiative intensity from a surface source.

Figure 4 :

 4 Figure 4: Double diffraction.The source is at infinite with incidence ϕ, the receiver at (x, y) and the diffracting points at (0, 0) and (0,L). The emission angle is θ.

λ 1 (

 1 0) = 2⇡I 0 D(', 0) + D(0, 0)D(', ⇡)/L 1 -D(0, 0)D(⇡, ⇡)/L 2 (24) λ 2 (⇡)=2 ⇡I 0 D(', ⇡)+D(⇡, ⇡)D(', 0)/L 1 -D(0, 0)D(⇡, ⇡)/L 2

Table 1 :

 1 Reverberation time (in s) for a rectangular room of size 20×30×5m 3 by standard formula (Sabine, Eyring, Millington) and the radiosity equation (

	14)

Acknowledgements

The authors sincerely acknowledge the Labex Celya (Lyon acoustic centre) in the framework of which this work was conducted.


Marie -Josée J Pac

S Giljean

Christophe Rousselot

Fabrice Richard

Patrick Delobelle
email: patrick.delobelle@univ-fcomte.fr

Microstructural and elasto-plastic material parameters identification by inverse finite

Keywords: Nitride films, nano-indentation, residual stress, inverse method, identifiability, finite elements, fracture toughness

evolution, a simple model taking into account the stiffnesses coefficients of TiN and AlN structures, the mean residual stress level and the variations of the lattice parameters in the two structure domains is proposed.

Knowing the elastic properties of these films, inverse finite elements analysis of the indentation curves considering a simple isotropic linear elasto-plastic behavior allows, as a function of the composition, the yield stress Y σ and the linear hardening coefficient * p H to be estimated. Y σ and * p H are in the ranges 4.2 to 6.8 GPa and 60 to 400 GPa, respectively. The maximum value of Y * p / H σ which characterizes the ability of these coatings to exhibit plastic strain hardening is maximum for x= 0.5 and 0.6. The quality of the estimation was discussed through a practical identifiability study and quantified using an identifiability index. Tip radius and elasticity of the Berkovich indenter are two very relevant parameters to improve identifiability and correctly extract the plastic parameters of the behavior law. Scratch crack propagation resistance shows a evolution similar to those of Y * p / H σ .

-Introduction

Among the metallic nitride used as protective coatings in order to improve the mechanical resistance of steel surfaces for various machining applications, Titanium and Aluminum nitrides have extensively been studied [1][2][3] and can be considered as model system. Hence, Ti 1-x Al x N ternary nitrides crystallize in the NaCl cubic structure (fcc) for the low Al contents and in the ZnS hexagonal Wurtzite structure (hcp) for the high Al contents. The transition between these two structures occurs around x = 0.65 depending on the elaboration conditions [2][3][4][5][6][7][8]. The hardness of Ti 1-x Al x N thin films increases with increasing Al content up to about the transition content and then a rapid decrease is observed [4][5][6], [9]. This is due to the structural decomposition and the formation of a mixed structure (fcc + hcp). Moreover, the replacement of Ti atoms by Al atoms increases the oxidation resistance at elevated temperature [2][3][4]. Therefore, the best mechanical properties combined to the best resistance to oxidation is obtained for the maximum Al content substituted to Ti in the fcc structure.

However, recent studies [6][7], [10][11][12] show that the micro and nano structures of films are much more complex. The low solubility of Ti atoms in hexagonal structure gives rise to segregation between Al-rich hexagonal domains and Ti-rich grain boundaries when the Al content exceeds 0.5 in the ternary system. Note that the role of this phase segregation under thermal treatments on the wear resistance has also been evidenced [3], [6][7], [10].

In our previous papers [8], [11][12][13], the structures at different scales (macro, micro and nano) of Ti 1-x Al x N films elaborated by physical vapor deposition radio-frequency sputtering have extensively been studied by X-ray diffraction (XRD), X-ray absorption fine structure (EXAFS) performed at both Al and Ti K edges, scanning and transmission microscopies (SEM and TEM). As a complement of these different papers the present investigation focuses on the mechanical properties of such deposited films. Hence, with respect to the Al content (0 < x < 1), the different Ti 1-x Al x N films deposited on Si and steel substrates were examined by X-ray diffraction to determine the microstructure parameters (growth directions, lattice parameters, FWHM (Full Width at Half Maximum) of diffraction peaks), by atomic force microscopy (AFM) to measure the root mean square roughness (R ms) and the morphology of the surface of the films, by four-probe electrical resistivity method, by mechanical deflection of cantilever beams to estimate the mean internal stress and by nano-indentation to determine the indentation modulus and the hardness of the films. The plastic deformation parameters (yield stress σ y and work hardening coefficient H * p ) will be determined for each Al content by inverse finite elements method of the indentation curves. The quality of the estimations will be discussed through a practical identifiability study and quantified using an identifiability index.. Note that the quantification of these two parameters has never been reported in the literature.

-Experimental procedures

Sample fabrication

Titanium aluminum nitride Ti 1-x Al x N films have been deposited by radio frequency magnetron reactive sputtering in Ar-N 2 reactive mixture from a composite metallic target [14] onto both Si (100) and stainless steel substrates [8], [14]. The steel substrates come from a thick plate made of high speed steel (X90WMoCrV6.5-4.2) which is commonly used in cutting tool manufacturing. After machining the blocks have been quenched and polished.

Target materials (diameter ∼ 50 mm) were made from metallic titanium and aluminum and powered by a rf generator (13.56 MHz) at 80 W. The target was divided into 12 even parts in order to change the S Al / (S Al + S Ti ) area ratio allowing the adjustment of the (x) Al content [14]. Thus seven compositions (0 < x < 1) have been elaborated (Table I). The substrate holder was water cooled at T s = 298 K. The substrates have been cleaned with acetone and alcohol and pre-sputtered in pure Ar atmosphere with a negative bias voltage of -220 V for 40 min. The pressure reached before the introduction of the gas was 10 -5 Pa and the Ar and N 2 mass flow rates were controlled and adjusted to keep a constant total sputtering pressure of 0.65 Pa. The Ar pressure was kept constant at 0.55 Pa using a mass flow rate controller and a constant pumping rate of 10 Ls -1 . The N 2 pressure was fixed at 0.1 Pa to run the process in the nitrided sputtering mode. More details on the deposition parameters are given in refs. [8][11].

With these experimental conditions and as reported in Table I, the growth rate V d is an increasing function of the Al content. V d increases from about 200 nmh -1 for TiN to about 680 nmh -1 for the AlN. According to the characterization methods, the deposition time has been adjusted in order to prepare three sets of films: thinner films of about 0.3 µm for XRD analysis, films of about 1 µm thick for the stress analysis and the AFM measurements and thicker films of about 3 µm for the nano-indentation experiments and the AFM measurements. Table I gives the exact values of the film thicknesses with respect to the composition and the nature of the substrate. The composition of the specimens has been determined from an energy dispersive spectrometer and the Al content is higher that the one given by the S Al /S Al +S Ti area ratio as shown by Rauch et al. [14] (Table I).

Experimental methods

The crystallographic structure of the films has been investigated by using a Panalytical X'pert PRO diffractometer equipped with a prefix hybrid device and the Cu Kα 1 radiation (λ = 0.154 nm). The X-ray diffraction patterns have been recorded at room temperature in ω/2θ parallel beam geometry with ω = 10°. A real time multiple strip-type detector has been used and the accumulation lasted about 900 s per step of 0.17°.

The exact thickness e f of the films (Table I) has been measured with a stylus surface profiler DEKTAC 3030.

The electric resistivity has been measured by the classic four-probe method at room temperature.

Si(100) cantilever structures correspond to a display with seven single beams of 0.5 mm width and different lengths in the range of 0.5 to 3.5 nm. This display allows us to measure with an optical microscope the deflection at the extremity of the beams due to the mean residual stress into the sputtered film. Note that the thickness e s of the Si beams has been determined by a direct observation of the cross section thanks to a scanning electronic microscope (SEM JEOL 2100FCs); 12 < e s < 16 µm.

The surface morphology and the R ms (root mean square) roughness of the films have been investigated with an atomic force microscope (AFM -PSIA XE 150) in a non-contact mode. Different areas of 1 x 1 µm 2 and 4 x 4 µm 2 were scanned at three different positions in the surface in order to make sure that the images and the roughness taken in AFM are representative of the film surface. Some TEM observations of the cross section of the films after focused ion beam milling (Philips FIB 200) have been realized.

Nano-indentation tests have been performed on the Ti 1-x Al x N films of 3 µm thick using a nano-indenter II s (Nano-Instrument) with a Berkovich tip. To accurately determine the elastoplastic properties of the deposited films thanks to inverse finite elements modeling a blunt Berkovich tip with a R radius of 550 ± 50 nm has been used [15][16][17][18][19][20][21][22][23][24][25]. Indeed, as further shown, the tip radius and the deformability of the indenter are very important parameters to accurately determine the inelastic properties of thin films thanks to inverse identification of nano-indentation curves. Note that, for a perfect self similar Berkovich tip (R<100 nm) only one point of the stress-strain curve could be identified [15][16][17][18][19][20]. The study was conducted following the continuous contact stiffness measurement procedure (CSM) [26] with a frequency of 45 Hz and an indenter vibration amplitude of 1 nm during the penetration of the tip into the sample. For each tested film the measurement sequence consists on 6 indents with a 50 µm space between them on three or four different zones of the film (18 or 24 indents), with a maximum load of P max = 160 mN. The penetration speed was not constant but increased with depth from 0.3 to 2.1 mNs -1 . Note that some complementary experiments have also been carried out using a CSM-NHT instrument equipped with a Berkovich tip. In this case the tip radius of the indenter is of the order of 100 ± 20 nm.

-Structural and electrical characterizations of the films

X-ray diffraction analysis

The diffraction patterns of these thin films (e f ∼ 300 nm) deposited on Si substrate have already been presented in the paper of Tuilier et al. [11]. From the 2θ values of the maximum of the diffraction peaks and the full width half maximum (FWHM) values Δ2θ of these peaks, the growth directions, the lattice parameters (Bragg's formula) and the structural domain size (Scherrer's formula) have respectively been deduced [START_REF] Cullity | Elements of X-ray diffraction[END_REF]. As a function of the Al content, , respectively [START_REF]Powder Diffraction File, JCPDS International Center for Diffraction Data[END_REF]). This point will be further confirmed (section 4.1). The narrowness of the (002) peaks shows the good crystalline quality of these films. For the intermediate compositions (x = 0.6 and 0.68) a distorted hexagonal structure is highlighted. The (002) reflection is shifted toward the low 2θ values (34.45° and 34.84°, respectively) and an additional peak at 2θ = 32.15° assigned to the (100) reflection is equally observed. The peaks are significantly broadened with respect to the AlN and Ti 0.14 Al 0.86 N films. The large lattice parameters increase suggests an incorporation of Ti atoms in the hexagonal structure [11]. Moreover, an increase of the background centered at about 2θ = 37° reveals a weak contribution of the (111) oriented cubic domains. Indeed, as shown by Tuilier at al. [11] for x = 0.68 (Al K-edge EXAFS and XRD), hexagonal domains have grown while at the grain boundaries disordered cubic Ti rich clusters has been formed.

AFM measurements

The R ms root mean square roughness and the surface morphology of the different films of about 1 µm and 3µm thick have been investigated on 1 x 1 µm 2 and 4 x 4 µm 2 areas. A very fine structure in agreement with the columnar growth of this kind of deposit has been observed [13]. Moreover, a coarse-grained structure appears with the increase of the Al content. However, as shown in Table I, the thickness of the films is not constant but increases with the Al content and the principal change in the surface morphology could be attributed to the rise of the column diameter with the film thickness which is a common property of the growth of sputtered films. The surface morphology on 1 x 1 µm 2 area of 3 µm thick films exhibits very coarse grains, excepted for x = 0.6, which shows that the column diameter greatly increases with the film thickness. Fig. 2a shows, as a function of the Al content, the evolution of the R ms roughness for different scan areas: 0.25 x 0.25 µm 2 , 1 x 1 µm 2 and 4 x 4 µm 2 for the 1 µm thick films. For these films the 1x1 µm 2 and 4 x 4 µm 2 areas give the same R ms . R ms is of the order of 2.5-3 nm in the fcc domain and of 3.5 nm in the hcp domain, in a qualitative accordance with the AFM surfaces morphologies. For x = 0.6, in the fcc + hcp range a very small roughness of about 1 nm is measured. Over very small scan area, 0.25 x 0.25 µm 2 , excepted for x = 0 and 0.6 where a very fine surface morphology can be observed and consequently R ms is independent of the scan area, the roughness is smaller than for the 1 x 1 µm 2 area. Indeed, at this scale and as it can be observed on the AFM pictures a very fine microstructure can be evidenced. For the 3 µm thick films, on 4 x 4 µm 2 scan area, the roughnesses are greater than those measured on 1 µm thick films and also present a minimum in the fcc + hcp range. Hence, R ms is an increasing function of the film thickness. Note that, the standard deviations of the nano-indentation measures are directly related to the R ms parameter, as further described in (4.2.1.1).

The mean column diameter at the surface of the films can be evaluated from the AFM pictures. This has been done on the 1 x 1 µm 2 and 0.25 x 0.25 µm 2 areas for the 1 µm thick films and the results are presented in Fig. 2b: the mean column diameter φ is plotted as a function of the Al content. Moreover, for the 0.3 µm thick films (XRD, Table I) the structural domain size (or the nano-crystal size) d* was calculated from the FWHM of the diffraction peaks using the Scherrer's formula. The results are given in Fig. 2b and it is shown that the evolutions of φ (1 x 1 µm 2 ) and d* are parallel and present two maxima for x = 0.5 and 0.86 in the fcc and the hcp domains, respectively. Of course the φ and d* values are different. This is due to the application domain of the Scherrer's formula which is not entirely verified in the present case and to the difference in the film thicknesses.

However, the values of the diameters of the very fine crystallite structures determined on the 0.25 x 0.25 µm 2 areas are close to those of d* and in the range 5 to 40 nm. This is in agreement with the HRTEM results reported by Girleanu et al. [13] on the same films: d* is in the 5-30 nm range. The φ and d* values in the hcp domain are higher than those in the fcc domain and present a minimum in the fcc + hcp region. This evolution is consistent to the one observed on the roughness (Fig. 2a). The values of the column diameters φ (1 x 1 µm 2 ) are in the range 60-80 nm in the fcc domain and 80-120 nm in the hcp region. These values are close to those obtained by TEM direct imaging after FIB sectioning of 1.2 µm thick films: φ = 80 nm for TiN and Ti 0,5 Al 0,5 N, and φ = 60 nm for Ti 0,14 Al 0,86 N (Fig. 2b).

Electric resistivity of the films

The electric resistivity of the films deposited on the Si substrate was measured at room temperature on two or three different zones. The results are given in Fig. 3. Note that TiN conducts electricity like a metal but that AlN is a good insulator. Hence, the electric resistivity slowly increases with the Al concentration in the bcc structure (conductor), but as soon as a small amount of hcp structure (insulator) appears (x > 0.58-0.59, see Fig. 1) the resistivity sharply increases from about 66 10 -4 Ωcm for x = 0.5 to 460-600 10 -4 Ωcm for x = 0.6. For x > 0.6 the electric resistivity cannot be determined with our experimental device.

-Mechanical characterization of the films

Mean residual stress measurements

The analysis of the bending of the bilayer cantilevers (Si substrate + films, Table I) as a function of the Al content gave valuable informations about the mean stress levels in the film.

Note that the nano-indentation results (indentation modulus and hardness) are directly affected by the mean stress level. This will be further quantified (4.2). For a given display corresponding to a fixed Al content the deflections δ m at the extremity of the seven cantilevers of different lengths L (0.5 mm to 3.5 mm by step of 0.5 mm and a width of 0.5 mm) have been measured thanks to optical microscope.

Several theoretical models regarding stress analysis of single layer of thin film deposited on an elastically isotropic substrate have been developed since the first works of Stoney in 1909 [START_REF] Stoney | [END_REF][30][31][32]. The first order development of the general solution gives a generalized writing of the Stoney's equation:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ε + + ν - δ = σ ) 0 ( s s f f 2 f m 2 s s 0 e E e E 4 1 L e ) 1 ( 3 e E ( 1
)
Note that in this study, the use of the Stoney's equation (first term in Eq.( 1)) would have introduced a significant under-estimation of the stress in the film since the value of the term

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + s s f f e E e E 4 1
is in the range of 1.5 to 2.2. The E f values have been determined from the nanoindention tests on the different films (see below) and E s =165 GPa. Following the Eq.( 1) the representations δ m = f(L 2 ) (Fig. 4) should be linear and the slope of each straight line permits the mean stress σ 0 to be determined. All the films present a compressive stress state excepted for the AlN film which is in a tensile one (Fig. 4). Figure 5 gives, as a function of the Al content, the evolution of the mean residual stress. Some results obtained in previous unpublished work have also been reported. In the fcc area the stress weakly decreases from -2.7 GPa to -2 GPa with the Al content which is in agreement with the results reported in the literature [4] [7]. On the contrary, in the hcp domain the stress greatly evolves from -1.9 GPa for x = 0.86 to +0.5 GPa for x = 1. In the intermediate domain where the two phases coexist the stress values are relatively low of the order of -0.8 GPa for x = 0.6. A minimum value close to zero could be expected for x ~ 0.64-0.66. Hence, the stress variation presents a very significant minimum in the fcc + hcp region. As shown in Figs 2a-b, the roughness and the estimated superficial mean column diameter φ also present a minimum in the vicinity of x = 0.64. The low residual stress values could be attributed to the stress relaxation due to the grain boundaries sliding of the un-deformable small grains. This is supported, as previously mentioned, by the Al K-edge EXAFS measurements reported by Tuilier et al. [11] which

show that grain boundaries of disordered cubic Ti rich clusters has been formed for x = 0.68.

Note that for TiN and AlN films the mean stress levels (compressive and slightly tensile, respectively), are in quantitative agreement with the shifts of the XRD peaks.

Nano-indentation results

As mentioned in appendix A1, nano-indentation tests allow the indentation modulus M <hkl> and the Berkovich hardness H B<hkl> to be determined. These two quantities depend on the direction <hkl> of the domain growth. As the indentation modulus of the Si substrate (M S =M Si ~173 GPa, Table II) is lower than those of the TiAlN films, due to the substrate effect, the measured indentation moduli continuously decrease with h/e f . Moreover, for h/e f < 3 %, due to the spherical shape of the Berkovich indenter tip (R~550 nm), the M <hkl>meas moduli are increasing functions of the penetration depth [33]. The true modulus of the films can be obtained by extrapolating the measured values in the neighbourhood of h/e f = 0. Among the numerous analytical expressions to take the substrate effect into account [33][34][35][36], the King's model [34] has been used as the indentation modulus M <hkl> of the film is the only one unknown in this model.

Hence, the measured indentation modulus M <hkl>meas is given by:

n n 2 c c 1 c f hkl S hkl meas hkl a h 56 , 24 A with A e exp M 1 M 1 M 1 M ∑ + = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ α - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + = - > < > < > < (2a) .
α is a parameter which depends on the projected contact area and for a triangular Berkovich indenter, [34] :

2 / 1 f c e A 4 3 2 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ≈ α (2b)
M <hkl> and M S are the indentation moduli of the film and the substrate, respectively.

The King's model predictions over the entire h/e f range have been plotted in Fig. 6a. From these predictions the true indentation modulus without substrate effect (horizontal part of the curves) is obtained for h/e f < 2%, thus for h< 60 nm. The values for h = 0 have been reported in Fig. 6b and correspond to those of the indentation modulus of the films. The three domains highlighted by the XRD analysis (section 3.1) are perfectly identified. The value determined for x = 0.5 with the steel substrate is close to the one measured with the Si substrate. Note that the results obtained by linear extrapolation to h = 0 of the curves of Fig. 6a (linearization on the range h/e f = 3.3 to 6.6 %) are also reported in Fig. 6b. In the fcc domain M <hkl> is in the range 400-440 GPa, in the range 300-340 GPa for the hcp domain and of the order of 240-380

GPa in the fcc + hcp region for x = 0.6 to 0.68. As a comparison, the Kutschej et al. [6] results

have also been plotted in this figure and are in a fairly good agreement with those of the present study. This point will be further discussed in the paper. From a quantitative point of view and as mentioned in appendix A1, knowing the stiffness coefficients C ij of TiN and AlN, the indentation modulus for a given direction can be calculated. TiN is weakly anisotropic [37] and the C ij values taken into account [38][39] are listed in Table II. The values M iso and ν iso of isotropic film and evaluated thanks to the Hashin and Shtrickman formulation [40] have also been reported in Table II. Applying the Eq.(A1-6) [41] for the growth directions of the TiN film the values of M <100> = 450 GPa and M <111> = 440 GPa are obtained.

Experimentally, a value of M = 400 ± 20 GPa has been determined, thus weakly lower than the theoretical one. However, this can be attributed to the poor crystallinity of the TiN films as observed on the X ray diffraction patterns [11]. The films with x = 0.38 and 0.5 which present a good crystallinity [11] exhibit an indentation modulus of 430 GPa and 440 GPa, respectively. These values are in a fairly good quantitative agreement with the theoretical one (Table II), especially for the <111> direction which has been identified as a growth direction (Fig. 1).

For the AlN, the C ij coefficients of sputtered thin films [42] are listed in Table II. As previously done, the application of Eq.(A1-7) allows to M <001> = 342 GPa and M <100> = 319

GPa. For the <001> direction which is a growth direction of the AlN film the M <001> value is in a very good agreement with the experimental one: M <001> = 340 ± 10 GPa. Note that the effect of the residual stresses on the measured value of the indentation modulus is quite negligible (<1%) as σ 0 /H B < 10% (see appendix A2) [43][44][45]. Moreover, as soon as the scattering of the measures is concerned, the relative standard deviation ΔM/M mean is related to the dimensionless parameter R ms /h by the following relation [46][47]:

m ms mean h R M M ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ α = Δ (3)
with α = 0.35+-0.1 and m = 0.64. The R ms values have previously been reported in Fig. 2a.

Excepted for h = 50 nm, there is a fairly good agreement between the experimental values of the scattering and those given by the relation (3). As an example, for h=100 nm, ΔM/M mean is in the range 4-6 % over the entire range of Al content. Of course, ΔM/M mean increases when the indentation depth decreases and ΔM/M mean could be greater than 8 % (25-35 GPa) for h = 60 nm which is the indentation depth below which the indentation modulus is independent of the substrate effect (h/e f < 2%).

The variations of the indentation modulus in the fcc and hcp domains have been connected with the variations of the lattice parameters in these two regions. However, as previously mentioned (section 3.1) and due to the compressive residual stresses the lattice parameters a i (i= fcc or hcp) calculated from the XRD patterns are overestimated. Knowing the mean value of the residual stress σ 0 (section 4.1), the values of the isotropic indentation modulus and Poisson's ratio, M iso and ν iso , respectively (Table II), the true values a i0 of the lattice parameters can be estimated by: The indentation modulus increases with the decrease of the lattice parameter, thus:

) hcp ( or ) fcc ( i and 2 / , ) 1 ( M sin 1 a a iso iso 2 0 i 0 i = ω - π = ψ ν - ψ σ - = (4)
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - β - = > < > < ) TiN ( a ) TiN ( a a 1 ) TiN ( M ) fcc ( M 0 c 0 c 0 c c hkl hkl (5a)
and a similar expression for the hcp domain:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - β - = > < > < ) AlN ( c ) AlN ( c c 1 ) AlN ( M ) hcp ( M 0 h 0 h 0 h h hkl hkl (5b)
or with the lattice parameter a h0 (AlN) as ( )

) hcp ( M %) x ( 1 ) fcc ( M %) x ( M hkl hkl hkl > < > < > < α - + α = (6a)
with α(x %) = 1 for x % < 0.58 and α(x %) = 0 for x % > 0.7. The choice of the α(x %) function is somewhat arbitrary and an exponential function which verifies the boundary conditions has been considered. So: Hence, knowing the stiffness coefficients C ij of the TiN and AlN crystals, the preferential growth directions <hkl> of the crystallites and the lattice parameters (a c , a h , c h ) determined by XRD, the mean residual stress σ 0 in the film, the indentation modulus values can be estimated thanks to the previous model whatever the x(%Al) content is. Moreover the precise knowledge of the elastic properties of these films is absolutely necessary to correctly identify the plastic coefficients thanks to inverse analysis.

n % x 7 . 0 58 . 0 % x exp %) x ( ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ > - < > - < - = α (

Berkovich hardness H B<hkl> :

Figure 7a shows the variation of the measured Berkovich hardness H B<hkl>meas as a function of the reduced indentation depth. As for the indentation modulus due to the rounding of the indenter tip the hardness increases with the indentation depth until h/e f ∼ 4-5 % [33] and then slightly decreases due to the substrate effect (H S = H B (Si) ∼ 13 GPa < H B<hkl> film). To quantitatively take this effect into account and to determine the true hardness of these films, the Bhattacharya and Nix's model [48] for hard film on softer substrate has been considered.

Thus:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ σ σ - - + = > < > < > < > < f S hkl S f S hkl B S hkl B S meas hkl B e h E / E / H / H exp ) H H ( H H (7a)
where H B<hkl> and H S are the hardnesses, σ f , σ S the yield strengths, E <hkl> , E S the Young's moduli of the film and the substrate, respectively. In the absence of information on the yield strengths it has been assumed that H B<hkl> /H s ~ δ σ f /σ S with 0.5 < δ < 1. Hence Eq.(7a) is rewritten as:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ δ - - + = > < > < > < f S hkl S hkl B S meas hkl B e h E / E 1 exp ) H H ( H H (7b)
To fit the decreasing of the experimental curves (Fig. 7a) knowing E <hkl> , E S (deduced from the M <hkl> values with ν = 0.2) and H S (H Si ∼ 13 GPa) a value of δ = 1/2 has been determined.

Thus the yield strength ratio is twice the hardness one. Note that, due to the confinement of the plastic deformation field the decrease of the hardness with the indentation depth is less pronounced than the one of the indentation modulus related to the elastic deformation field (see Fig. 8). As a function of the Al content, the measured hardnesses H B<hkl>meas at h = 100, 300 and 600 nm have been reported in Fig. 7b. The extrapolated values for h = 0 correspond to the true hardness of the films. The films deposited on steel substrate (x= 0.5) appreciably present the same hardness as those deposited on Si substrate. As a comparison the results issued from the Kutschej et al. [6] work have also been reported.

The complete analysis of the indentation curves thanks to inverse finite elements modeling is presented in the next section.

5 Inverse finite elements analysis of the indentation curves

Method

In this section the identification of the plastic material parameters of the films thanks to the different indentation curves is described. The analytical solution for indentation of an elastoplastic work-hardening solid does not exist [15][16][17][18][19][20][21][22] and the analysis is based on inverse numerical simulation.

The indentation test has therefore been modeled by a 2D axisymmetric representation with the ANSYS finite element software. The semi-vertical angle of the conical indenter ϕ = 70°3

gives the same contact area to depth ratio than a perfect Berkovich pyramid [15-22][49]. The cone tip is smoothed by a sphere of radius R (R = 550 nm and 100 nm depending on the used indenter). A film of 3 µm thick on a substrate has been modeled (see Fig. 8). Indenter size and substrate thickness have been determined to be insensitive to far-field boundary conditions: about 40 times greater (∼ 25 µm) than the maximum indentation depth (h max ∼ 600 nm).

The augmented Lagrangian method runs the contact algorithm between the indenter and the film (Coulomb friction coefficient equal to 0.2 [50]). In each computation and at the maximum load, the minimum number of contact elements in the contact zone is not less than 16. Note that for this indenter geometry the friction coefficient does not have any influence on the numerical responses [17] [50]. An implicit method with non-linear geometric formulation has been used with higher order 2D elements. The mesh density has been refined in order to the smallest elements are localized in the contact region. A convergence analysis has been performed in order to avoid any mesh influence.

The indenter (E i = 1100 GPa, ν i = 0.07) and Si substrate (E s = 165 GPa, ν s =0.23, see Table II)

have been considered as isotropic linear elastic materials. The case of rigid indenter has also been considered. For the film behavior model and for simplicity an isotropic linear elastoplastic model has been considered. A very important point: the Young's modulus of the films E(x) is deduced from the M <hkl> (x) experimental indentation modulus values considering a mean Poisson's ratio of ν = 0.21 (Table II). The considered loading surface is: , an inverse method has been used. The details of the procedure (Eq.A3-1), the method to calculate the uncertainties (Eq.A3-4) and the quantification of the quality of the estimations thanks to an identifiability index ) θ ( Î (Eq.A3-5) are given in Appendix A3. In order to illustrate the tip rounding, the deformability of the indenter and the substrate effects on identifiability, the following tests have been performed. ℑ has been calculated for x = 0.5 with R=500 nm, but also for three virtual cases: sharp and blunt indenters, rigid and elastic indenter associated with bulk and thin specimen (film thickness = 40 h max ). Table III shows that the identifiability is significantly improved when compared with sharp indenter (12%), rigid indenter (24%), but that for the considered thickness the substrate effect does not affect the index. So, the elasticity effect and the tip radius of the Berkovich indenter are shown to be two very relevant parameters to correctly identify the plastic parameters of the behavior law.

p * p s y s VM d H dR with R f ε = σ - - σ = (8).

Results

All sets of randomly chosen initial input parameters

Complete identification procedure for the set of Al contents has been conducted for rigid and elastic indenters. Figures 8 and9 give, for a deformable indenter, examples of elastic and plastic von Mises strain fields (x=0.5) and of the fit of the loading experimental curves by the numerical model (x= 0, 0.5, 0.86 and1), respectively. The elastic deformation of the indenter on one hand and the substrate and the confinement of the plastic deformation for the maximum loading (P max =160 mN) on the other hand are clearly highlighted in Fig. 8.

Simulation results (examples in Fig. 9) show a very good agreement with the experimental results, excepted at x = 0 and x = 1 where a fairly good agreement is obtained (see the errors in Table IV), probably due to local cracking in these less ductile materials. Figure 10 H values which are about 3 to 4 times lower than those calculated with an elastic indenter. However, their evolutions with the Al content are quite similar. These different points are discussed in the following section.

About the relationship between plastic parameters, hardness and rigidity of the indenter:

Considering the generalized Tabor's relation [START_REF] Tabor | The hardness of metals[END_REF], the corresponding stress σ (Δεpm) should be equal to:

3 / H H 0 B pm * p Y ) pm ( = ε Δ + σ = σ ε Δ (9)
where H B0 is the true hardness at zero stress (see appendix A2-1). If the indenter is rigid and the Berkovich tip perfect, the representative point of the stress-strain curve determined by inverse method should be at Δε pm = 3.3 % [15][16][17][18]. However, due to the tip rounding and the indenter stiffness, the representative point evolves with these two parameters [16] [START_REF] Tabor | The hardness of metals[END_REF]. To approximately determine the mean value of the representative strain Δε pm for the used blunt Berkovich tip (R=550 nm), a perfect correlation between the Tabor's stress and the flow rule is argued (Eq.9). Knowing for the set of composition, mean values of Δε pm ∼ 4.2 ± 0.5% and 1.1 ± 0.2% have been obtained for rigid and elastic indenters, respectively.

Note that the first value is in the range 3.3 % to 6.8 % corresponding to a perfect Berkovich tip and a spherical cap at the tangent point to the conical shape (

R / h 2 2 . 0 P ≈ ε Δ ),
respectively. For the two considered cases, as Y σ is practically independent of the indenter stiffness, the ratio ~3.8 between the two representative strains is approximately equal to the ratio between the hardening coefficients (~3 to 4).

For these very hard coating films, the elasticity of the diamond of the indenter greatly affects the values of the hardening coefficients * P H determined by inverse analysis and the representative plastic strain for our blunt Berkovich tip is about 1%, thus three to four times lower than the expected one for a perfect rigid tip.

-Discussion and analysis with regards to the literature

An interesting result which has never been directly evidenced is illustrated in Fig. 11 with the Y * P / H σ ratio, is the capacity as a function of the Al content of the film to exhibit plastic strain hardening. This ratio is maximum for x= 0.5 and 0.6:

70 60 / H Y * P - ≈ σ
. The same tendency has been obtained with a rigid indenter but in a ratio of 3-4 (Fig. 11). Indeed, this observation should be linked to the film toughness. The toughness of 3 µm thick films deposited on steel substrate by the same sputtering procedure as the one of the present study have recently been reported by Henry et al. [START_REF] Henry | [END_REF]. On these films, the "Scratch Crack Propagation Resistance" (CPRs) criterion proposed by Zhang et al. [57] has been evaluated and the values reported in Fig. 11.

Y * P / H σ and CPRs follow similar evolutions. Contrary to the Al rich films which present a weak resistance to the crack initiation, particularly for the AlN, the films with fcc structure exhibit a better fracture toughness. Note that the same conclusion could be obtained from the fracture toughness K 1c evaluated thanks to radial crack length measurements after Vickers indentations [56][58-59]. TEM studies [12] on crosssection of TiN milled by focused ion beam revealed that the microstructure of TiN is made of small grains clustered in vertical columns whereas Al rich films like for x = 0.86 consists of long fibers well textured along the <002> direction of hexagonal structure. Hence, the propagation of cracks within the film is thus easier along fibers in hcp domain whereas it is better impeded by a structure composed of small grains like in fcc region. The aptitude of the films to support plastic deformation follows the variation of the toughness.

Concerning the indentation modulus M and the hardness H, the presented results are close to those reported by Kutschej et al. [6] for dc magnetron sputtered films at relatively low temperature (T d ~350°C) and agree with those of Hörling et al. [7] after an annealing at 1050°C of the as deposited film to suppress the stress build-up created during the film deposition by the lattice point defects. The critical composition of the sharp decreasing of the mechanical properties is located in the 0.62 to 0.69 range (Kutschej et al. [6]: 0.67 -0.69, Hörling et al. [7]: 0.66 -0.67 and the present study: 0.62 -0.66) and are in a fairly good agreement with the theoretical one [2]: x= 0.65 ± 0.05.

-Conclusions

The mechanical properties of Ti 1-x Al x N (0 < x < 1) films of different thickness deposited by r.f. reactive magnetron sputtering on Si <100 > and high speed steel substrates have been investigated. The coatings with x < 0.58-0.59 present cubic structure whereas for x > 0.7 a hexagonal structure is observed. Between these two limits the two structures coexist. The p / H σ which characterizes the ability of these films to exhibit plastic strain is obtained in the fcc domain and fcc domain with a small amount of hcp phase, for x = 0.5 and 0.6. The quality of the estimations was discussed through a practical identifiability study and it has been shown that tip rounding and elasticity of the indenter offer informations which combined themselves improve the identifiability (+24%).

As the conclusion and according to the different published works on the mechanical properties of these ternary coatings, the optimum properties are located in the 0.5 to 0.6 composition range.

Appendix A1

In nano-indentation procedure, the M <hkl> modulus measured on a single crystal (or on a well textured material with very fine grains) whose normal to the surface has director cosine α i is given by [60][61]:

[ ] 1 k ij 1 km m 2 hkl d ) ( C 16 M 2 0 - - > < ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ γ α γ β α ∫ π = π with dh dP M and E 1 M 1 M c A 2 r 1 ind 2 r hkl β π - > < = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - - = (A1-1a)
For a Berkovich indenter the projected contact area A c is given by:

S P h h with ) h a ( h 56 , 24 ) h ( A c n / 1 c n 4 1 n 2 c c c ε - = ∑ + = = (A1-1b)
In these relations P is the applied load, S = dP/dh the unloading stiffness measured at the depth of penetration h, β = 1.034 and ε = 0.72 for the Berkovich'tip [19]. (E/(1-ν 2 )) ind is the reduced modulus of the diamond indenter. γ is the angle in surface plane and β mk a very complicated matrix functions of the elastic stiffnesses C ij (γ) [60] [61]. No general exact analytic solution of Eq.(A1-1a) exists. Nevertheless, for solids with cubic symmetry, Vlassak and Nix [60] showed that, using a numerical modeling, the following equations can be written:

M - = - + = β ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ν - β = > < > < > < (A1-2)
(E/1-ν 2 ) iso is the isotropic reduced modulus calculated thanks to the Hashin and Shtrickman [40] relations for homogenized materials. Coefficients a, c, f o and n depend on the Poisson ratio ν <100> and on the crystal orientation <hkl>. f ani is the anisotropic coefficient. Note that for isotropic materials, f ani = 1 which makes β <hkl> = 1 and, as a result,

. 1

E M iso 2 hkl ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ν - = > <
This formulation is usually used in the literature and Eq.(A1-1a) is recasted as:

1 ind 2 r iso 2 E 1 M 1 1 E - ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ν - - = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ν - (A1-3)
For orthotropic materials, good analytical approximated solutions are given by Delafargue and Ulm [62]. Hence, for the hexagonal structure the indentation moduli in the <001>, <110>

and <100> directions are written as: Based on the works of Delafargue and Ulm [62] and for cubic and hexagonal structures, Delobelle et al. [41] have proposed a general approximated solution applicable whatever the orientation is. Hence,

> < > < > < > < - = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + - = 001 33
[ ] γ α γ α π = > < > < > < = ∫ π > < d ) ( C 2 1 C with C C M M k mk ij m * ij * 33 * 11 iso hkl 2 0 (A1-5)
As for the Eq.(A1-2) M iso = (E-1/ν 2 ) iso is the isotropic reduced modulus and > < * ij C represent the mean values in the indentation plane (integration following the γ angle) of the elastic stiffnesses written in the indentation direction.

TiN is weakly anisotropic [37] and the C ij values taken into account [38][39] are listed in Table II. The values of M iso and ν iso evaluated thanks to the Hashin and Shtrickman formulation [40] have also been reported in Table II. Applying the Eq.(A1-5) for the growth directions of the TiN film with:

> < + + = > < + + = > < > < + + = > < = > < (A1-6)
the values of M <100> = 450 GPa and M <111> = 440 GPa are obtained. The application of the Vlassak et al. model [60] (Eq.A1-2) appreciably gives the same values.

For the AlN, the C ij coefficients of sputtered thin films [42] are listed in Table II

> < + + + = > < = > < > < = > < = > < (A1-7)
and M iso =304 GPa allows to M <001> = 342 GPa and M <100> = 319 GPa. The application of the Delafargue and Ulm's model [62] (Eq. A1-4) appreciably gives the same values.

Appendix A2

As shown by different authors [43][44][45] the measured values of the hardness H B and of the indentation modulus M depend on the level of the residual stress beneath the indented surface.

To take this effect into account the relation (A2-1) and which has been verified by F.E. analysis [45] has been proposed by Suresh et al. [34]: poor identifiability [START_REF] Gujarati | Basic Econometrics, Economics Series[END_REF]. For a value between 2 and 3, acceptability depends on the expected accuracy of the results.

0 if H sin 1 H f 0 if H 1 H f : with H f H H 0 1 0 B 0 0 B o 0 1 0 B o 0 B o 0 B o 0 B B < σ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ α σ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ σ > σ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ σ + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ σ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ σ = - - (A2-
Considering the indentation analysis a 2.5 value appears to be a critical threshold below which the identification can be considered as acceptable [25]. Nevertheless, whatever the critical threshold, the variations of ) θ ( Î can be studied to quantify the effects of the geometrical parameters of the F.E. model on the identifiability improvements due to the tip rounding, the elasticity of the indenter, the substrate effect and the data used on the indentation curve (loading/unloading parts). An example is given in the 5.2 section of the text. f ani 0.854 [38] 0.794 [39] 1.56 --M iso , ν iso 444 , 0.197 [38] 447 , 0.203 [39] 172.8 , 0.23 M iso (GPa) 304 M <100>

Figure captions

(GPa) 450 [38] 456 [39] 166 M <001> (GPa) 342 M <111>

(GPa) 439 [38] 440 [39] 177 M <100> (GPa) 319

Fig. 1

 1 Fig.1summarizes the evolution of the lattice parameters (a c : lattice parameter of the cubic

4. 2 . 1 . 1

 211 Figure 6a represent for each tested Al content the evolution of the mean measured values (18 or 24 indents) of the indentation modulus M <hkl>meas versus the reduced indentation depth h/e f .

.

  As an example, for TiN, with a c<100> = 4.27 o A , σ 0 = -2.7 GPa, M iso = 444 GPa, ν iso = 0.2 and ω = 10°, the previous relation gives a c0 = 4.239 o A , which is close to the tabulated value a c0 = 4.24 o A . For each composition the a c0 values have been calculated.

  β h are the two model parameters. For the intermediate compositions the two structures coexist (fcc + hcp) and a simple mixture law has been considered. Hence:

  6b) where < x > = x H(x) with H(x) the Heaviside function. The exponent n is a model parameter. The complete model for the evolution of the indentation modulus as a function of the Al content is constituted by Eqs. 5a-b and 6a-b. The identified model with β c = 1.0, β h = 9, n = 2, M iso (TiN) = 444 GPa (as <100> and <111> orientations have been observed), M <001> (AlN) = 342 GPa, a c0 (TiN) = 4,24 o A and c h0 (AlN) = 4,98 o A has been plotted in Fig. 6b.

  Mises stress and the equivalent von Mises plastic strain increment, s R the radius of the loading surface (isotropic hardening variable), Y σ the yield stress and * P H the linear hardening coefficient. Note that the isotropic elastic behavior of the films is perfectly justified for fcc phase but less for the hcp phase. Hence, as the Young's moduli and the Poisson ratio are known, the two only unknowns of the model are the yield stress Y σ and the linear hardening coefficient * P H . It is well known that indentation methodologies suffer from the non-uniqueness and sensitivity of the determined material properties [18-24][51-53]. It has been shown that the phenomena are not independent, rather non-uniqueness is an extreme case of sensitivity [24]. In this work, the quality of the estimations is evaluated through an identifiability index and it has been shown that the use of a non-rigid indenter in the model and tip rounding improve the sensitivity of the determined material parameters.To identify the optimal set of parameters )

  0 θ calculated by the genetic global algorithm using the loading curve lead to the same θ solution determined by the local algorithm[54]. Hence, the numerical uniqueness is considered as guaranteed. Whatever the Al content (for R = 550 nm and 100 nm) the practical identifiability analysis shows that the mean value of the identifiability index ) θ ( ℑ does not decrease during the unloading part of the indentation curves. Therefore the unloading part of the curves is not necessary as the Young's moduli of the films are known. On the contrary, the mean value of )

  shows, as a function of the Al content and for rigid and deformable indenters, the evolution of the optimal values of Y σ and * P H . For the legibility of the Fig.10 the uncertainties θ Δ calculated using Eqs.(A3-1 and 4) (elastic indenter) determined from the experimental loading curves and corresponding to the tests carried out with R = 550 nm and 100 nm are given in

  with a sharp indenter). Such high uncertainties in the determined material parameters have been reported by previous researchers [24][51-53]. Although the identified curves perfectly match the experimental ones (see the errors in table IV), the uncertainties on the plastic coefficients are relatively large as they are linked to the low convexity of the function to minimize. For rigid indenter the values of Y σ and * P H have also been reported in Fig.10. The Y σ values are appreciably independent of the indenter stiffness, contrary to the * P

  cubic structure shows a (200) c and (111) c preferred crystallization orientation whereas a(002) h orientation is observed for the hcp structures. The roughness and the superficial mean column diameters φ depend on the Al content and on the film thicknesses. A minimum is observed in the vicinity of x = 0.6. φ is in the 60-120 nm range whereas the crystallite diameters d* estimated by XRD and AFM are in the range 7 to 35 nm. As soon as small amount of hcp structure appears the resistivity sharply increases. Excepted for the AlN coating the residual stresses are compressive. As for the other measured properties the residual stress presents a minimum at the neighborhood of x~0.64. The indentation modulus M <hkl> and the Berkovich hardness H B<hkl> values greatly depend on the composition and a transition occurs when the two structures coexist for 0.58 < x < 0.7. To describe the M <hkl> = f(x %) evolution a simple model taking into account the C ij stiffnesses of TiN and AlN, the variations of the lattice parameters with the composition and the residual stresses, and the preferred crystallite orientations has been proposed. Inverse finite elements analysis of the indentation curves considering the tip rounding, the elasticity of the indenter and a simple isotropic linear elasto-plastic law whose elastic parameters are known (Young's modulus and Poisson's ratio) allows the yield stress Y σ and the hardening coefficient * p H to be extracted: Y σ and * p H lie in the ranges 4.2 to 6.8 GPa and 60 to 400 GPa, respectively. However, due to the low convexity of the functions to minimize, for a tip radius close to 500 nm the uncertainties on these two parameters are of the order of 40 to 70%, but drastically increase (up to 200 to 400%) for very sharp indenter (R~100 nm). Moreover, due to the great hardness of these films and the elasticity of the diamond of the indenter, the representative plastic strain is obtained at about 1.1% thus three times lower than the expected one. The maximum value of Y *

Fig. 1 :

 1 Fig. 1: Lattice parameters determined by XRD as a function of the Al content. Determination of the structural domains: fcc (a c < 200> and a c <111>), hcp (a h <002>, c h <002>) and fcc + hcp (a c < 200>, a h <002>, c h <002>).

Fig. 2 :

 2 Fig. 2: a)-Roughness R ms as a function of the Al content for different scan areas (4 x 4 µm 2 , 1 x 1 µm 2 and 0.25 x 0.25 µm 2 ) and two films thicknesses (e f ∼ 1 µm and 3 µm). b-Mean diameters φ of the top of the columns and mean diameters d* of the crystallites deduced from the AFM and SEM pictures and the FWMH of the XRD patterns.

Fig. 3 :

 3 Fig. 3: Electrical resistivity as a function of the Al content.

Fig. 4 :

 4 Fig. 4: Deflection δ m as a function of the square length L 2 of the cantilevers for each elaborated materials.

Fig. 5 :

 5 Fig. 5: Calculated residual stresses σ o as a function of the Al content. Results obtained in previous unpublished work (sputtering at T d = 350°C) have also been reported.

Fig. 6 :

 6 Fig. 6: a)-Evolution of the measured indentation modulus M <hkl> as a function of the reduced depth h/e f for the different studied films. Modeling with the King's relation Eq.(2). b)-Indentation modulus M <hkl> of the different tested films. Modeling thanks to the Eqs.(5-6). The Kutschej et al. [6] results have also been reported. Fig. 7: a)-Berkovich hardness H B<hkl> for the different tested films as a function of the reduced depth h/e f . Modeling with the Bhattacharya's relation Eq.(7b).b)-Hardness as a function of the Al content for three indentation depths (h=100, 300 and 600 nm). Modeling with the Battacharya's equation and determination of the true hardness H B<hkl> at h = 0. The results obtained with a steel substrate and those of Kutschej et al.[6] have also been reported.

Fig 8 :

 8 Fig 8:Example of elastic and plastic von Mises strain fields for R=550 nm, P max =160 mN and x= 0.5.

Fig 9 :

 9 Fig 9:Examples of experimental and calculated curves thanks to inverse method for x= 0, 0.5, 0.68 and 1.

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Yield stress σ Y and linear hardening coefficient H P * determined thanks to inverse finite elements analysis as a function of the Al content. Cases of rigid and elastic

  ) ( θ ℑfor x=0.5. The reference is the column present study TiN

Table IV .

 IV Note that no value of these two parameters are reported in the literature. The yield stresses

	Y σ and the hardening coefficients * P H lie in the range 4.2 to 6.8 GPa and 60 to 400 GPa,
	respectively. With R = 550 nm, the mean value of the relative uncertainty on the yield stresses
	calculated using Eq.(A3-4) is about	Δ	σ	Y	/ Y σ	≈	40	%	(	16	%	<	Δ	σ	Y	/	σ	Y	<	71	%	) and the
	relative uncertainty on the hardening coefficient is approximately 1.7 times greater:
	H Δ	* P	/	H	* P	≈	% 70	(	% 25	<	H Δ	* P	/	H	* P	<	% 150	). However, with R = 100 nm and x=0.6,
	Δ	σ	Y	/ Y σ	≈	90	%	and	H Δ	* P	/	H	* P	≈	% 300	which corresponds to a great loss of identifiability
	(see Table III for the

Table I

 I Thickness of the films e f (nm)

					Si<100> substrate	Steel substrate
	Composition X(Al%) S Al /(S Al +S Ti )	V d	XRD Stress	Nanoind	Nanoind.
				(nmh -1 )		AFM	AFM.
	Ti N	0	0	213	210	613	~3000±100	-
	Ti 0.62 Al 0.38 N	0.38	0.17	248	260	992	~3000±100	-
	Ti 0.5 Al 0.5 N	0.5	0.25	344	245	1012	~3000±100 ~3000
	Ti 0.4 Al 0.6 N	0.6	0.33	380	380	1007	~3000±100 ~2000
	Ti 0.32 Al 0.68 N	0.68	0.42	405	320	1565	~3000±100	-
	Ti 0.14 Al 086 N	0.86	0.67	522	270	1850	~3000±100	-
	Al N	1	1	670	310	1582	~3000±100	-

Table I :

 I Composition, deposition rate and thicknesses of the studied films

Table II Table II :

 IIII Stiffnesses C ij and indentation moduli M <hkl>

Table III Present study Virtual 1 Virtual 2 Virtual 3 Virtual 4

	Tip rounding (R=550 nm)	YES	YES	YES	NO (sharp) NO (sharp)
	Film thickness (3 µm)	YES	NO (bulk)	YES	YES	NO (bulk)
	Deformable indenter	YES	YES	NO (rigid)	YES	NO (rigid)
	) ( θ ℑ	(loading)	1.7	1.7	2.1	1.9	2.0
			(reference)	0%	+24%	+12%	+18%

Table III :

 III Tip rounding, deformable indenter and substrate effects on the identifiability index

).

In these two expressions the σ 0 /H B0 ratio is the argument of the f function. In this study due to the small values of this argument (σ 0 /H B0 < 10%) the corrections on the indentation modulus and the hardness are appreciably negligible (<2%).

Appendix A3

To identify the optimal set of parameters

⎣ is the maximum indentation load which makes ω dimensionless.

- θ and + θ are the minimum and maximum values of the parameters.

The constrained minimization process combines a global (genetic algorithm) [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF][64] and a local (Levenberg-Marquard algorithm) [START_REF] Levenberg | [END_REF][66] algorithms implemented in a Matlab toolbox MIC2M (mic2m.univ-fcomte.fr). An exploration of the field

) is performed using the global genetic algorithm to find different sets of initial input parameters 0 θ required in the local algorithm. Hence, the numerical uniqueness is considered as guaranteed if all of these 0 θ parameters lead to the same solution θ using the local algorithm. The determination of the uncertainty θ Δ of the