J Briottet

L Pierro

Labarre

Nicolas Lalanne

Jean-Claude Krapez
email: krapez@onera.fr

Christophe Le Niliot

Xavier Briottet

Jean Pierro

Luc Labarre

Development and validation of a numerical tool for simulating the surface temperature field and the infrared radiance rendering in an urban scene

Keywords: infrared thermography, building, urban, remote sensing, thermal bridge

come

Development and validation of a numerical tool for simulating the surface temperature field and the infrared radiance rendering in an urban scene

Introduction

Human activity in urban area strongly influences the regional and local environment.

For instance a practical consequence in terms of urban climatology is the development of urban heat islands. Urban climate is the result of heat and water flow processes between the urban surfaces and the environment. Infrared remote sensing is a powerful tool for measuring the thermal influence of these interactions on the urban surfaces at different spatial scales. As a matter of fact the infrared images acquired by airborne sensor inform about the thermal balance of urban components (buildings roofs and walls, roads, sidewalks,…) and, by merging all these effects at a larger scale, about the urban heat island [START_REF] Lagouarde | Experimental characterization and modeling of the nighttime directional anisotropy of thermal infrared measurements over au urban area: case of study of Toulouse (France)[END_REF]. Another application of infrared thermography is, at a higher resolution, the remote search for thermal disruptions on buildings envelope. This includes local flaws in the insulation layer, water ingress and thermal bridges [START_REF] Asdrubali | A quantitative methodology to evaluate thermal bridges in buildings[END_REF][3] [START_REF] Martı Ń Ocaña | Thermographic survey of two rural buildings in Spain[END_REF].

The radiance images are however the result of complex phenomena both on the thermal level (i.e. how is a particular temperature field produced?) and on the radiation level (i.e. what is the link between previous temperature field and the recorded radiance image?); therefore a theoretical model is often necessary for interpreting these images.

Later on this model is also necessary in the inversion process whose aim it to retrieve from the infrared images, first the true surface temperature field [START_REF] Oltra-Carrió | Analysis of the Performance of the TES Algorithm Over Urban Areas[END_REF] [START_REF] Datcu | Improvement of building wall surface temperature measurements by infrared thermography[END_REF], then some thermal parameters of the buildings walls/roofs [START_REF] Chaffar | Thermal characterization of homogeneous walls using inverse method[END_REF] [8] [START_REF] Albatici | Infrared thermovision technique for the assessment of thermal transmittance value of opaque building elements on site[END_REF] [START_REF] Monchau | Diagnosis of insulated building walls using passive infrared thermography and numerical simulations[END_REF], and the local heat losses [START_REF] Hénon | High resolution thermo-radiative modeling of an urban fragment in Marseilles city center during the UBL-ESCOMPTE campaign[END_REF] if they can ever be measured [START_REF] Vavilov | How accurate is the IR thermographic evaluation of heat losses from buildings?[END_REF]. This model has to take into account the energy budgets at the different surfaces (roofs, roads, walls, etc.) including the radiative interactions in a spectrum extending from visible to far infrared. It should then solve the dynamic heat transfer problem in the solid materials of the urban scene at a scale compatible with the sensor resolution, which can be of the order of a tenth of a meter when the purpose is, for example, to detect insulation defects or to characterize the status of thermal bridges. Heat transfer through a wall may be multidimensional.

Generally speaking, 2D transfer is observed at each corner of the building envelope (Lshape), and at thermal bridges like the floor/wall junctions (T-shape). 3D transfer is observed for example at a corner of a floor/wall junction. Finally, 1D transfer is observed over the building envelope far from previous thermal singularities. The computed surface temperature field has finally to be combined with the emissivity field in a radiative transfer module for getting the at-sensor infrared radiance field in the part of the [3-15µm] spectrum corresponding to the IR sensor bandwidth.

The energy exchanges between cities and the atmosphere can be simulated at a scale between the mesoscale (whole city) and the building scale; such models are for example TEB [START_REF] Masson | A physically-based scheme for the urban energy balance in atmospheric models[END_REF] and SM2-U [START_REF] Soux | A model to calculate what a remote sensor 'sees' of an urban surface[END_REF]. A second category of models includes the thermo-radiative models establishing the energy budget at a scale significantly lower than the building size (typically one meter or less). These models estimate all terms of the energy budget, in order to compute the temperature. A list of such models is described in table 1. However they only compute surface temperature or the total radiance field in the direction of a virtual sensor obtained by applying the Stephan-Boltzmann law.

The thermo-radiative model SOLENE [START_REF] Miguet | A daylight simulation tool for urban and architectural spaces -application to transmitted direct and diffuse light through glazing[END_REF] solves the heat conduction problem in walls with a nodal method based on a preliminary evaluation of the radiative heat fluxes on the boundaries. A wall is modeled with one or two layers and discretized with five nodes at most. The consequence of this simplification is that the model does not correctly capture the behavior of complex walls, of heavy inertial materials, and of the ground [START_REF] Hénon | High resolution thermo-radiative modeling of an urban fragment in Marseilles city center during the UBL-ESCOMPTE campaign[END_REF]. For this reason we retained SOLENE only for computing the solar heating of the urban surfaces during one or diurnal cycles for alimenting our thermal software.

The temperature field can be obtained by analytical or numerical methods. The 1D thermal quadrupole is an analytical method particularly well suited for multilayer walls [START_REF] Maillet | Thermal quadrupoles: solving the heat equation through integral transforms[END_REF]. However its extension to multidimensional problems is restricted to simple geometries. In case of thermal bridge geometries, numerical methods are more appropriate such as the Boundary Element Method [START_REF] Tadeu | Simulation of dynamic linear thermal bridges using a boundary element method model in the frequency domain[END_REF], the Finite Difference Method as used in EPS-r [START_REF] Clarke | Energy simulation in building design[END_REF], the Finite Volume Method [START_REF] Déqué | Effect of 2D modelling of thermal bridges on the energy performance of buildings Numerical application on the Matisse apartment[END_REF] [START_REF] Martin | Problems in the calculation of thermal bridges in dynamic conditions[END_REF] or the Finite Element Method [START_REF] Depecker | Design of buildings shape and energetic consumption[END_REF]. Numerical methods would however be excessively time consuming when applied on an urban scene fragment. We thus propose a hybrid method combining the finite element method and the quadrupole method. The former is applied on walls and roofs in the vicinity of the thermal bridges and the latter is applied elsewhere on the buildings envelope and on the ground.

