Joël Puibasset
email: puibasset@cnrs-orleans.fr

.-M Pellenq
email: pellenq@crmcn.univ-mrs.fr

A Grand Canonical Monte Carlo simulation study of water adsorption on Vycor-like hydrophilic mesoporous silica at different temperatures

Keywords: number: 68.43.De, 61.43 Gt, 61.43.Fs

A Grand Canonical Monte Carlo simulation study of water adsorption on Vycor-like hydrophilic mesoporous silica at different temperatures

Joël Puibasset, Roland

INTRODUCTION

The water/silica system is very important in a wide range of research fields, geophysics, pharmaceutics, industry environment. Most of the evoked processes, water infiltration, depollution, phase separation, catalysis, etc. involve nanometric scale physics. This explains the intense theoretical, experimental and simulation activity aiming at getting insight on the water nature at the silica interface. The thermodynamic, dynamic and structural characteristics are actually expected to differ considerably from bulk ones as interfacial properties dominate the physical processes at work at the mesoscale. [START_REF] Gregg | Adsorption, Surface Area and Porosity[END_REF][START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF][START_REF] Korb | [END_REF][4][5][6][7][8][9][10][11][12][13] The capillary condensation is a striking example of confinement effect on the thermodynamic properties of a fluid. Structural modifications have also been observed by neutron scattering experiments. For instance, it has been shown that water structure is distorted in the vicinity of silica surfaces; distortion that may explain why water freezes in cubic ice structure when confined in Vycor, a mesoporous silica glass. 14 More recently, Bruni et. al. 15,16 have shown how neutron diffraction experiments with isotopic substitution are able to give improved sitesite radial distribution functions for water confined in Vycor glass. The results show that water is still hydrogen bonded, but with a distorted network compared to bulk. Molecular Dynamics simulation was used to help understanding these results. Previous simulation works on water confinement focused on simple silica substrate geometry, slit-like 7,8 or cylindrical 9,10 pores, supposed to mimic real porous solids. However, such ideal pores are a rough modeling of real pores, since they introduce only one characteristic pore size and miss the pore size dispersion and interconnections between pores, or in other words the topology of the system. The present work aims at addressing both the thermodynamics and the structure of water confined in a realistic mesoporous silica glass, presenting both a morphological (pore size distribution) and a topological (interconnections) disorder. Another originality of this work lies in the statistical ensemble that we explore (Grand Canonical Ensemble) which differs from previous studies of confined water based on Molecular Dynamics.

DESCRIPTION OF THE MESOPOROUS SILICA SUBSTRATE

An off-lattice reconstruction algorithm has been used to generate the atomistic vycor-like mesoporous structure considered in this study. Portions of an initial cubic crystal of cristobalite (106.95 Å of edge) are cut out by applying the off-lattice functional representing the gaussian field associated to the volume autocorrelation function of the studied porous structure. [17][18][19] Periodic boundary conditions are applied to the gaussian field, so as to generate the periodic structure shown in figure 1. It has been shown 20 that this numerical sample reproduces quite well micro-and meso-textural properties of real Vycor, and constitutes an improvement in the topological description of silica mesoporous glasses, because of the presence of disordered and interconnected pores. To have a realistic description of the surface chemistry, all silicon atoms in an incomplete tetrahedral environment are removed, and all oxygen dangling bonds are saturated with hydrogen atoms placed 1 Å from oxygen perpendicular to the surface. The porosity of the sample is 0.28, the density 1.56 g/cm [START_REF] Korb | [END_REF] , the specific surface area 210 m 2 /g, the mean pore diameter 36 Å (from chord distribution), and OH surface density 7 OH/nm 2 , all in agreement with experimental values: 0.28, 1.50 g/cm [START_REF] Korb | [END_REF] , 200 m 2 /g, and around 7 OH/nm 2 if out-gassed around 400°C 21,22 , except for the pore size slightly smaller in our pseudo-vycor compared to real one (40 Å). It has been shown that this numerical sample reproduces quite well micro and meso textural properties of real vycor. [17][18][19]

INTERMOLECULAR POTENTIALS

The SPC model [START_REF] Berendsen | Intermolecular Forces[END_REF] has been used to describe the water-water interaction, because it is a fast computable model well suited for very large systems (the completely saturated system contains more than 10000 water molecules). This model reproduces well the thermodynamic and structural properties around ambient temperature, like vapor pressure, enthalpy of vaporization [START_REF] Jorgensen | [END_REF] , and radial distribution functions. In this model, the water molecule is rigid with the geometry defined in fig 2 . The molecules interact through electrostatics with the partial charges given in fig 2, and through usual quantum interactions modeled by a Lennard-Jones like term between the water molecules. Since the main contribution comes from the oxygen atoms, the Lennard-Jones potential is calculated from the oxygen-oxygen distance, and with the parameters given in fig2.

The interaction of water with Si, O and (surface) H atoms forming the vycor structure is assumed to remain weak, in the physisorption energy range. In this work, we have used a TrAZ form of the original PN-type potential function as reported for adsorption of rare gases and nitrogen in silicalite-1 25 . The PN-TrAZ potential function is based on the usual partition of the adsorption intermolecular energy restricted to two body terms only: it includes a dispersion interaction term, a repulsive short range contribution and an induction term. Three body dispersion (XYY) terms (where X is the adsorbate and Y an adsorbent species) can be obtained in principle from the same approach 26 but are neglected in the present study since they are expected to be significant for water molecules in contact with the silica surface (which is a small number of the total adsorbed amount). The choice of this particular model to describe the water/adsorbent potential was motivated by the good degree of parameter transferability. Indeed, in a previous study 31 , we found that using a set of potential parameters previously derived for adsorption in silica zeolite 25 augmented to take into account hydroxyl groups, the TrAZ model allows reproducing both low coverage experimental adsorption isotherm (amount adsorbed versus pressure at constant temperature) and isosteric heat curve with no further adjustment. In the TrAZ model, the adsorbate-surface energy (

A k i u ∈ ) of the k th
atom of molecule A at a position i in the simulation box, is given by

{ } 2 i E k H , Si , O j n 2 ij kj , n 2 5 3 n n 2 ij r kj b kj A k i 2 1 r C f e A u α -         - = ∑ ∑ ∈ = - ∈ (1)
Where r ij is the distance between the matrix species j and the position i of the k th atom of a given molecule, f 2n are so-called damping functions that depend on interspecies distance and on a repulsive parameter (see below) and finally E is the electric field created at adsorbate position i due to the inorganic species ; all other symbols are physical parameters that depend on electronic properties of interacting species only. The sums runs over all (j) atomic sites in the matrix that are oxygen, silicon and hydrogen atoms. The first term in the sum is a Born-Mayer term representing a two-body form of the short range repulsive energy due to finite compressibility of electron clouds when approaching the adsorbate at very short distance of the pore surface. There is one such a term per pair of interacting species. The repulsive parameters (A kj and b kj ) are obtained from mixing rules of like-atoms pairs (see below). The second term in the above equation is a multipolar expansion series of the dispersion interaction in the spirit of the quantum mechanical perturbation theory applied to intermolecular forces. [START_REF] Stone | The Theory of Intermolecular Forces[END_REF] It has been shown that two (and three body) dispersion C 2n coefficients for isolated or in-condensed phase species can be obtained from the knowledge of the dipole polarizability and the effective number of polarizable electrons N eff of all interacting species 25,26 which are closely related to partial charges that can be obtained from ab initio calculations. The f 2n terms in the above equation are damping functions of the form:

( ) ij r kj b n 2 0 m m ij kj n 2 e . ! m r b 1 f - = ∑         - = (2)
The role of these damping functions is to avoid divergence of the dispersion interaction at short distance where the wave functions of the two species overlap (i.e. when the interacting species are at contact) [START_REF] Tang | [END_REF] . They allow taking into account the possible electronic exchange which has a non zero probability at short distance even for two close shell structures. For each pair of interacting species, they are parametrized with the single b kj repulsive parameter. The damped dispersion multipolar expansion can be seen as a convenient way to make the perturbation theory valid at short inter-atomic separations. The last term in equation ( 1) is the induction interaction as written in the context of the quantum mechanical perturbation theory applied to intermolecular forces. [START_REF] Stone | The Theory of Intermolecular Forces[END_REF] It represents an attractive energy arising from the coupling of the polarisable electronic cloud of the adsorbate atom of polarizability α k at position i with the electric field E i induced by the charges carried by framework species (O, Si and H) that result from the bonding process within the inorganic matrix itself. In the case of water interacting with vycor, one has to parametrize six different adsorbate/adsorbent-species potentials; all atomic parameters and coefficients are given in Table 1 and2. Note that repulsive parameters for like-pairs are taken from a previous work on the simulation of rare gases in silicalite 25 using Bohm and Ahlrichs 29 combination rules based on ab initio results at the Hartree-Fock level approximation which only describes the repulsive energy term in the case of non charged and polar species. This type of potential function based on the PN-TrAZ parametrization method was used in various studies of molecular and covalent fluids at interfaces from open surfaces 30,31 to microporous zeolites 32,33,34,35,36 and more recently in the case of mesoporous vycor-like materials 20,37,[START_REF] Puibasset | Proceedings of the 6 th Int. Symposium on the Characterization of Porous Solids[END_REF]31,39,40,41 .

Minimal image convention is adopted to calculate all interactions. Long range mean-field correction to dispersion terms cannot be precisely calculated because the density is not uniform. Actually, the size of the box is large enough (107 Å) so that large distance contributions are negligible. The electrostatic contribution is evaluated by summing on neutral subgroups of atoms of highest symmetry (tetrahedral silicon SiO 4 and surface hydroxyls OH).

Implementation of Ewald summation procedure 42 has proven to be of little improvement, probably due to the large box size, absence of isolated charges, and negligible total dipolar moment of the substrate (due to isotropy in normal orientation on Vycor surface).

Calculations show that the largest part of the substrate contribution to electrostatics comes from the surface dipoles (hydroxyls), the tetrahedral SiO 4 giving rise to weak quadrupolar term only.

In order to save computing time, the water/substrate interactions have been calculated in advance and saved on a grid. This is possible since the substrate is kept rigid in this study.

The grid mesh has been chosen equal to 1 Å. Due to memory constraint it could not be made smaller. However, close to the silica surface, this mesh is not fine enough for small dipolar molecule like water. This is why the contribution of the substrate has been divided in two parts, inside and outside a sphere of 6.5 Å in radius, centered on the grid cell to be evaluated.

Only the contribution from outside is saved on the grid, the rest being calculated directly during the course of the simulation. A neutrality constraint has been added so as to avoid fast oscillations of electrostatic contributions due to isolated charges coming in and out of the 6.5 Å radius sphere: for instance O and H atoms from the same surface hydroxyl group are taken together. The choice of the radius is a compromise between CPU time saving and grid smoothness. Three different grids are necessary: two dispersion-repulsion grids for both O and H atoms of water, plus one for electrostatic potential. Since there are eight corners for each cubic cell, and since the 6.5 Å radius sphere has to be centered on the cell, the total number of grid used is 8 time 3, and interpolation in each cell is made between 8 grids. All grids are defined everywhere in the simulation box due to the short distance cutoff. The result obtained with this method for the whole interaction is always within 4% of the value obtained by a completely direct calculation with all the atoms of the simulation box. Our method is then very accurate, and very fast (few tens of atoms taken into account instead of the 60 000 species that are present in the silica texture of our vycor numerical sample). This procedure is well suited to study adsorption properties of polar molecules and allows considering large systems.

The water-water interaction calculations become excessively slow at low temperature when the substrate is almost completely filled (more than 10000 water molecules). It was then decided to implement a hierarchy method 43 . The simulation box is divided in a 3-dimensional array of 5 A wide cells. The total interaction energy of a water molecule with the others is the sum of the individual contributions of the molecules contained in the neighboring cells (within a distance of 15 A) and the multipole representation for the other cells. The uncertainty introduced by this method do not exceed few percent.

WATER ADSORPTION ISOTHERMS

The water adsorption isotherms have been calculated using the Grand Canonical Monte Carlo Simulation [START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF] (GCMC), which mimics a real adsorption experiment where the temperature, the volume and the chemical potential (related to the pressure of the gas) are kept constant. An equal number of translation, rotation, creation or destruction of molecules has been chosen, around 2.10 4 per water molecule at full saturation, up to 10 5 at low coverage. The isotherms have been calculated for four different temperature, 300 K, 350 K, 400 K and 500 K. For each temperature, the pressure of the vapor in the particle reservoir is normalized to the saturating vapor pressure for SPC model [START_REF] Errington | [END_REF]46 given in Table 3. The adsorbed quantity is normalized to the known geometric surface of the numerical sample (401 nm 2 ) 37 , and expressed in µmol/m 2 . The result is shown in figure 3. As can be seen, the isotherms show strong adsorption at low pressure as expected for hydrophilic surfaces. Above a relative pressure of 0.7 the adsorption increases faster: this is the so-called capillary condensation. However, as can be noted, the full saturation is generally reached above the saturating vapor pressure. This point will be discussed later.

The energetics of the adsorption has been analyzed by calculating the isosteric heat of adsorption q st on the silica substrate from cross-fluctuations in energy (U) and adsorbed quantity (N) according to the formula 44