Once the temperature field is known, the radiance images at the sensor level can be calculated. A radiative transfer code is required to simulate the multiple reflections between the surfaces of the scene and the atmosphere contributions. Atmospheric radiative transfer codes like Modtran [START_REF] Berk | MODTRAN4 radiative transfer modeling for atmospheric correction[END_REF] and MATISSE [START_REF] Labarre | An overview of MATISSE-v2.0[END_REF] are based on the radiative transfer equation and can calculate the spectral contributions of the atmosphere along the line of sight (transmission, self-emission, molecular/aerosol scattering…). They are however limited to a simplified geometry (plane or spherical ground). 3D radiative transfer codes were thus developed in order to take into account radiative interactions between 3D elements of a complex scene like urban scenes. In the infrared domain, DART [START_REF] Gastellu-Etchegorry | 3D modeling of satellite spectral images, radiation budget and energy budget of urban landscapes[END_REF] and TITAN [START_REF] Fontanilles | Thermal infrared radiance simulation with aggregation modeling TITAN an infrared radiative transfer model for heterogeneous 3D surface -application over urban areas[END_REF] are such codes. DART is based on a 3D discretization of the urban scene and of the atmosphere in cubic voxels and the associated model DARTEB which allows to obtain the energy budget simulation and the 3D distribution of temperature uses the same thermal model as TEB. TITAN is an accurate urban radiative transfer tool for performing fine phenomenological analysis of remote sensing measurements [START_REF] Fontanilles | Thermal infrared radiance simulation with aggregation modeling TITAN an infrared radiative transfer model for heterogeneous 3D surface -application over urban areas[END_REF] but it is limited to the MWIR/LWIR and supposes that the surface temperature is known. Thereafter, the new radiative transfer code MOHICANS (MOdélisation Hyperspectrale d'Images en entrée Capteur pour l'ANalyse et l'inversion du Signal -Model of at-sensor hyperspectral images for analysis and inversion of the signal) is proposed which results of the combination of AMARTIS [START_REF] Thomas | AMARTISv2: 3D radiative transfer code in the [0.4; 2.5 µm] spectral domain dedicated to urban areas[END_REF] for the visible to 2.5 µm spectral band and of TITAN for the 2.5 to 14 µm spectral band. The interest of this code is to be able to simulate the signal of a hyperspectral sensor in the whole [0.4 -14 µm] spectral band.

As a consequence of what was presented before, we developed the thermal code SUSHI (Simulation in Urban Scene of Heat dIffusion) and coupled it with Sketchup, SOLENE and MOHICANS for getting a software tool able to simulate the whole thermal and radiative chain of processes at the origin of infrared radiance images of an urban scene in the thermal IR spectral domain [3-15 µm] and at a decimeter resolution.

The numerical tool we developed is aimed at simulating the dynamic thermal behaviour of walls presenting areas where heat flow is either 1D or 2D. For getting a simulating IR image close to the one measured in a particular condition it is necessary to provide to the thermal software input data that finely describe the environment heat forcing (solar flux, atmospheric IR flux, internal and external convection). These external/internal conditions should be recorded for a sufficiently long period for enabling the dynamic thermal analysis to capture all thermal inertia effects.

In the first section we will describe the models, in particular SUSHI and MOHICANS. The first two points are common to all dimensional studies (1D or 2D) namely the hypothesis regarding the time evolution and the boundary conditions. In particular we will describe the modules used from SOLENE for computing the solar and infrared flux absorbed by each surface of the urban scene. Then we present the 1D thermal scheme based on the quadrupole method and the 2D scheme based on a FE method. The characteristics of the spectral model MOHICANS for the infrared rendering are finally recalled.

The second section is devoted to computation results. We will start by analyzing the validity of the hypothesis about the periodic character of the environment heat forcing. A short comparison of 1D SUSHI results with those obtained with SOLENE thermal model will then be presented. A detained sensitivity analysis regarding a typical 1D wall will allow sorting the thermal and radiative parameters according to their impact on the external surface temperature.

The third section will present some results obtained during the BATIR 2013 measurement campaign performed with two ground-based cameras and an airborne camera on two close buildings during the heating period for the validation of our numerical tool. The 1D temperature results will be compared to thermocouple readings and the 2D infrared radiance field will be compared with LWIR images recorded over a whole façade presenting different types of thermal bridges.

Models description

First we will describe the main aspects of the thermal code SUSHI and its coupling with SOLENE and MOHICANS. A series of programs are actually used as described in the chart in figure 1: • SketchUp is used to build-up the 3D-scene and subdivide each wall in areas where heat flow is either 1D or 2D;

• SOLENE is used to mesh each area, to compute the view factors and then the incident solar radiation after multiple reflections and shadow tracing at each time step;

• SUSHI (Simulation in Urban Scene of Heat dIffusion) is used to model the infrared radiation energy budget, the heat convection flux and then to solve the 1D and 2D heat transfer equation inside the walls and the soil elements for getting the surface temperature;

• MOHICANS is finally used to get the spectral infrared radiance incident to the sensor from the temperature field computed by SUSHI and from the atmosphere radiative parameters computed by MATISSE [START_REF] Labarre | An overview of MATISSE-v2.0[END_REF].

Thermal model: SUSHI

Time evolution

Temperature and fluxes are assumed to be periodic. This allows applying the Fourier transform which, thanks to the Fast Fourier Transform tool, leads to semi-analytical solutions and yields the results in a very short time. However, for simulating the façade temperatures at best, we have to consider real meteorological data for the input, i.e. solar radiation, air temperature, etc… These data can be seen as random variables; the 24h periodic character is merely an approximation. When compared to the transient approach, for which an initial condition has to be specified (i.e. the initial temperature distribution in the walls and in the soil), the periodic approach relies on an initial condition which is actually the result of the repetition of an infinite number of periods of the input fluxes. For approaching the true aperiodic case, one should either reject the initial time far in the past for the transient approach, or take into account a long period for the periodic approach. In both cases it means recording meteorological data for a long duration before the time at which temperature has to be computed. We will call it the "pseudo-period". The question was then how to choose the pseudo-period aimed at representing the aperiodic phenomena. Of course it should be 24h at least. If it is too short, let us say only 24h, the computed temperature would be in error due to fact that the atmospheric conditions were actually different any previous day. One should thus extend the duration of the pseudo-period. In this way, the first mismatch between the modelled atmospheric conditions (i.e. through the periodicity assumption) and the real atmospheric conditions would be rejected further in the past. The error on the computed temperature would then be lowered at the expense of a longer computation time. The optimal choice for the pseudo-period duration will be discussed later after completing the presentation of the thermal model.

Boundary conditions (BC)

Each element of the mesh can exchange heat fluxes with the environment. In the general case, these fluxes are of convective and radiative nature. The radiative fluxes are subdivided into solar and infrared components. Both solar and infrared radiation can be partly absorbed and reflected by the surfaces. The solar radiation may be split into a direct and a diffuse component in order to give hourly shadow distribution. SOLENE is used to compute these solar components. In this software, the direct component follows the Perrin de Brichambaut formula, and the atmospheric diffuse component follows the "all weather model" of Perez et al. [START_REF] Perez | All weather model for sky luminance distribution -preliminary configuration and validation[END_REF]. This provides spectral-integrated fluxes depending on place and hour, and on sky clearness and brightness for the diffuse component. The energy budget through multiple reflections is computed by the progressive refinement radiosity method [START_REF] Cohen | A progressive refinement approach to fast radiosity image generation[END_REF] which first requires the calculation of the view factors between all surfaces of the scene including the "sky facets" (all urban surfaces are assumed to be Lambertian for the heat radiation calculation step). SOLENE however doesn't apply the progressive refinement radiosity method for evaluating the infrared radiative net fluxes, it rather considers a blackbody assumption. For avoiding this approximation we chose to perform a more rigorous infrared radiative exchange calculation by implementing the same progressive refinement radiosity method as for the solar radiation.

Let us now consider the convective flux on the external building surfaces. It may be modelled through a convective coefficient which depends on wind velocity and orientation, on air temperature, on local and global geometry, etc. A 3D fluid dynamics solver would provide a precise modelling of this coefficient [30] [31]. However, for keeping a reasonable computational time, we preferred to use one among the numerous empirical formulas giving the convective coefficient vs. wind velocity in the building thermal energy context [START_REF] Palyvos | A survey of wind convection coefficient correlations for building envelope energy systems modeling[END_REF]. The coefficient is thus time dependant:

( ) t h ex .
The thermal balance at any facet of the outdoor building surface can be expressed as follows (for an inward oriented z axis):

( ) ( ) ( ) env ex net IR net VNIR air ex ex ex ex T T t t t T T t h z T , , ) ( ) ( ϕ ϕ λ + + - - = ∂ ∂ - (1)
It can be considered as a second order BC where the specified heat flux corresponds to the right-hand side of the equation. T . This dependence is nonlinear which requires solving the problem through an iterative process. The time dependence of ( ) t h ex would also require iterations even though an alternative exists which avoids performing iterations but requires more involved algebra [START_REF] Krapez | Modélisation analytique de la variabilité spatiale et temporelle de la température de sols nus naturels[END_REF].

For speeding up the iterative process, the second order BC can be replaced by a third order BC where the linear heat exchange is given an arbitrary constant value, let us call it the fictitious transfer coefficient fic h . The BC now expresses as:

( ) ( ) ( ) t T t T h z T sa ex fic ex - - = ∂ ∂ -λ (2)
where ( ) t T sa corresponds to a modified sol-air temperature defined by:

( ) ( ) ( ) ( ) ( ) [ ] env ex net IR net VNIR air ex ex ex fic ex sa T T t t t T T t h h t T t T , , ) ( ) ( 1 ϕ ϕ + + - - + = (3)
In other words, the thermal excitation associated with the 3 rd order BC is:

( ) ( ) ( ) ( ) env ex net IR net VNIR air ex ex ex ex T T t t t T t h T t h , , ) ( ϕ ϕ + + + ∆ - (4)
where

( ) ( ) fic ex ex h t h t h - = ∆ .
A first iterative process identified by ( 1

1D configuration

As explained before, the 1D transfer model is used for the wall facets that are far from thermal singularities like façade corners, wall/floor junctions, … . Its semi-infinite version is used for soil facets. The walls may be composed of many layers. It is assumed that each layer is homogeneous, isotropic with constant thermal properties.

The thermal problem is solved by applying the quadrupole technique after performing a Fourier transform to the heat equation and to the boundary conditions [START_REF] Maillet | Thermal quadrupoles: solving the heat equation through integral transforms[END_REF]. The following backward transfer relation relates the vectors of temperature and conduction heat flux density at the external side and the internal side of the wall:

        =         in in ex ex T T ϕ ϕ M ( 5
)
The tilde represents the Fourier transform, and M is the global quadrupole matrix obtained by multiplying the N individual quadrupole matrices, N being the number of layers of the wall [START_REF] Maillet | Thermal quadrupoles: solving the heat equation through integral transforms[END_REF]:

( ) ( ) ( ) ( ) ∏ ∏ = =           =         =         = N i i i i i i i N i i i i i j j j b j j b j D C B A D C B A 1 1 cosh sinh sinh 1 cosh ω ξ ω ξ ω ω ξ ω ω ξ M (6)
For each layer i,

i i i a l = ξ
is the square root of the diffusion time through the layer, i a is diffusivity, i l is thickness and i b is effusivity. ω is the pulsation and j the imaginary variable 1 -.

From eq. ( 5) one can express the external temperature as a function of external and internal net flux densities at iteration k :

in k ex k ex Q Z Q X T , , + = (7)
with the external transfer function ( )

B h A X in ′ + ′ Λ = -1 (8)
the internal transfer function:

              - Λ = ∑ = - ω ξ j Z N i 1 i 1 exp (9)
the internal net flux density:

air in in in T h Q = (10) D h C B h h A h in in fic fic ′ + ′ + ′ + ′ = Λ (11)
and the external net flux density is defined iteratively through:

net k IR net VNIR air ex ex k ex ex k ex T h TF T h TF Q 1 , 1 , , ) ( ) ( ~-- + + + ∆ - = ϕ ϕ (12)
A', B', C' and D' are the reduced quadrupole coefficients, i.e. as obtained after factorizing each quadrupole matrix by the exponential of positive argument

( ) ω ξ j i exp
. This factorization and the subsequent simplifications in the transfer functions X and Z avoid round off errors otherwise occurring when the magnitude of the argument of the hyperbolic functions in eq. ( 6) reaches high levels ( [START_REF] Pailhes | Thermal quadrupole method with internal heat sources[END_REF], [START_REF] Krapez | Three methods based on the thermal quadrupole approach for an improved computation of temperature and flux in multilayered materials with internal heat sources[END_REF]).

The temperature is then obtained after an inverse Fourier transform of k ex T , ~ by applying the IFFTW algorithm [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. The term

2D configuration

When the one-dimensional assumption for heat conduction is not valid, like in the vicinity of thermal bridges, we substitute the 1D transfer functions X and Z by 2D transfer functions that are pre-calculated by a finite element method. For this purpose we use the software CAST3M. The objective is to get periodic transfer functions corresponding to the response of a thermal bridge when an elementary flux is applied over each one of the characteristic surfaces of the thermal bridge. Thus, when considering the symmetric L-shape in figure 3-left, a uniform flux is applied on surface 1 S , then on surface 3 S (due to symmetry, the response is the same when the perturbation is applied on 1 S or on 2 S ). Similarly, for the symmetric T-shape in figure 3-right a uniform flux is first applied on 1 S , then on 2 S . In all cases, a unitary step heating is applied for a relatively long duration, typically several days. The transient response of the system is then subtracted to itself after a delay of one time step. This yields the system response to an elementary pulse. N replica of this response delayed by multiples of 24h are then added together for getting the periodic response to an elementary pulse (N should be sufficiently long for approaching the true periodic response). One is essentially interested in the external temperature induced by a uniform heat flux applied through periodic pulses on each of the surfaces

∑ = = 3 1 ) ( ). , ( ) , ( ~l k l k p l k ex Q x G x T ω ω ω (13)
where ) ( ~k l Q ω is the Fourier transform of the net heat flux density averaged over surface l S and x is the position on the external surface. As before, for returning back to the temporal space an inverse fast Fourier transform is then performed with the IFFTW algorithm [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. Present method involving integrated Green functions could be compared with the weighting factors method in [START_REF] Depecker | Design of buildings shape and energetic consumption[END_REF] when conversely, the conduction flux is expressed as function of surface temperatures.

The lateral size of the 2D model is chosen sufficiently high (typically 1m from the thermal bridge centre) so that the response at its edges is very close to the 1D response. Then we retain only the 2D nodes where the response is significantly different from the 1D response.

Spectral model: MOHICANS

The software MOHICANS is a numerical tool developed at ONERA which aims at evaluating all radiative components contributing to the optical signal received by a multispectral or a hyperspectral sensor over a 3D scene. It is actually a fusion of two previously developed codes: AMARTIS [START_REF] Thomas | AMARTISv2: 3D radiative transfer code in the [0.4; 2.5 µm] spectral domain dedicated to urban areas[END_REF] for the visible to near infrared range and TITAN [START_REF] Fontanilles | Thermal infrared radiance simulation with aggregation modeling TITAN an infrared radiative transfer model for heterogeneous 3D surface -application over urban areas[END_REF] for the midwave to longwave infrared range. Both of them rely on atmospheric radiative contributions which are computed by the ONERA code MATISSE [START_REF] Labarre | An overview of MATISSE-v2.0[END_REF]. The spectral domain finally covers the 0.4 to 14 µm range.

MOHICANS is organized in four calculation steps: (i) computation of the sky irradiance hemisphere and the atmosphere transmission terms by calls to MATISSE, (ii) computation of geometrical terms corresponding to sun visibility, facet to facet visibility and facet to pixel visibility, (iii) computation of total irradiance, corresponding to the sum of solar direct irradiance, the atmosphere irradiance and the environment irradiance, (iv) computation of the radiance of the radiative flux reaching each pixel of the imaging sensor. A schematic of the radiative exchanges between two facets of a wall-soil angle is presented in figure 4.

Computation results

Temporal configuration

The first analysis concerns the validation of the time-periodic model. For simulating the temperature field over a one day time interval [t 0 -24,t 0 ] we assume that the radiative and convective fluxes are periodic with a pseudo-period D and that their variations correspond to those observed experimentally during the time interval [t 0 -D,t 0 ] with D=kx24h; k≥1. As explained earlier, due to weather changes, the experimental fluxes unavoidably show random variations. Despite an apparent 24h periodicity they are not perfectly periodic. One way to circumvent this problem consists in taking a longer interval D, i.e. several days, and considering it as the new baseline for periodicity. This amounts to rejecting further to the past the unknown initial condition. As a matter of fact, the longer the pseudo-period D, the lower the error on temperature over the interval [t 0 -24,t 0 ] due to the wrong assumption about periodicity. The minimum number of days k needed for getting an error lower than a specified value was determined by considering both an elementary part of a typical wall (see properties in table 2) with an albedo of 0.5 and an emissivity of 0.947 and a semi-infinite ground with an effusivity of 1000 and 2000 J.s -1/2 .m -2 .K -1 with same optical parameters. Both were submitted to the real atmospheric fluxes observed during the BATIR campaign for a horizontal surface (see last section). A long pseudo-period D of 12 days was considered for getting the reference solution. Temperature was computed for progressively shorter pseudo-periods D and compared to this reference. The maximum difference and the RMS error value over the day of interest (i.e. the last one) are plotted in figure 5 for the wall and for the ground. The relative difference in figure 6 is obtained by dividing previous errors by the day-night temperature difference observed during the last day, which is 13.3 K for the wall, 18.2 K for the low effusivity soil and 13.2 K for the high effusivity soil. These simulations show that a pseudo-period of one day is too short both for the wall and the soil temperature simulations: the maximum error is between 0.3 and 0.6 K. A pseudoperiod of two days is long enough for the wall: both maximum error and RMS error are less than 0.01 K. On the other side a pseudo-period of five days is necessary for both soils when setting the maximum error criterion to 0.1 K. A soil with higher thermal inertia needs a longer pseudo-period for getting the same error. For the remaining we will perform the thermal simulations with a pseudo-period of 48h. The error is guaranteed to be low for the walls, but one has to be aware that it may be up to about 0.8 K (i.e. 6% of day-night amplitude) for high effusivity soils, depending on the environmental fluxes irregularity.

Software validation and comparison with SOLENE

Analytical validations of SUSHI model performed on basic configurations like a semiinfinite layer submitted to a sinusoidal boundary flux with or without convective losses showed a perfect match. A comparison with SOLENE results showed the limitations of this nodal model for simulating soils and very thick materials. Comparisons with SOLENE were also performed for standard walls: in the case of a single two layer slab both models give similar results. However, SUSHI permits to model walls made of more than two materials as opposed to SOLENE.

Sensitivity analysis

The sensitivity of the front temperature of a three layers wall was analysed with respect to the thermal, convective and radiative parameters. A series of 14 parameters was considered: the effusivity and the square root of the diffusion time of each layer (see table 2), the albedo and the emissivity of the external surface, the external and internal convective coefficients, the external and internal air temperature, and the VNIR and IR incident flux density. Their nominal values correspond to those observed for the façade of the building LC4 during the measurement campaign BATIR (see fig. 8).

We defined two groups of parameters, according to their time dependency. The external convective coefficient, the external and internal air temperature and the incident radiative flux densities all vary with time. For static parameters β we introduced, as usual, a specific relative variation with respect to their reference value . In both cases the reduced sensibility to a parameter i β is approximated by the finite difference:

( ) K T T T S i i i i i i β β β β β - ≈ ∂ ∂ ≡ ∆ + (13)
In the case of a time-dependent parameter ( ) t β the application of a systematic relative variation has the disadvantage of minimizing the sensitivity to ( ) t β when it comes close to 0. We thus introduced another type of analysis, based on adding random variations to the nominal time evolution ( )

t ref β according to: ( ) ( ) ( ) 1 2 1 , 0 - + = u K t t ref ref β β β (14)
where ref β is the time average of ( )

t ref β , 1 << K
and u 0,1 is the uniform distribution on the support [0,1]. A series of stochastic simulations has to be performed for a given K value and the temperature scatter is then compared with the nominal temperature evolution. The reduced sensitivity to ( ) t β is then given by the ratio of the standard deviation of the simulated temperature to the standard deviation of the uniform distribution of the considered parameter, i.e.:

( ) ( ) K t t S T i 3 σ = ( 14
)
We plotted the reduced sensitivity to the thermal parameters of the three layers in figure 6. The thermal behavior of any layer can be described with two parameters taken from the list (

) i i i i b c R ξ , , ,
where i R is thermal resistance

i i l λ and i c is thermal capacity i i i l C
ρ . We started with the pair ( )

i i b ξ
, and noticed that for layers 2 and 3 these parameters are highly correlated. Choosing ( )

i i c R ,
for these layers shows that temperature is mostly sensitive to i R and nearly not to i c (the sensitivity to the capacity of the third layer is not plotted because it is less than 10 -3 K in absolute value). In figure 6 a reduced sensitivity of 1 K to a given parameter means that a deviation of 10% of this parameter yields a front temperature deviation of 0.1 K. We can observe that the sensitivities to effusivity and square root of diffusion time of the first layer are decorrelated: they are opposite in sign and present a delay. Maximum absolute values of sensitivity are with respect to the front layer effusivity; this sensitivity is negative at day and positive at night; it reaches about 0.3 K per 10% effusivity variation. The sensitivity to the middle layer resistance 2 R is rather stable and it is a bit lower (in magnitude) than the sensitivity extremes with respect to 1 b and 1 ξ . The sensitivity to the third layer resistance 3 R is still one order of magnitude lower.

We can see in figure 7-a the sensitivity of the front temperature to the convective coefficients and to the external and indoor air temperature. For the timevarying parameters, we present both the stochastic sensitivity as obtained through eq. ( 14) (curves with "stoc" label) and the sensitivity obtained after systematic flux alteration (curves with "syst" label). Regarding the external parameters, the air temperature is most critical to the results when the convective coefficient is high. For the two days presented here, the sensitivity to the air temperature reaches about 240 K for the systematic variation analysis, respectively 100 K for the stochastic analysis. It means that an error of 0.5% on the air temperature (in Kelvin), i.e. around 1.4 K, yields an uncertainty of around 1.2 K, resp. 0.5 K, on the front surface temperature. The sensitivity to the internal air temperature is from 10 to 100 times lower. Similarly, the sensitivity is lower to the internal convection coefficient than the external convection coefficient.

The optical parameters and the incident flux densities define the radiative front excitation. In figure 7 As a conclusion on the sensitivity analysis we can state that for the present wall example the front temperature is mostly sensitive to the thermal and geometrical parameters of the first layer and to the thermal resistance of second layer, to the outdoor convective parameters (convective coefficient and air temperature), to the IR incoming flux and to emissivity. Therefore these parameters have to be measured with a high accuracy for getting the absolute temperature of the wall with a small error: for example the uncertainty has to be lower than about 2% for the IR incident flux, or 0.6 K for air temperature, to get an error lower than 0.5 K. Incidentally, this analysis shows that the third layer and the indoor heat transfer coefficient would be extremely difficult to evaluate from the external surface temperature through a data inversion.

Experimental campaign

Description

The experimental measurement campaign BATIR took place in winter 2013 on the military airbase 701 at Salon-de-Provence, France. The thermal behaviour of two buildings was studied for a period of 13 days, from January 19 th to 31 st 2013 (see figure 8). Two meteorological stations were used: one station ( of the temperature and of the heat flux at particular locations of the wall. This wall is made of three layers, consisting of a reinforced concrete layer, an isolating PS layer and a BA13 plaster board. The respective layer thicknesses and thermal properties are given in table 2. The effusivity of the 1 st and 3 rd layers was measured with a Thermicar sensor by implementing the asymmetric hot plane method [START_REF] Krapez | Caractérisation thermique non intrusive au moyen d'un plan chaud et d'un fil chaud asymétriques[END_REF]. Their conductivity and other missing thermal data were taken from literature.

Airborne thermography was also performed with a FLIR A325 microbolometer camera onboard of a motor glider. The camera was flown at about 40 m.s -1 speed and about 400 m altitude. It was aiming at nadir. During the first thermography period, two flights were performed short after sunrise, at midday, and during the second period, three flights were performed short after sunrise, at midday and shortly before sunset.

Experimental validation of SUSHI

The first experimental validation consisted in comparing our 1D thermal model results with the thermocouple measurements on the north facade of LC4 building. It was located far from the corners and from any thermal bridge so that the 1D hypothesis for heat flow in the wall was justified. The outside convective coefficient was calculated from wind speed by the empirical formula of Croiset: [START_REF] Soux | A model to calculate what a remote sensor 'sees' of an urban surface[END_REF] with wind speed measured between the two buildings. Wind speed data were available from two meteorological stations, one in an open area, and the other between the buildings. The first data series was available for whole month and the second one only for the last day of interest (started at 10:00). One day of data was thus missing for allowing an analysis with a pseudo-period of 48h as described before. Data from the first station were thus used for the missing day by applying a correction factor that was established by comparing the wind data recorded during the second, i.e. common, day. This extrapolation is justified by the fact that the wind direction didn't change during this couple of days. The outside air temperature was given by the remote station. The inside convective coefficient was arbitrarily fixed to 10 W.m -2 .K -1 , with air temperature measured at an approximate distance of 5 cm from the surface. Let us mention here that only the net radiative solar flux was obtained from the temperature 3D field estimated by the SOLENE thermal code; the nested iterations aimed at updating the convective heat flux term and the radiative infrared flux term were thus performed by SUSHI (see figure 2). Temperature measured with a given thermocouple is compared with the simulated temperature of the corresponding facet. These thermograms are shown in figure 9. The uncertainty (at a confidence level of 95%) of the thermocouple is 0.14°C, and the day-average uncertainty of the simulated temperature is estimated at 2.2°C, with maximal of 2.92°C, as illustrated by the +2σ and -2σ curves in figure 9. As the sensitivity study has shown a small systematic error on the external air temperature and on incident infrared radiation is dramatic for the precision of the simulated temperature. This is why the uncertainty of the simulated temperature was conservatively estimated by considering a systematic error (and not a "stochastic" error) of 5% for all constant parameters, for external convective coefficient, for the incident solar and infrared radiation, of 0.5°C for external air temperature, and of 0.1°C to internal air temperature.

( ) 605 . 0 1 . 4 1 . 3 V t h ex + =
Notwithstanding the relatively important value for the uncertainty of the simulated temperature, the RMS difference between the simulation and the measurements is only 0.33°C (with a maximum difference of 0.78°C). The observed deviation between 03:00 and 09:00 may be explained by the presence of clouds at the sunrise. In the end, the thermal simulator provides satisfactory results; a significant effort has however to be undertaken for precisely evaluating the most influent environment parameters. The infrared cameras are then used for a spatial validation.

Experimental validation of the full chain of simulation

The façade of the building called "LC4" is used for the validation of the full chain of programs. The figure 10 The image presented in fig. 10-left was recorded at 05:00 PM on January 29th.

For the computation, a pseudo-period of two days was used to take into account the thermal dynamic of the scene. The scene was triangulated with about 10000 facets. Some areas with known thermal bridges were meshed more densely, namely the Tshape bridges corresponding to floor/façade junctions (they are clearly seen on the IR image in Figure 10-left). For convenience, the ground, the roof, the windows and the door contained few facets and we discarded the 45° T-shape bridges corresponding stair/façade junction. At present, semi-transparency isn't modelled; this is why temperature and radiance of the windows aren't realistic. For the example presented here we didn't model the L-shape bridge corresponding to the building corners. The image of the relative difference shows that the radiances of 1D-wall parts are correctly estimated with a relative difference less than 2.5%. One should nevertheless have in mind that relative RMS difference to about 12%. Again, only with adapted meteorological sensor can we hope getting simulation results close to the real infrared sensor data. The effect of the thermal bridges at floor junctions on the radiance map is relatively well retrieved by the simulation despite the size of the facets used here.

Differences observed for the ground and the roof can be explained by the high "roughness" of these surfaces (they were considered flat in the 3D model). The difference reaches 9% at the soffit (the underside of the roof): the radiance of this part is higher than revealed by the simulation because the internal temperature can be higher than the measured temperature and the material is unknown. From this part some warm air could flow from the attic against the façade, which also could explain the difference just below the roof.

Conclusion

In this work a numerical tool was developed for the simulation of the temperature field and the infrared rendering in an urban scene. The radiation energy budget at each facet of the modelled environment takes into account the multiples reflections of the incident solar flux (by SOLENE) and the IR radiative coupling between facets (by the new code SUSHI). The 1D response factors are calculated interactively by the quadrupole method whereas the 2D response factors are extracted form a database that was built from finite element calculation results. The model is based on a periodic approach. With real meteorological data, a pseudo-period of two days is sufficient for reducing the error induced but the periodicity hypothesis to satisfactory levels: less than 0.03 K for a typical wall and less than 0.8 K for the ground. The sensitivity study indicates the importance of the thermal properties of the first two layers of a wall, and the influence of measurement errors about external air temperature, emissivity and infrared incident flux. Validations of this tool were performed by comparing the output with thermocouples and infrared thermography measurements obtained during the experimental campaign BATIR. The final results show relatively little difference between measured and simulated radiance for the wall and the considered thermal bridges. This should not hide the fact that when combining all input uncertainties, the simulated radiance uncertainty may be significant, 10% and more, highlighting the requirement for a precise monitoring of the atmospheric thermal and radiative contribution.

Higher differences are observed for the roof (because of the geometry of the tiles) and the ground (because of its roughness) and at some locations for which boundary conditions were not well known (soffit and crawl space level). Further studies will be performed for analysing the airborne data.

Similar geometrical computations are made by SOLENE and by MOHICANS for solar interactions in the scene. The next step will thus be to unify there computations and to use a unique module for the form factor matrix evaluation. In the future, the adding of a semi-transparent material model in SUSHI will enable us to consider more general situations with glazing surfaces and vegetation.

k (W.m -1 .K -1 )
ρ.Cp (kJ.m -3 .K -1 ) b (J.s -1/2 .m -2 .K -1 ) ξ (s 1/2 ) concrete 0.17 1.8 [START_REF] Giaconia | On the reliability of ASHRAE conduction transfer function coefficients of walls[END_REF] 2112 1950 (1) 184.1 EPS 0.027 0.037 [START_REF] Wyart | Polystyrène expansé ou PSE[END_REF] 21 [START_REF] Wyart | Polystyrène expansé ou PSE[END_REF] 27.9 20.3

Plaster (BA13) 0.013 0.31 825 [START_REF] Kuznik | Optimization of a phase change material wallboard for building use[END_REF] 504 (1) 21.3 (1) Values given by on-site measurements

ϕ

  air ex T is the external air temperature, net VNIR ϕ is the net radiative flux in the visible and near infrared domain (i.e. the solar spectrum, is the net radiative flux in the thermal infrared domain ( µm 3 > λ ); it depends on the atmospheric IR radiation, on the temperature of the surfaces exchanging radiation with it, all these environmental contributions being symbolically expressed through env T , and on the temperature of the surface of interest itself, ex

  ) in figure2is implemented by updating only the pseudo-convective flux term step (2) in figure2) and the global iterative process is repeated until a second global convergence criterion is reached. The iterative process is split into two steps because updating the radiative flux terms is the most expensive. For the indoor BC, it is assumed that the environment is homogeneous and it is characterised by an indoor air temperature constant radiato-convective transfer coefficient in h .

  runs on the surfaces 1 S and 2 S for the L-shape and on the surface 1 S for the T-shape (the responses ) as integrated Green functions). The Fourier transform of these periodic responses 3 for future computations of the external temperature according to:

  -b we plotted the reduced sensitivity (in absolute values) to the radiative parameters and to the radiative incident fluxes. Without surprise, the sensitivity to albedo and to a systematic variation of the incident solar flux are opposite and correlated in module. Sensitivity to a systematic variation of the solar heat flux is naturally maximum around midday: it reaches about 0.1 K for a 10% flux variation. It doesn't drop to 0 during night hours simply because the flux increase occurring during sunlit hours has also a positive impact on night temperatures. The sensitivity to albedo is homothetic to the sensitivity to incident solar flux by a ratio of about 3.6, a value which can be explained by the actual value of albedo (i.e. 0.78). The sensitivity to the stochastic variations of the solar flux is lower: it is roughly 0.03 K for a 10% random flux variation. The fact that it is lower is mainly due to the inertial character of the wall: rapid variations of the solar flux (at the 15 min considered time step) are strongly attenuated with respect to a slow 24h flux variation as considered in the systematic variation analysis. The noticeable result is that the front temperature is mostly sensitive to the infrared atmospheric flux (from 0.8 to 2.8 K per 10% systematic variation). As a consequence, this flux must be measured with a great care for expecting a good precision on the simulated temperature. The sensitivity to the external surface emissivity is from 3 to 10 times lower than the sensitivity to IR flux, in accordance to the fact that emissivity drives both the absorbed flux and the self-emitted flux (nominal emissivity was 0.947).

  Campbell CR23X) was placed in an open area next to the runway and the other one (Campbell CR3000) was placed between the two buildings. Meteorological data were collected in order to estimate outdoor air temperature, wind velocity, global solar irradiance and IR irradiance. Air temperature inside the rooms was measured by probes (EL-USB-1-PRO) suspended close to back side of the façade wall of interest. The spectral reflectance of the participating surfaces was measured with an ASD spectroreflectometer in the [0.4-2.5 µm] spectral domain and a SOC400T spectroreflectometer in the [2.9-13.3 µm] spectral domain. The albedo was calculated for each surface from the first spectral reflectance by weighting it with the solar spectral radiation provided by the atmosphere radiative transfer model MATISSE. Similarly the average IR emissivity was calculated from the measured IR spectral reflectance by weighting it with the Planck's law corresponding to a surface at 278 K. The thermography analysis was realised during two separate periods of 24 hours. Two 512x640 IR cameras were used: a band II camera (FLIR SC6000) and a band III camera (FLIR SC7000). The latter camera has been provided with a narrow filter. The halfheight spectral bands are respectively [3-5 µm] and [8.4-8.7 µm]. The cameras were periodically aimed at an extended blackbody for drift correction and for establishing a relation between the output signal and the IR radiance in the respective bandwidths. For the first 24h period thermographic measurements were performed by aiming the south façade of building LC3 from the north façade of building LC4. For the second period we performed the opposite. During this second period, five thermocouples and a fluxmeter were fixed on the north façade of building LC4 for getting independent values

  presents three images of this façade: (a) is the infrared image measured by the camera SC7000 (in the [8.4-8.7 µm] spectral band), (b) is the simulated image obtained at the outlet of the simulation chain, (c) is the relative difference between the real and the simulated image and (d) is the vertical profile of radiance along the right side of the façade. The uncertainty (at a confidence level of 95%) of the measured radiance is about 0.11 W.m -2 .sr -1 . The main source of error for simulated radiance is the uncertainty of the transmitted-emitted radiance of the wall, which can be analytically computed. For an uncertainty of 1.1°C for the wall temperature and of 0.03 on the product emissivity-transmittance, the uncertainty of the simulated radiance is about 0.21 W.m -2 .sr -1 .

Figure 2 .Figure 3 .

 23 Figure 2. Flow diagram of the parameters/variables computed by SOLENE and by SUSHI

Figure 4 .

 4 Figure 4. Radiative components computed by MOHICANS and contributing to the radiance at sensor level.

Figure 5 .

 5 Figure 5. RMS error and maximum absolute error (Max) for surface temperature depending on the number of days considered for the pseudo-period in the case of a wall and a ground of effusivity 1000 SI or 2000 SI. The reference is obtained with a pseudoperiod of 12 days. (a) absolute error and (b) relative error to the min/max temperature contrast obtained for current day.

Figure 6 .

 6 Figure 6. Relative sensitivity of external temperature to the thermal properties of the three-layer wall ( i b , i ξ , i R and i c are the effusivity, the square root of diffusion time, the thermal resistance and the thermal capacity of layer i).

Figure 7 .

 7 Figure 7. Same with respect to the convective parameters (a) and to the radiative parameters (b).

Figure 8 .

 8 Figure 8. Analysed scene during the BATIR campaign (a) and airborne IR image (b).

Figure 9 .

 9 Figure 9. Comparison between the temperature simulated by SUSHI and the temperature measured by a thermocouple on the façade of LC4 building.

Figure 10 .

 10 Figure 10. Comparison between the radiance measured by thermography (a) and the radiance simulated by SOLENE/SUSHI/MOHICANS (b). Radiance contrast for the image (c) and radiance profile along the façade of LC4 building (d).

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Flow chart of the simulator

Figure 4 .

 4 Figure 4. Radiative components computed by MOHICANS and contributing to the radiance at sensor level.

Figure 5 .

 5 Figure 5. RMS error and maximum absolute error (Max) for surface temperature depending on the number of days considered for the pseudo-period in the case of a wall and a ground of effusivity 1000 SI or 2000 SI. The reference is obtained with a pseudoperiod of 12 days. (a) absolute error and (b) relative error to the min/max temperature contrast obtained for current day.

Figure 6 .Figure 7 .

 67 Figure 6. Relative sensitivity of external temperature to the thermal properties of the three-layer wall ( i b , i ξ , i R and i c are the effusivity, the square root of diffusion time, the thermal resistance and the thermal capacity of layer i).

Figure 8 .

 8 Figure 8. Analysed scene during the BATIR campaign (a) and airborne IR image (b).

Figure 9 .

 9 Figure 9. Comparison between the temperature simulated by SUSHI and the temperature measured by thermocouple on the façade of LC4 building.

Figure 10 .

 10 Figure 10. Comparison between the radiance measured by thermography (a) and the radiance simulated by SOLENE/SUSHI/MOHICANS (b). Radiance contrast for the image (c) and radiance profile along the façade of LC4 building (d).

Table 1 .

 1 Main existing thermo-radiative urban models.

Table 2 .

 2 Wall characteristics.Figure1. Flow chart of the simulator

Table 1 .

 1 Main existing thermo-radiative urban models.

	Model: autors,	Groleau et al.	Krayenhoff, Voogt	Asawa, Hoyano	Kastendeuch, Najjar	Yang, Li
	(NAME), year	SOLENE [21] (2003)	TUF-3D [18] (2007)	[17] (2008)	[16] (2008)	MUST [20] (2013)
	3D scene	From SketchUp	3D cells organized in cluster geometry	3D-CAD with trees	3D	3D cells organized in cluster geometry
	Direct solar	Measure or Perrin de Brichambaut model	Iqbal's model (83)	Bouguer's law	Measure + position by Reda Afshin	Bouguer's law
	Diffuse solar	Measure + "all weather" Perez sky	Iqbal's model (83)	Nagata's formula	Measure + "all weather" Perez sky	Berlage's relation
	Atmospheric radiation	Measure + sky view factor	Prata's formula (96)	Brunt's formula	Measure + homogeneous	Model of Berdahl and Martin (84)
	Shape factor	By subdivision in triangular facets	Exact plane parallel analytical equation	Deduced from sky view factor	Pianykh's method	By a method similar to Monte Carlo
	Multiple reflections VNIR	Progressive refinement	Two reflections minimum	One	Progressive refinement	Gebhart factor
	Multiple reflections IR	Progressive refinement (not with sky)	Isotropique sky radiance	CN + other surface at Tair	Progressive sky) refinement (not with	Gebhart factor
			Function of surface			Function of surface
	Sensible heat flow Imposed h	roughness and	Jürges's formula	h by Ito (78)	roughness and
			effective wind speed			effective wind speed
	Wall thermal model	Nodal (5 nodes)	Finite differences + iterative process with IR radiation	Finite differences (5 days)	Finite differences + IR radiation iterative process with	Finite differences
	Main limitation	Limited number of nodes	Simplified geometry (3D Cartesian)	Simplified reflection of VNIR and IR fluxes	3D scene and wall discretization	Simplified geometry (3D Cartesian)

Table 2 .

 2 Wall characteristics.

	Layer	e (m)

Acknowledgements

The authors thank PACA Region for its financial support, the laboratory CERMA for the provision of SOLENE software, the laboratory IUSTI for the provision of the IR MWIR camera and the ONERA staff who participated to the measurement campaign BATIR, in particular M. P.-Y. Foucher for his measurements of spectral reflectance.


G Guillonneau
email: gaylord.guillonneau@empa.ch

G Kermouche

J Teisseire

E Barthel

S Bec

J.-L Loubet

Is the second harmonic method applicable for thin films mechanical properties characterization by nanoindentation?

Keywords: nanoindentation, films, second harmonic, numerical simulation, contact depth '

The second harmonic method is a dynamic indentation technique independent of the direct indentation depth measurement. It can be used to determine near-surface mechanical properties of bulk materials more precisely than classical dynamic nanoindentation. In this paper, the second harmonic method is extended to the measurement of the mechanical properties of thin PMMA layers deposited onto silicon wafers. It is shown that this new technique gives precise results at small depths (less than 100nm), even for films with a thickness lower than 500nm, which was not possible to achieve with the classical CSM method. However, experimental and numerical results obtained both with classical nanoindentation and second harmonic methods differ at high indentation depth. Using FE simulations and AFM measurements, it is shown that the contact depth calculation with classical models can explain this difference.

Introduction

About nanoindentation measurement

Because the nanoindentation is a non-destructive technique which measures the mechanical properties at the nanometer scale, it is particularly adapted to probe surface properties of very thin films. The principal properties measured with this apparatus are the reduced elastic modulus E '* and the hardness H of the tested material [1,2]:

c A P H = (1) 1 1 1 i c E E E - = (2)
With:

c c A S E π 2 '* = ( 3
)
Where P is the load, A c is the projected contact area, S is the harmonic contact stiffness, '* i E is the reduced elastic modulus of the tip, and '* c E is the reduced contact modulus between the sample and the tip. The contact area A c is directly linked to the contact depth h c by the following geometrical equation (for a perfect pyramidal indenter) [3]:

( ) 2 2 tan c c h A θ π = ( 4
)
Where θ is the equivalent semi-angle of the tip. This expression can be changed in order to take into account the non-perfect geometry of the indenter [3]. The sample stiffness S can be measured by fitting the unloading curve [3], permitting to obtain E and H, at maximum load, or with the dynamic CSM (Continuous Stiffness Measurement) technique. The CSM method consists of adding a small oscillation to the load signal in order to measure the dynamic contact stiffness as a function of the indentation depth [4,5]. This dynamic technique is particularly used for thin films characterization, because the mechanical properties of the layer can be measured as a function of the indentation depth with only one test, which is impossible with the quasi-static method.

About indentation of thin films and related difficulties

Thin film indentation differs from bulk material characterization, and some difficulties have to be considered. First, the substrate properties have to be taken into account. Indeed, during thin film indentation, the measured properties are a combination of the film properties and the substrate properties. To account for the substrate effect, several literature models, more or less complex, are proposed. Nix proposes a simple model to determine the film modulus, dependent on an adjustable parameter [6]. Rar et al. developed a more complex model, requiring the knowledge of the Poisson ratio of the tested materials [7,8]. Perriot and Barthel proposed another elastic model to take into account the substrate effect. This model can be expressed by the following equation [9]:

( ) n c f s f h tx E E E E         + - + = θ tan 1 0 '* '* '* '* (5)
Where

( ) ( )         +         + = - '* '* 1 '* 1 tan 2 tan 2 1 1 s c f c E h t E h t E θ π θ π (6)
This model gives good results for a large range of materials and is not dependent of adjustable parameters.

The second difficulty, common to bulk materials, comes from the contact depth calculation. The contact depth value cannot be directly measured, but can be calculated with models (in order to model pile-up or sink-in phenomena) using the displacement measurement h. The most popular is the Oliver and Pharr model which links the contact depth h c to the indentation depth h by the following equation [3,5]:

S P h h c ε - = (7)
Where ε is a parameter which depends on the tip geometry (ε=0.75 for a Berkovich indenter). Another expression proposed by Loubet et al., can be used [12,13]:

      + - = 0 h S P h h c α (8)
Where α=1.2 for a Berkovivh indenter and h 0 is a coefficient calculated to take into account the tip defect [13,14].

The third difficulty comes from the uncertainties linked to the displacement measurement. As explained in a previous paper, the displacement measurement is influenced by uncertainties like contact point detection, tip defect, thermal drift, or sample roughness [14,15]. These uncertainties are particularly important at very small penetration depths, resulting in less precise mechanical properties measurement at very small indentation depths, which can be problematic when the film thickness is lower than 100 nm.

What could we expect from the use of the second harmonic method?

The second harmonic method can be used to measure the mechanical properties of bulk materials at indentation depths lower than 100 nm with better accuracy than classical CSM method because the technique is independent of the direct displacement measurement [16].

Consequently, the method is expected to give precise values of mechanical properties at small indentation depth for thin films. In this paper, an application of the second harmonic method on thin PMMA layers deposited on silicon wafer is presented. In the following part, the expressions used to calculate thin films properties with the second harmonic method are detailed. In the third part, the experimental tests, the samples and the apparatus are described.

The results are shown and discussed in the fourth part, and we discuss the difference between results at high indentation depths in the last part.

Second harmonic method: substrate effect

In a previous paper, it was shown how to measure the mechanical properties of bulk materials with the second harmonic method [16]. In the case of thin film indentation the same approach can be used, but the calculation of dS/dh (equation (11) of [16]) has to consider the reduced contact modulus variation as a function of the indentation depth, because of the substrate's elasticity influence. The Sneddon relation links the stiffness to the reduced contact modulus and the contact depth by the following equation:

( ) c c h E S θ tan 2 '* = (9)
Considering the variation of '* c E and h c with the indentation depth, and using Equations ( 3) and ( 9), the derivative of Equation ( 9) with respect to the displacement gives:

( ) ( ) dh E d SE dh dh E dh dS c c c '* '* '* / 1 tan 2 - = θ (10)
The expression depends on the indenter geometry, the reduced contact modulus, dh c /dh, and d(1/E '* )/dh. The derivative of the contact depth with respect to the indentation depth can be determined by a simple local derivative of the h c -h curve. However, the term d(1/E '* )/dh is not known. We propose the use of a literature model to calculate this expression. Because of its simplicity, and also because it does not depend on adjustable parameters, the Bec et al. model was used to calculate d(1/E '* )/dh [10,11]. After derivation of Equation ( 6), the term d(1/E '* )/dh can be obtained:

( ) ( ) ( ) dh dh h t E E h t dh E d c c f s c 2 '* '* 2 '* tan 2 1 1 1 tan 2 / 1         +         - = θ π θ π (11)
In this expression, all terms can be measured or calculated. Finally, with Equations ( 10) and ( 11), the expressions of '* c E and H can be obtained:

( ) ( ) dh E d S dh dh dh dS E c c '* '* / 1 tan 2 - = θ (12) 2 '* 2 4 c E S P H π = (13)
And E '* is calculated using Equation (2).

Experimental details

Tested samples

Experimental tests were performed on poly(methyl methacrylate) (PMMA) (Sigma-Aldrich, St Quentin Fallavier) [Mw=15 kg/mol, entanglement molecular weight Me = 10 kg/mol].

Thin layers were deposited by spin-coating onto clean silicon wafer Si(100) (Neyco). Five samples were tested with different PMMA film thicknesses t, adjusted by dilution in methylisobutylketone (185 nm, 348 nm, 631 +/-3 nm, 982 +/-3 nm and 1914 +/-3 nm, thicknesses measured by AFM). The residual stresses and residual solvent in the film were minimized by annealing at 170 °C for 20 min and cooling down to room temperature in 4 h. The samples were cut, then glued with cyanoacrylate on an aluminum sample holder, and cleaned with ethanol using a standard paper towel.

Apparatus

A SA2® Nanoindenter equipped with a DCM head was used in the experiments. The resolutions in force and displacement are 1nN and 0.2pm respectively. The apparatus and the second harmonic measurement method were described in more details in a previous paper [14]. A sharp diamond Berkovich tip was used in the experiments (h 0 =5 nm, calculated using the Loubet et al. method [13]). For each sample, 10 tests were done to obtain a mean value of the properties. The maximum load was 10 mN. Because the PMMA is a material whose properties are time dependent, a constant strain rate was applied with P P / & =0.01s -1 [17,18].

The CSM technique was applied with an oscillation amplitude varying between 2 and 4 nm following the indentation depth value [19], and a frequency of 31 Hz.

Mechanical properties are computed using both the classical CSM method and second harmonic method and are plotted as a function of the ratio a/t (contact radius/layer thickness).

Results for depths lower than 25 nm are not plotted in the paper, the tip defect being too much important at lower depths to measure precise mechanical properties. The theoretical values of the reduced elastic modulus used in Equation ( 11 GPa for silicon [10,14]. The contact depth used in Equation ( 2) is calculated with the Loubet et al. model [12,13].

Results

Hardness

From Figures 1 and2, the hardness of PMMA layers is approximately constant for small a/t values (a is the contact radius) with the classical CSM method and second harmonic method, and is in good agreement with the PMMA hardness [20]. With the second harmonic method, the hardness curves are superposed at very small a/t ratio. It is not the case with the classical CSM method, as some scattering in hardness value is observed at low a/t ratio. This scattering at low depths is caused by of the uncertainties related to the tip displacement, especially the contact point detection and the tip defect. This confirms that the second harmonic method is particularly adapted to the measurement of mechanical properties at small indentation depths, given that it is not dependent of the displacement measurement. An increase of the hardness can be observed for CSM method from a/t = 0.8, showing the substrate effect on the measurement. The scattered observed in the curves is due to the stiffness measurement which is scattered at this range. With the second harmonic method, the hardness increase begins for smaller a/t ratio (a/t=0.2), for samples with thickness between 631 and 1914nm. This difference will be explained in section 5.2.

Elastic modulus

On Figures 3 and4, the experimental reduced elastic modulus determined with the two techniques is compared to values calculated with the Bec et al. and the Perriot and Barthel models. Good correlation between experimental and numerical results is observed for small a/t ratios in both methods. Furthermore, as obtained on hardness results, the elastic modulus measured with the second harmonic method is more precise than the modulus determined with the classical CSM technique for the same reason than in the previous paragraph.

However, for higher a/t ratios, the experimental elastic modulus is higher than calculated values. With the CSM method, the curves differ from a/t=1, and with the second harmonic method, the curves differ from a/t=0.4. Moreover, the elastic modulus calculated with the second harmonic method is twice higher than the modulus determined with the CSM method.

A first explanation of this divergence can come from the contact depth calculation. Indeed, with their numerical model, Perriot and Barthel show a strong difference between the real contact depth and h c calculated with the literature models when the

'* '* / f s E E
ratio is very high or very low [9]. The contact depth is calculated supposing a homogeneous material (which is not the case for thin films indentation) and Perriot and Barthel explained that this hypothesis is the cause of this divergence.

The CSM method is based on the calculation of h c , and the second harmonic method is based on h c and dh c /dh determination. Consequently, if h c is not correctly calculated, its derivative also can be not correctly calculated. This explains the large difference between experimental and numerical results. In the following paragraph, a Finite Element simulation of the PMMA/Si indentation and AFM measurements are presented in order to further investigate the contact depth calculation.

Finite Element analysis

Finite Element model

Calculations were performed with Systus/Sysweld [21,[START_REF]Systus/Sysweld. User's Manual[END_REF] using axisymmetric elements and a large displacement / large strain option (updated Lagrangian formulation). Here the commonly-used assumption in which the Berkovich tip can be replaced by a cone of semi angle 70.32° is considered. The Si thickness has been chosen to be one hundred times the maximum contact area. To ensure plastic incompressibility, four node quadrilateral isoparametric elements with a selective reduced integration scheme are used in the plastically deformed area. The plastic flow is described via a plastic von Mises stress. The loading is achieved by imposing a quasi-static displacement of the indenter.

One important difficulty related to the Finite Element analysis of the indentation of soft thin films over hard substrate is the large mesh distortion of the film at high penetration depth and the necessity to use very small elements at low penetration depth. The only way to bypass this issue is to use an adapted remeshing procedure. For that purpose, we developed an automatic remeshing procedure based on the work of Kermouche et al. in the case of scratch testing of coating/substrate system [START_REF] Kermouche | [END_REF]. The remeshing strategy is based on the definition of several zones in which the mesh density is different. A very fine mesh is defined near the contact. In the other zones, the greater the distance to the fine mesh zone, the lower the mesh density. During a remeshing stage, these zones are calculated starting from the coordinates of the indenter and the parameters of the contact area (contact radius, contact depth). The meshing algorithm is such that the quantity of nodes in contact is maintained constant during the calculation and thus does not depend on the penetration depth of the indenter. The node quantity in the film thickness is computed from the size of the contact area. Indeed, when the thickness of the film is negligible compared to the contact radius, it is not necessary to describe finely what happens in the film. In the opposite case, when the contact radius is lower than the film thickness, an important elements quantity is required, because the deformation of the film may play a very significant role. Therefore, it is possible to accurately determine the contact area and the contact stiffness whatever the indentation depth by superimposing small oscillation amplitude to the imposed displacement in order to simulate the Continuous Stiffness Measurement method.

The diamond and silicon are supposed to follow a linear isotropic elastic behavior, and the PMMA layer is supposed to be elastic-perfectly plastic. Note that the well-known timedependence of PMMA mechanical properties [18,21,24] is not taken into account to simplify the model and the analysis. This will be a future issue of this work. The elastic modulus of diamond is * i E =1150 GPa. The PMMA yield stress is y σ =100 MPa [21,25].

Results

On Figure 5, the composite reduced elastic modulus determined by Finite Element simulation is in good agreement with values obtained by literature models, at both small and high a/t ratios. It is satisfying because the modulus was calculated only from the basis of computational data and was not affected by the use of any model to estimate the contact depth.

On Figure 6, the h c /h ratio determined by Finite Element simulation is compared with values calculated with the Loubet model, the Oliver and Pharr model, and the Perriot and

Barthel model [3,9,12,13]. These three models were computed from the contact stiffness, indentation depth and load extracted from the Finite Element analysis. The h c /h ratio calculated with the Loubet and the Oliver and Pharr models are almost independent of the penetration depth, contrary to the h c /h ratio calculated from the Finite Element simulation.

From a/t>1.5, this h c /h ratio is significantly higher than values obtained with classical models, which can explain why the experimental elastic modulus determined with CSM method is overestimated. Moreover, the h c /h ratio is higher than 1 from a/t>1.5, which means that pileup occurs around the indenter as shown in Figure 8. This is confirmed by the residual print morphology measured with the help of an Atomic Force Microscope (Figure 9).

Consequently, the contact depth cannot be calculated correctly with the Oliver and Pharr model. It can also be observed that the h c /h ratio obtained with the Perriot and Barthel model is the lower one. This result is not surprising because this model is based on the assumption that both film and substrate follow a linear isotropic elastic behavior. Note however that its variation follows the same trend of the Finite Element results.

As the second harmonic method depends also on the derivative of the contact depth with respect to the indentation depth, the variation of dh c /dh versus a/t ratio was plotted in the combination of the errors in the computation of h c /h and dh c /dh makes the second harmonic method more inaccurate than the CSM method in the case of film/substrate system for high a/t ratio. More precisely, the ratio a/t below which the second harmonic method can be applied is lower than the one corresponding to the CSM method as observed on Figures 3 and4.

All these results show the importance of the correct calculation of the contact depth for thin film indentation. However, this error in the calculation of h c and dh c /dh is not sufficient to totally explain the difference. Another phenomenon linked to the dependence of elastic properties of PMMA to the hydrostatic pressure can be considered [26,27]. A Finite Element simulation of the PMMA/Si indentation with elastic properties of PMMA dependent on the hydrostatic pressure could be a perspective complementary to this work.

Conclusion

In this paper, the second harmonic method was applied to thin PMMA layers deposited onto silicon wafer. To calculate correctly the mechanical properties, it was necessary to recalculate the derivative of the contact stiffness with respect to the displacement in order to account for the variation of the elastic modulus versus indentation depth. This calculation was performed using the model of Bec et al. [10,11]. Experimental application shows precise properties measured at small indentation depths with both CSM and second harmonic methods, with a better accuracy for second harmonic method at low a/t ratio. At high indentation depth, a significant difference between experimental and numerical results was observed. A Finite Element investigation coupled with AFM observation of indentation print evidenced that classical contact models fail to compute accurately both h c /h ratio and dh c /dh leading thus to an overestimation of the mechanical properties.

Finally, one important conclusion of this paper is that both the CSM method and the second harmonic method are not accurate enough to measure the mechanical properties of film/substrate system over a large range of a/t ratio using standard contact models. With the second harmonic method, it is possible to measure accurately the mechanical properties at low a/t values, where the CSM method is subjected to various measurement uncertainties. On the contrary, with the CSM method, it is possible to obtain accurate results over a larger range of a/t ratio (until a/t=1). The reason why the CSM method fails at larger indentation depth may be linked to the bad estimation of the contact area but it can also be related to the possible hydrostatic pressure dependence of the elastic properties of PMMA when it is confined between the tip and the silicon substrate. This last point will be considered in a future paper.

Figures

  ) are '* f E = 5 GPa for PMMA and '* s E =180

Figure 7 .

 7 Figure 7. The same trends are observed. The three contact models fail to compute accurately

Figure 1 .

 1 Figure 1. Hardness versus a/t ratio (contact radius/film thickness) for PMMA/Si samples, calculated with the CSM method. Each curve corresponds to the mean value from 10 indentation tests.

Figure 2 .

 2 Figure 2. Hardness versus a/t ratio (contact radius/film thickness) for PMMA/Si samples, calculated with the second harmonic method. Each curve corresponds to the mean value from 10 indentation tests.

Figure 3 .

 3 Figure 3. Composite reduced elastic modulus versus a/t ratio (contact radius/film thickness) for PMMA/Si samples calculated with the CSM method. Each curve corresponds to the mean value from 10 indentation tests. Comparison with modulus calculated with literature models.

Figure 4 .

 4 Figure 4. Composite reduced elastic modulus versus a/t ratio (contact radius/film thickness) for PMMA/Si samples, calculated with the second harmonic method. Each curve corresponds to the mean value from 10 indentation tests. Comparison with modulus calculated with literature models.

Figure 5 .

 5 Figure 5. Composite reduced elastic modulus versus a/t ratio (contact radius/film thickness) for PMMA/Si sample, determined by Finite Element simulation. Comparison with modulus calculated with literature models.

Figure 6 .

 6 Figure 6. h c /h ratio versus a/t ratio (contact radius/film thickness) for PMMA/Si sample obtained by F.E. simulation. Comparison between the value calculated from F.E. simulation and h c /h calculated with literature models using the simulation data as input.

Figure 7 .

 7 Figure 7. dh c /dh versus a/t ratio (contact radius/film thickness) for PMMA/Si sample obtained by F.E. simulation. Comparison between the value obtained with F.E. simulation and dh c /dh calculated with literature models.

Figure 8 .

 8 Figure 8. Von Mises stress distribution near the contact during indentation on PMMA layer from EF simulation.

Figure 9 .

 9 Figure 9. a) AFM image (error mode) of an indentation print on a PMMA/Si layer (t=982 nm). b) Topography profile along the line.