RT N N N U UN q 2 2 st + - - - = (3)
The values for the four temperatures are given as a function of water adsorbed amount in figure 4. The bulk SPC water latent heat 46 corresponding to each temperature is also represented as an horizontal line (see Table 3). As can be seen for all temperatures the heat of adsorption is significantly higher than the latent heat which is characteristic of an hydrophilic surface. The initial heat of adsorption obtained for very low pressure reach 80 kJ/mol at 300 K, which is 1.8 times the bulk latent heat of water, and 65 kJ/mol at 500 K which is 2.1 times the latent heat of water at the same temperature. These values are in agreement with experimental determinations 47 , which validate the intermolecular potential used in this work.

Inspection shows that the isosteric heat of adsorption has reached its bulk value as soon as the monolayer coverage is reached (around 16 µmol/m 2 at 300 K). The second layer and subsequent ones are then expected to have properties very similar to the bulk. Furthermore, it can be inferred from these data that the sites are energetically heterogeneous. This point will be reconsidered during the BET analysis.

The Brunauer Emmett and Teller (BET) 48,1,2 model has been applied to the water adsorption isotherms previously calculated to evaluate the monolayer coverage and energetics of the system. In this model the first layer adsorption is supposed to take place on energetically uniform substrate sites with an energy ε close to the heat of adsorption, and the next layers are adsorbed on the previous ones with a different energy: L (supposedly closer to the latent heat of the bulk fluid). Assuming this hypothesis, it is shown that the adsorbed amount n (for instance in µmol/m 2 ) and the relative pressure 0 /P P can be cast into the following equation:

0 0 0 0 0 / 1 1 ) / 1 ( / P P C n C C n P P n P P - + = - (4)
where 0 n is the monolayer capacity expressed in the same units as n (µmol/m 2 ), and

kT L e C / ) ( - = ε
is related to the difference between heat of adsorption on the substrate and bulk latent heat. Figure 5 shows the corresponding left hand quantity as a function of reduced pressure. As can be seen the curves show a linear portion between relative pressure of 0.1 and 0.4. However, the low pressure adsorption do not follow the BET model. This was expected since the isosteric heat of adsorption do not show a sharp variation between low pressure, where it is close to the adsorption site energetics, and higher pressure where the isosteric heat of adsorption is close to the bulk fluid latent heat. The linear regression of the curves between relative pressure of 0.1 and 0.4 have been determined and drawn in figure 5. Their analysis gives the monolayer capacity 0 n and the difference between adsorption energy and bulk latent heat. The results are summarized in Table 3. As can be seen, the monolayer capacities at different temperatures are close to each other, except for the highest temperature which probably corresponds to conditions for which the localized adsorption assumption at the heart of the BET model is no longer valid. The obtained value for the monolayer capacity is around 7.5 µmol/m 2 , which corresponds to 22 Å 2 per water molecule. This value is much larger than the commonly admitted surface of 10.5 Å 2 per molecule deduced from bulk water density, or even more refined and larger values proposed by different authors (12.5 Å 2 49 , or 14.8 Å 2 50,51 ), to account for the low specific surface area obtained by BET analysis of water isotherms. Inspection of water configurations shows that the water molecules do not cover the whole surface as already observed for xenon in the same vycor numerical sample [START_REF] Pellenq | [END_REF] . As a matter of fact, the molecules form disconnected patches of different sizes, living some regions completely dehydrated (not so hydrophilic sites). The BET surface is a thermodynamic surface in the sense that it is a determination of the area of substrate presenting energetically favorable sites. The previous result has then to be interpreted in terms of heterogeneity of the silica surface, with hydrophilic zones measured by the BET surface and almost as much hydrophobic regions not taken into account in the BET surface. This heterogeneity was previously noticed while analyzing the isosteric heat of adsorption. It seems then that the BET method for specific surface determination fails to be applicable to our water/silica system. It is proposed that this heterogeneity might explain the relative underestimation of the amount of water that should be adsorbed compared to experiments. 47 The presence of hydrophobic sites also induces a delay in capillary condensation which ends slightly above the saturating vapor pressure. This issue will be reconsidered while analyzing the water monolayer structure.

STRUCTURE ANALYSIS

The structure of water adsorbed on vycor has been analyzed in terms of site-site radial distribution functions. The data have been averaged in the Grand Canonical Ensemble, which is expected to give results identical to micro-canonical ensemble Molecular Dynamics. The site-site radial distribution function has been carefully normalized by taking into account the excluded volume effect. 16 In this structure analysis, we only consider the adsorbed film on the pore surface. Previous work already focused on the filling effect, showing strong structure distortion compared to bulk only for monolayer coverage. Actually, the structure of water is expected to be very close to bulk for multilayer coverage. In this analysis, the adsorbed film is supposed to be completed at 15 µmol/m 2 , which corresponds, by visual inspection to a nearly uniform coverage of the internal surface of vycor. This corresponds to a relative pressure of about 0.6. However, the higher the temperature, the more diffuse the layer, which is actually not so well defined at 500 K.

The oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen radial distribution functions are given for three different temperatures in figures 6, 7 and 8. Inspection shows that the O-O first peak is situated around 2.8 Å, as for bulk water. The position of the maximum is slightly shifted to larger values as the temperature increases, as expected. The larger the temperature the broader the peak, due to thermal agitation. The second peak is situated around 5 Å, slightly displaced compared to the bulk (4.5 Å). It corresponds to the formation of chains of water molecules each having two neighbors. Such a displacement was first observed by Brovchenko et. al 53,54 and also reported by Gallo et. al. 11 . Furthermore it has been shown 54 that the intensity of this displaced second peak depends on the strength of the water-substrate interaction. The more intense the interaction, the stronger the peak. This suggest that the second peak at 5 Å is related to the quasi bi-dimensional character of the adsorbed film. For temperatures above 400 K, this second peak disappears. The O-H correlation function also shows well defined peaks at positions close to bulk values. By contrast, peak intensities are lower for confined water. As for O-O correlation function, the peaks decrease, broaden and shift to larger distances as the temperature increases. The third peak disappears above 400 K due to thermal agitation. The same remarks hold true for the H-H radial distribution function.

In a neutron scattering experiment, the total site-site correlation function is a weighted sum of the partial atom-atom correlation functions (according to the atomic cross sections). Since the hydrogen cross-section is unfavorable, the hydrogen species are generally replaced by deuterium in real experiment, and in this case the total correlation function is written as: 55

) ( 090 . 0 ) ( 421 . 0 ) ( 489 . 0 ) ( r g r g r g r l D D D O O O - - - + + = (5)
This quantity is shown in figure 9 (see figure 13) show, for low temperature, a broad peak around 2.6 Å and a small one around 4.4 Å. These peaks are reminiscent of the first and third peaks in bulk water positioned at 2.4 and 4.6 Å respectively. It is expected that the second peak at 3.8 Å in bulk water be absent because the surface hydroxyls contain only one hydrogen. These data then infers that the H-H average distance is larger than in bulk. In conclusion, it seems that the hydrogen bond is elongated, and might not meet the energetic criterium stating that a hydrogen bond has a characteristic energy around 10 kJ/mol (analysis underway). These strong structural modifications might be understood by considering that the electric field in the vicinity of the surface hydroxyls differs significantly from that of a water molecule. For instance the dipole moment of the hydroxyls is directed along the O-H direction, which is not the case for water.

The partial charges carried by the oxygen and hydrogen species are not identical (see Table 1). Furthermore, the silicon tetrahedron gives rise to quadrupole moment in the vicinity of the silica surface. The figure 14 shows a molecular configuration of water adsorbed on the silica surface. Only two tetrahedral SiO4 are represented for clarity, and two adsorbed molecules. Inspection shows that the water molecules and surface hydroxyls are not inserted in a tetrahedral structure as mentioned before.

CONCLUSION

A thermodynamic and structure analysis of water / silica system at the mesoporous scale has been carried out in the Grand Canonical Ensemble using the Grand Canonical Monte Carlo simulation. The numerical mesoporous glass has been obtained by off-lattice reconstruction method. It reproduces quite well the geometric complexity of real vycor, like its specific surface area, pore size, chord distribution and correlation peak in small angle neutron scattering spectra. Its surface chemistry is made realistic by hydroxylating the surface. The

  figures 11 and 12 are strongly distorted compared to bulk. As the temperature decreases, a main peak appears around 2.8 Å instead of two in bulk water. It means that the dissymmetry of the distances between on one hand the oxygen of one molecule and hydrogen of the other

  adsorbed water is described by the SPC model, and water/silica interactions are calculated in the framework of the PN-TrAZ model. The thermodynamic quantities like water adsorption isotherms and isosteric heat of adsorption have been calculated for different temperatures. The model has proven to reproduce quite well the first steps of water adsorption. Since the

FIGURE 6 .

 6 FIGURE 6. Water oxygen-oxygen pair correlation function g(r): at 300K (circles), 400K (squares) and 500K (triangles).

FIGURE 7 .

 7 FIGURE 7. Water oxygen-hydrogen pair correlation function g(r): at 300K (circles), 400K (squares) and 500K (triangles).

FIGURE 8 .

 8 FIGURE 8. Water hydrogen-hydrogen pair correlation function g(r): at 300K (circles), 400K (squares) and 500K (triangles).

FIGURE 9 .

 9 FIGURE 9. Total water-water l(r) (see text).

FIGURE 10 .

 10 FIGURE 10. Water oxygen -vycor oxygen pair correlation function g(r): at 300K (circles), 400K (squares) and 500K (triangles).

FIGURE 11 .

 11 FIGURE 11. Water hydrogen -vycor oxygen pair correlation function g(r): at 300K (circles), 400K (squares) and 500K (triangles).

FIGURE 13 .

 13 FIGURE 13. Water hydrogen-vycor hydrogen pair correlation function g(r): at 300K (circles), 400K (squares) and 500K (triangles).

FIGURE 14 .

 14 FIGURE 14. Molecular configuration showing two water molecules adsorbed on a portion of silica surface. Yellow: Si, red: oxygen, grey or white: hydrogen.

Table 1 :

 1 PN-TrAZ parameters for water and vycor species

	a 0 =0.529177Å, 1E h =3.1578 10 5 K

Table 2 :

 2 Dispersion and repulsion parameters obtained in the framework of the PN-TrAZ model with Bohm and Ahlrichs combination rules for repulsive interactions.

Table 3 :

 3 SPC water thermodynamic quantities at different temperatures: P*: bulk saturating vapor pressure[START_REF] Errington | [END_REF]46 , L: bulk latent heat[START_REF] Errington | [END_REF]46 , qst: isosteric heat of adsorption on 3D representation of the numerical sample of Vycor used in this study. The porosity is represented in gray. The mean pore size is 36 Å. FIGURE 2. SPC model. FIGURE 3. Adsorption isotherms of water in Vycor at various temperatures.FIGURE 4. Isosteric heat of adsorption at various temperatures.FIGURE 5. BET plots for water adsorbed in Vycor at various temperatures.

	Vycor, n 0 : monolayer capacity, C parameter, and e-L (see text).
		300 K	350 K	400 K	500 K
	P*	4.8e-2	0.52	3	33
	L (kJ/mol)	45	41.4	37.8	30.6
	qst (kJ/mol)	80	76	70	65

geometric surface is known, the BET analysis allows to calculate the monolayer capacity, or the inverse which is the mean surface occupied by a single water molecule of the monolayer.

The value obtained (22 Å 2 ) is significantly larger than the value deduced from bulk water density (10.5 Å 2 ). It is qualitatively in agreement with experimental observations. On the other hand, the structure of the adsorbed film prior to capillary condensation has been studied in the same statistical ensemble. It shows a strong distortion of the pair correlation functions compared to bulk water due to the influence of the silica surface. The hydrogen-bond is significantly elongated, and angles distorted. In a further study the surface hydroxyls will be made free to rotate, which was not the case in this work for the sake of simplicity. However it is thought that this will probably greatly influence the monolayer structure, which might reaccommodate some bulk water features.

ACKNOWLEDGMENTS:


Ultra-broadband contactless imaging power meter

Knowledge 1

of the spatial and temporal distribution of heat flux is of great interest for the quantification of heat sources. In this work, we describe the development of a new ultra-broadband contactless imaging power meter based on electromagnetic to infrared technology. This new sensor and the mathematical processing of images enable the reconstruction of both spatial and amplitude distributions through a wide spectral range of sources. The full modeling of the thermoconverter based on 3D formalism of thermal quadrupoles is presented first before deriving a reduced model more suitable for quick and robust inverse processing. The inverse method makes it possible to simultaneously identify the heat losses and the spatial and temporal source distribution for the first time, to the best of our knowledge. Finally, measurements of multispectral sources are presented and discussed, with an emphasis on the spatial and temporal resolution, accuracy and capabilities of the power meter.

Nomenclature

There are two main types of flux sensors: thermal and photonic sensors. Photonic sensors are based on the photoelectric effect and are predominantly used at higher frequencies in the visible, ultraviolet, and X-ray spectral ranges [START_REF] Flores-Fuentes | Comparison between different types of sensors used in the real operational environment based on optical scanning system[END_REF]. In contrast, thermal detectors are mainly used at longer wavelengths [START_REF] Pradere | Photothermal converters for quantitative 2D and 3D real-time terahertz imaging[END_REF][START_REF] Müller | Novel detectors for traceable THz power measurements[END_REF].

In this study, the choice of hardware is based on the electromagnetic infrared (EMIR) technology developed by Balageas et al. [START_REF] Balageas | Characterization of electromagnetic fields using a lock-in infrared thermographic system[END_REF]. In their work, Balageas et al. reported a hyperspectral sensor using a very thin and homogeneous carbon film called a hyperspectral thermoconverter. This device has the ability to absorb radiation over a very wide spectral range (from visible to radio waves) with different sensitivities (depending on the wavelength) [START_REF] Pradere | Photothermal converters for quantitative 2D and 3D real-time terahertz imaging[END_REF][START_REF] Romano | Broadband sub-terahertz camera based on photothermal conversion and IR thermography[END_REF]. Electromagnetic waves are absorbed by the carbon film and converted into heat, which is transported by conduction in the thin film before being re-emitted as infrared (IR) light. Thus, this carbon thermoconverter enables imaging of ultra-broadband heat flux when used in combination with an IR camera. However, before obtaining a quantitative image of the heat flux or the thermal source detected by the sensor, thorough mathematical processing based on thermal inverse methods is required.

Depending on the nature of the heat transfer, i.e., conduction, convection, or radiation, there are several methods for estimating the heat source [START_REF] Gidik | Development of a radiative heat fluxmeter with a textile substrate[END_REF][START_REF] Jim | Estimating heat flux transmission of vertical greenery ecosystem[END_REF][START_REF] Yang | A temperature prediction-correction method for estimating surface soil heat flux from soil temperature and moisture data[END_REF]. In the present study, the reconstruction of the heat fields is essentially based on the knowledge of conductive transfer. Garderein et al. [START_REF] Gardarein | Miniaturized heat flux sensor for high enthalpy plasma flow characterization[END_REF] developed a point sensor system based on a thermocouple to estimate local fluxes by analytical inverse thermal methods. Another study from Zeribi et al. [START_REF] Zribi | Design, fabrication and characterization of thin film resistances for heat flux sensing application[END_REF] reported the fabrication of a 2D non-imaging heat flux sensor based on the spatial temperature gradient method. Image reconstruction was addressed by Groz et al. [START_REF] Groz | Bayesian inference for 3D volumetric heat sources reconstruction from surfacic IR imaging[END_REF][START_REF] Groz | Three-dimensional reconstruction of thermal volumetric sources from surface temperature fields measured by infrared thermography[END_REF], who reported a method to reconstruct deep heat sources using analytical models and two inversion methods (statistical and deconvolution by Toeplitz). Another 3D reconstruction method was recently developed by Burgholzer et al. [START_REF] Burgholzer | Three-dimensional thermographic imaging using a virtual wave concept[END_REF][START_REF] Waters | 3D reconstruction of tilted cracks using infrared thermography and the virtual wave concept[END_REF]. This method combines IR thermography and the concept of virtual waves. However, in these studies, the quantitative goal of estimating the heat flux amplitude was not achieved because the location of the sources is assumed to be unknown.

Finally, Nortershauser et al. [START_REF] Nortershauser | Resolution of a three-dimensional unsteady inverse problem by sequential method using parameter reduction and infrared thermography measurements[END_REF][START_REF] Nortershauser | Estimation of moving heat sources with a three-dimensional unsteady inverse method[END_REF] addressed the development of a source reconstruction method (spatial distribution and flux amplitude) using numerical models for inversion. The constraint of this method is mainly related to the calculation time.

However, the authors proposed the use of a statistical estimator in the cosine transformed space to reduce the calculation time.

From all of these seminal works, it was demonstrated that both The heat absorption depends on the wavelength. For example, 111 the absorbance is 100% in the NIR to far-IR range and only 112 61% in the millimeter range [START_REF] Romano | Broadband sub-terahertz camera based on photothermal conversion and IR thermography[END_REF]. Heating the film increases its

                           1 a r ∂ T r (x ,y ,z,t) ∂t -1 λ r =2 Y r =2 (x , y , z, t) -∇ 2 T r (x , y , z, t) = 0, r = 1, 2, 3, -λ r ∂ T r (x ,y ,z,t) ∂ x x =± L x 2 = 0, -λ r ∂ T r (x ,y ,z,t) ∂ y y =± L y 2 = 0, λ r =1 ∂ T r =1 (x ,y ,z,t) ∂z z= j = -λ r =2 ∂ T r =2 (x ,y ,z,t) ∂z z= j = -h c T r =2 (x , y , z = j , t), -λ r =2 ∂ T r =2 (x ,y ,z,t) ∂z z=k = λ r =3 ∂ T r =3 (x ,y ,z,t) ∂z z=k = h c T r =2 (x , y , z = k, t), T r =1 (x , y , z = j , t) = T r =2 (x , y , z = j , t), T r =2 (x , y , z = k, t) = T r =3 (x , y , z = k, t), T r (x , y , z, t = 0) = 0, ( 1
)
where r is the layer number, i.e., 1 for air, 2 for the thermocon-146 verter, and 3 for air, as depicted in Fig. 2.

147

B. Assumptions and Solution 148

In this part, the solution of the three-layer system is described. cosine Laplace spaces. The equivalent electrical network of the system is shown in Fig. 3(a), in which each layer is composed of three impedances. In this case, it is considered that the contact between the layers is perfect (no contact thermal resistance between the layers), and the internal source is uniformly distributed through the entire thickness of the thermoconverter.

To solve the three-layer system, the following assumptions are made:

• Adiabatic boundary conditions on each side of layers • Semi-infinite medium configuration: the purpose of choosing a three-layer system as a model is to take into account the influences of adjacent media (mainly air) on the thermoconverter, so the impedances Z 1 3 and Z 1 2 of layer 1 and the impedances Z 3 1 and Z 3 3 of layer 3 are replaced by two semi-infinite impedances, Z 1 ∞ and Z 3 ∞ , respectively [Fig. 3(c)].

• Thermally thin body: the thermoconverter was previously described as a thermally thin body (no temperature gradient in the thickness). This leads to the suppression of the impedances Z th 1 and Z th 2 , and only the impedance Z th 3 is necessary to evaluate the temperature of the thermoconverter [Fig. 3(d)].

Figure 3 shows the reduction steps according to the assumptions made, as well as the equivalent electrical networks obtained at each step. Based on the electric current conservation law at the node of the thermoconverter layer in Fig. 3(d), the temperature expression of the thermoconverter in the Laplace cosine transformed space is written as follows:

θ Th = 1 Z 1 eq + 1 Z Th 3 + 1 Z 3 eq -1 × Y , (2)
with

Z 1,3 eq = h + 1 Z 1,3 ∞ -1
, where Z 1,3 ∞ = (λγ ) -1 and

γ = p a + α 2 n + β 2 m . ( 3
)
The capacitive impedance of the thermoconverter is written as follows:

Y (α n , β m , p) = +∞ 0 L x 0 L y 0 Y(x , y , t) × exp(-pt) cos(α n x ) cos(β m y )dtdx dy . (5)
Finally, the temperature of the thermoconverter in the 192 Laplace cosine transformed space is obtained as

193 θ Th (α n , β m , p) = +∞ 0 L x 0 L y 0 T Th (x , y , t) × exp(-pt) cos(α n x ) cos(β m y )dtdx dy . (6)
Equation ( 2) can be rewritten simply as follows:

194 θ Th (α n , β m , p) = H(α n , β m , p) × Y (α n , β m , p), H(α n , β m , p) = 1 Z 1 eq + 1 Z Th 3 + 1 Z 3 eq -1 . ( 7
)
To compute the spatial temperature field, one inverse 195 Laplace in time [START_REF] Toutain | Numerical inversion of Laplace transform for time resolved thermal characterization experiment[END_REF] and two inverse cosine transformations 196 are necessary.

197 Table 1. Thermophysical Properties of the Different Layers Given in Ref. [START_REF] White | Series in Mechanical Engineering[END_REF] for Air Air Thermoconverter

Thickness (m) ∞ 37 ± 1 × 10 -6 λ (W m -1 K -1 ) 0.
C p = 1572 ± 8 J • K -1 • kg -1 .
The thermophysical properties of air were 5 obtained from the literature [START_REF] White | Series in Mechanical Engineering[END_REF] (see Table 1).

C. Validation of the Assumptions

Two studies have been carried out, as described in this section, to validate the assumptions made previously. The first study consists of a comparison between the temperature fields of the thermoconverter calculated by means of the four different models described in Fig. 3. The spatial shape of the source used is concentric circles, and a Dirac pulse δ(t) was used for temporal excitation. The value of the heat loss was set to must be equal to 100% at all times (conservation of energy). In Fig. 5(a), we can see that at short times (t < 2 s) the influence of the air layers appears very early and is important enough to be neglected (5% at t = 0.17 s, 10% at t = 0.8 s). Moreover, the heat flux dissipated by convective heat loss intervenes later than the air (5% at t = 1 s, 10% at t = 2 s) to then reach a more important rate at long times (35% at t = 10 s). In conclusion, it should be noted that the dissipation of the flux in the semi-infinite impedances of the air occurs rapidly and is quite Based on Eq. ( 7), the output temperature can simply be writ-253 ten as a space-time convolution product of the source and the 254 impulse response of the source point: 255 T Th (x , y , t) = Y(x , y , t) ĥ(x , y , t), [START_REF] Balageas | Characterization of electromagnetic fields using a lock-in infrared thermographic system[END_REF] where ĥ(x , y , t) is the impulse response in real space-time.

Y 0 × F(α n , β m ) = θTh (α n , β m , t) × Ĥ (α n , β m , t) -1 . ( 10
)
Inverse thermal problems are known to be ill-posed problems 268 [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Differential Equations[END_REF]. This is essentially due to the condition of instability of

Y 0 × F(α n , β m ) = θTh (α n , β m , t) × Ĥ (α n , β m , t) | Ĥ (α n , β m , t)| 2 + µ| D(α n , β m )| 2 , ( 11
)
where D is a derivation matrix [START_REF] Gray | Toeplitz and circulant matrices: a review[END_REF] in the cosine transformed space, and µ is the regularization coefficient [START_REF] Picard | Sur un theoreme general relatif aux equations integrales 613 de premiere espece et sur quelques problemes de physique mathe-614 matique[END_REF]. In the end, to retrieve the spatial distribution of the source, two inverse cosine transformations are necessary.

B. Heat Losses Estimation

As shown in Section 3.C, the heat losses by convection need to be measured along with the spatial distribution of the source in order to adapt to the external environment. Thus, based on the fact that the spatial average of a 3D temperature field leads to the one-dimensional (1D) temperature field [START_REF] Maillet | Thermal Quadrupoles, Solving the Heat Equation Through Integral Transforms[END_REF] and the fact that the relaxation of the temperature field due to a Heaviside excitation corresponds to a response to an amplified Dirac excitation, a method of heat loss estimation based on a linear least-squares minimization between the normalized spatial average of the relaxation experimental temperature field (see Fig. 6) and the normalized spatial average of temperature field obtained from the model Eq. ( 2) (for a Dirac time excitation) is described in this section.

Assuming that the measured temperature field of the thermoconverter is T mes (x , y , t), the spatial average of this temperature field is calculated for each time step and normalized by its maximum to arrive at T * mes (t).

Then, a variable change over time is applied, t * = tt 0 , where t 0 corresponds to the time of the beginning of relaxation Finally, the minimization is achieved using the Nelder-Mead simplex method [START_REF] Nelder | A simplex method for function minimiza-616 tion[END_REF],

Fig. 6. Square-shaped pulse time excitation to be used for the estimation of the heat losses as well as the heat source.

Table 2. Thermophysical Properties of the Foam [START_REF] Bondareva | The influ-618 ence of external temperature and convective heat exchange with an 619 environment on heat transfer inside phase change material embed-620 ded brick[END_REF][START_REF] Placido | Thermal properties 622 predictive model for insulating foams[END_REF] Thickness (m) λ (W m -1 K -1 ) ρC p (J K -1 m -3 ) a (m To simultaneously and continuously estimate the convective heat loss and the excitation flux of the source, a square-shaped pulse is used for temporal excitation of the internal source.

Figure 6 shows the square-shaped pulse time excitation as well as the temperature profile of the thermoconverter in the presence of the resulting heat losses. The choice of square-shaped pulse time excitation is justified by the fact that there is a constant level where the internal source is switched on (temperature rise of the thermoconverter), which is used to identify the source, and another constant level where the internal source is switched off (relaxation of the thermoconverter), which is used to estimate the heat loss. The estimation of the heat loss makes it possible to reintroduce this loss into the model described in Section 3.B and thus makes it complete.

RESULTS AND DISCUSSION

A. Contact Validation by Joule Effect

One of the most robust ways to control the power dissipated by a source is to use an electrical resistance heated by the Joule effect. This is used to validate the method of reconstruction of the source. To do this, a complex-shaped resistance with an internal ohmic resistance of 73 is used. The setup used is described in • K -1 . The estimated value of this coefficient is coherent, given that one side of the thermoconverter is insulated by the foam, and the average increase in the temperature is low (on the order of 0.4 K at t = 1 s). where S pixel is the area of one pixel with a value of 8 × 10 -8 m 2 .

The estimation represents a relative error of 0.45%.

B. Contactless Application for Different Optical Sources

In this section, the use of our appliance as a power meter is demonstrated.

Near-IR Laser Source

In this first application, a NIR laser diode (λ = 980 nm) with a power of P = 280 mW is used as an excitation source. The time of the square-shaped pulse is 40 ms. A dichroic mirror is placed between the IR camera and the thermoconverter, as shown in Fig. 9, and enabled to reach the temperature of the thermoconverter at the front side. As the thermoconverter is a thermally thin body, the temperature at the front side is equal to the temperature at the back side (no temperature gradient along the thermoconverter thickness). The thermoconverter is placed at a distance of 40 cm from the IR camera. • K -1 . This value is very large compared to that estimated for the Joule heated resistance. This is due in small part to the fact that the thermoconverter is not insulated and in large part to the fact that the average temperature of the thermoconverter is very high (4 K at t = 33 ms). The power estimated using Eq. ( 13) was found to be 285.3 ± 3.7 mW with a relative error of 1.9%.

IR Blackbody Source 401

In this second application case, the source used comes from an • K -1 . This value is lower than that in the NIR laser case because the average temperature of the thermoconverter is also lower (2.5 K at t = 1 s).

The power estimated using Eq. ( 13) is 1.357 ± 0.015 W. To

Gigahertz Source

In this third application case, a Terasense gigahertz source ( f = 100 Ghz, λ = 3 mm) with a power of P = 400 mW is used, and the absorbance (A) of the thermoconverter in this wave is equal to 61% [START_REF] Romano | Broadband sub-terahertz camera based on photothermal conversion and IR thermography[END_REF]. The thermoconverter is placed at a distance of 5 cm from the source. The gigahertz source is synchronized with a waveform generator that allows delivery of a square-shaped pulse of 1.25 s.

The estimated convective exchange coefficient is

h = 20.24 W • m -2 • K -1 .
This value is higher than that in the blackbody case and lower than that in the NIR laser case. This is related to the value of the average temperature of the thermoconverter, which is between those of the two cases (6 K at t = 1 s). Figures 11(b) and 11(c) show the reconstruction of the spatial distribution and power density of the source after inversion. In this case, the source power is estimated by modifying Eq. ( 13) to take into account the absorbance of the thermoconverter at the emission wavelength of the source as follows:

P estimated = Y(x , y ) × S pixel × 1 A . ( 14
)
The power estimated using Eq. ( 14) is 404.8 ± 5 mW with a relative error equal to 1.19%.

Radio-Frequency Source

In this fourth application case, an ultra-high-frequency (UHF) • K -1 . This value is the lowest of the convective exchange coefficient values estimated for the optical sources (except for the joule heat resistance). This is because the average temperature of the thermoconverter is lower than that in the previous cases [34] (1.3 K at t = 1 s). The power of the source estimated using Eq. ( 14) is 6.5 ± 0.17 × 1 A W. This

j

  Temperature in real space, K θ Temperature in the cosine transformed space θ Temperature in cosine and Laplace transformed space h c Convective exchange coefficient, W m -2 K -1 ĥ Source point impulse response in real space H Source point impulse response in cosine and Laplace transformed space Ĥ Source point impulse response in cosine transformed space Ĥ Source point Heaviside response in cosine transformed space Y Internal source in real space, W m -3 Y Internal source in cosine and Laplace transformed space δ Dirac function Heaviside function λ Thermal conductivity, W m -1 K -1 a Thermal diffusivity, m 2 s -1 ρ Mass density, kg m -3 C p Specific heat, J K -1 kg -1 φ 0 , φ e Excitation flux in boundary condition (respectively, input and output fluxes), W m -2 e Thickness of the thermoconverter, m (Table continued) L Lateral dimension of the thermoconverter, m t Time, s x , y , z Spatial coordinates Z i Thermal impedance in Laplace and cosine transform space, i layer number, j Thermal impedance number considerable amount of attention in the research community in recent years. Whereas many works have 25 focused on imaging sensors for visible light [1-4], few tools exist 26 for other wavelength ranges, particularly the millimeter to meter 27 range [5,6]. Knowledge of the heat flux generated by a multi-28 spectral source is of prime interest in numerous experimental 29 setups, i.e., in building sciences, aeronautic industries, optical 30 applications, and heat transfer. A quantitive image of radiative 31 heat flux (or thermal source) would enable the performance 32 of heat balance to serve as input data for models or to control 33 online processes. In this context, the development of a new 34 sensor implies, first, developing the hardware (the component sensitive to the heat excitation) and, second, determining the mathematical modeling of the hardware response to obtain quantitative data.

Fig. 1 .

 1 Fig. 1. A multispectral source is used to optically illuminate

  110

Fig. 1 .

 1 Fig. 1. Schema of the experimental setup.

  Thermal Model of the Thermoconverter 126 The thermoconverter can be considered a thin homogeneous 127 medium (a thermally thin body due to its small thickness, 128 L = 5 cm and e = 37 µm) subjected to convective/radiative 129 heat loss and containing an internal heat source. Due to the very 130 large ratio between the lateral dimensions and the thickness 131 (L e ), the temperature gradient across the thickness of the 132 thermoconverter is neglected. Thus, it is advisable to consider 133 the surrounding environment [21]. The influence of the air 134 around the thermoconverter is not negligible and considerably 135 affects the temperature of the thermoconverter. This can be 136 explained by the presence of effusive exchanges between the 137 thermoconverter (thermally thin body) and the surrounding 138 environment. Therefore, the analytical model that accurately 139 describes the heat transfer in the thermoconverter must include 140 the surrounding layers. In Fig. 2, there are three layers: air-141 thermoconverter-air. The evolution of the temperature fields 142 in the system described in Fig. 2 is determined by the three-layer 143 3D heat conduction equation with an internal volume source 144 located in the thermoconverter layer: 145

149Fig. 2 .

 2 Fig. 2. Complete thermal model of the thermoconverter (i, l, physical boundary of the air layers; j, k, plans of the interfaces between the thermoconverter and the air layers).

1 and 3 :Fig. 2 (

 32 Fig. 2 (absence of φ 0 and φ e ).

Fig. 3 .

 3 Fig. 3. Reduction steps of the thermoconverter model: (a) complete model, (b) adiabatic boundary conditions model, (c) semi-infinite configuration model, and (d) thermally thin body model.

  026 1.414 ± 0.041 ρC p (J K -1 m -3 ) 1313 2.83 ± 0.03 × 10 6 a (m 2 s -1 ) 1.98 × 10 -5 5.0 ± 0.1 × 10 -7 All of the parameters of the thermoconverter needed for the construction of the model have been estimated. The thickness measured by a micrometer was 37 ± 1 µm. The thermal diffusivity was estimated by the flying spot technique [25] to be a = 5.0 ± 0.1 × 10 -7 m 2 • s -1 . Using conventional thermal characterization techniques, the remaining thermophysical properties of the thermoconverter were measured: (i) the density was calculated by means of a helium pycnometer to be ρ = 1800 ± 10 kg • m -3 , and (ii) the specific heat was estimated by a Setaram 131 differential scanning calorimeter to be

Fig. 4 .

 4 Fig. 4. Comparison of the temperature determined by the four models: (a) diffusion of the source obtained by the complete model, (b) calculated error (in %) between the source diffusion obtained from the complete model and that obtained from the adiabatic boundary conditions model, (c) calculated error (in %) between the source diffusion obtained from the complete model and that obtained from the semi-infinite configuration model, (d) calculated error (in %) between the source diffusion obtained from the complete model and that obtained from the thermally thin body model, and (e) diffusion of one pixel versus time determined by the four models.

h c = 20 WFigure 5

 205 Figure 5 illustrates the flux balance.Figure 5(a) shows the rate of the transient heat flux (delivered by the internal source) dissipated by the volume of the thermoconverter, air, and convective heat loss.Figure 5(b) shows the sum of all of the fluxes, which

Fig. 5 .

 5 Fig. 5. Heat flux balance.

256, ( 9 )

 9 The source Y(x , y , t) can be decomposed into a product of a 257 spatial function F(x , y ), amplitude Y 0 , and a temporal func-258 tion (t). After passing through the space transformed cosine 259 base, only the temporal convolution remains. In the case of 260 a Heaviside-type temporal excitation and after applying the 261 Laplace transform on time, we have 262 θ Th (α n , β m , p) = Y 0 × F(α n , β m ) × 1 p × H(α n , β m , p) H where H represents the response of the source point to the 263 Heaviside temporal excitation in Laplace cosine transformed 264 space. By applying the Laplace inverse transform to Eq. (9), the 265 source can be estimated in the cosine transformed space via the 266 following relation:267

  269 the solution obtained by inversion. To remedy this condition, 270 the inversion is carried out by constructing a Wiener filter. This 271 filter is based on Tikhonov's regularization method [28], which 272 is used and applied in the Cosine transformed space as follows: 273

  ( T * mes (t 0 ) = 1); this new time base t * is used to calculate the numerical temperature field of the thermoconverter (for a Dirac time excitation) from the model described in Section 3.B for different values of h c . The temperature field obtained is spatially averaged and then normalized by its maximum to obtain the following: T * model (t * ).

Fig. 7 (

 7 Fig. 7(a). An electrical current generator supplied the resistance during 1.5 s. The voltage is verified at the edge of the resistance by a voltmeter to be U = 4.93 V. This theoretically corresponds to the dissipation of a power equal to P = 333 mW. The thermoconverter is attached on the resistor, which is insulated by foam (polyurethane foam). Finally, an IR camera is used for the acquisitions. In this case, only the parameters of the model to be used for the inversion are modified. The proposed model allows changing the parameters easily to match the real experimental configuration. One of the semi-infinite impedances of the air is then replaced by a semi-infinite impedance of the insulated foam, and the heat loss by convection is neglected on this side.

Figure 7 (

 7 Figure 7(b) shows the normalized spatial average of the measured 3D temperature field. This field allows the identification of the t 0 and the construction of the new time base t * , which is used to estimate the heat loss.Figure 7(c) shows

Figure 7 (

 7 c) shows the result of the minimization described in Section 4.B, and the convective exchange coefficient is estimated to be 4.53 ± 0.3 W • m -2

Figures 8 (

 8 Figures 8(a) and 8(b) show the reconstruction of the spatial distribution and power density of the source after inversion.

Fig. 7 .Figure 8

 78 Fig. 7. (a) Experimental setup, (b) normalized spatial average of the 3D temperature field measured, and (c) comparison between the experimental data and the model using h c = 4.53 W • m -2• K -1 .

Fig. 8 .

 8 Fig. 8. (a) Image of the estimated source, (b) surface of the estimated source, and (c) estimated flux density for one pixel.

Figures 9 (

 9 Figures 9(b) and 9(c) show the reconstruction of the spatial distribution and energy density of the source after inversion. The diameter of the spot after the flux estimation is equal to 3 mm, which is in agreement with the optical considerations. The estimated convective exchange coefficient is h = 31 W • m -2• K -1 . This value is very large compared to

Fig. 9 .

 9 Fig. 9. (a) Experimental setup, (b) image of the estimated source, and (c) surface of the estimated source.

Fig. 10 .

 10 Fig. 10. (a) Experimental setup, (b) image of the estimated source, and (c) surface of the estimated source.

  Prisma instruments, whose emis-403 sivity is estimated to have an uncertainty of 0.5% to ε BB = 0.98. 404 The blackbody temperature is adjustable in the range [500 • C; 405 1200 • C]. In our case, the blackbody temperature is set to 406 T BB = 500 • C [see Fig. 10(a)]. A square-shaped pulse of 4.5 s is 407 applied using a chopper synchronized with the IR camera.

Figures 10 (

 10 Figures 10(b) and 10(c) show the reconstruction of the

  validate the value of the estimated flux, a short calculation based on the Stefan-Boltzman law and the notion of the form factor in thermal radiation to allow modeling the radiative exchanges between two black disks (ε ≈ 1) separated by a perfectly transparent medium is performed. This calculation predicts that the thermoconverter should theoretically receive a flux of P = 1.15 W and shows that the estimated flux is in the order of the magnitude of the expected power with a relative error equal to 15.25%.

  radio wave antenna source ( f = 500 Mhz, λ = 0.6 m) in an anechoic chamber is used. The thermoconverter used has a surface area of 1 m 2 , and it is placed at a distance of 1 cm from the source. The power delivered by the source and the absorbance of the thermoconverter at this wavelength are unknown, so the estimated flux is proportional to the absorbance of the thermoconverter. Figures 12(b) and 12(c) show the reconstruction of the spatial distribution and power density of the source after inversion. The estimated convective exchange coefficient is h = 10.14 W • m -2

Fig. 11 .

 11 Fig. 11. (a) Experimental setup, (b) image of the estimated source, and (c) surface of the estimated source.

Fig. 12 .

 12 Fig. 12. (a) Experimental setup, (b) image of the estimated source, and (c) surface of the estimated source.

Table 3 .

 3 Performances of the Flux Sensor

Joule Effect by Resistance NIR Laser Source IR Blackbody Source Gigahertz Source Radio-Frequency Source

	Spatial resolution (µm 2 )	283 × 283	180 × 180	290 × 290	280 × 280	1400 × 1400
	sensitivity (µW/m 2 )			1.815 × TI(µs) A
	pitch size, minimum sensitivity of 20 mK, and variable accord-

J.-W. Wu, W.-F. Sung, and H.-S. Chu, "Thermal conductivity of poly-

urethane foams," Int. J. Heat Mass Transfer 42, 2211-2217 (1999).

thank ADEME for the support of 516 this work through the IGAR project with Arcelor-Mittal.

underlying the results presented in this paper are 519 not publicly available at this time but may be obtained from the authors upon 520 reasonable request.

ing to the integration time), and, on the other hand, by the 479 photothermal effect within the thermoconverter (absorbance, 480 7 diffusion of the heat, and inversion of the heat transfer). The