Figure 1 .

 1 Figure 1. Hardness versus a/t ratio (contact radius/film thickness) for PMMA/Si samples, calculated with the CSM method. Each curve corresponds to the mean value from 10 indentation tests.

Figure 2 .

 2 Figure 2. Hardness versus a/t ratio (contact radius/film thickness) for PMMA/Si samples, calculated with the second harmonic method. Each curve corresponds to the mean value from 10 indentation tests.

Figure 3 .

 3 Figure 3. Composite reduced elastic modulus versus a/t ratio (contact radius/film thickness) for PMMA/Si samples calculated with the CSM method. Each curve corresponds to the mean value from 10 indentation tests. Comparison with modulus calculated with literature models.

Figure 4 .

 4 Figure 4. Composite reduced elastic modulus versus a/t ratio (contact radius/film thickness) for PMMA/Si samples, calculated with the second harmonic method. Each curve corresponds to the mean value from 10 indentation tests. Comparison with modulus calculated with literature models.

Figure 5 .

 5 Figure 5. Composite reduced elastic modulus versus a/t ratio (contact radius/film thickness) for PMMA/Si sample, determined by Finite Element simulation. Comparison with modulus calculated with literature models.

Figure 6 .

 6 Figure 6. h c /h ratio versus a/t ratio (contact radius/film thickness) for PMMA/Si sample obtained by F.E. simulation. Comparison between the value calculated from F.E. simulation and h c /h calculated with literature models using the simulation data as input.

Figure 7 .

 7 Figure 7. dh c /dh versus a/t ratio (contact radius/film thickness) for PMMA/Si sample obtained by F.E. simulation. Comparison between the value obtained with F.E. simulation and dh c /dh calculated with literature models.

Figure 8 .

 8 Figure 8. Von Mises stress distribution near the contact during indentation on PMMA layer from EF simulation.

Figure 9 .

 9 Figure 9. a) AFM image (error mode) of an indentation print on a PMMA/Si layer (t=982 nm). b) Topography profile along the line.