A Chopard

Q Cassar

J Bou-Sleiman

J P Guillet

M Pan

J B Perraud

A Susset

Patrick Mounaix

J Bou Sleiman

Terahertz waves for contactless control and imaging in aeronautics industry

Keywords: Terahertz, Imaging, Multilayer painting, FMCW radar, Time domain spectroscopy

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Terahertz systems are demonstrating an increasing appeal through a wide panel of applicative fields thanks to their ability to provide contactless and non-destructive sensing solutions for in depth inspection. From security applications [START_REF] Cooper | Penetrating 3-D imaging at 4-and 25-m range using a submillimeterwave radar[END_REF] to biomedical sensing [START_REF] Cassar | Pilot study of freshly excised breast tissue response in the 300 -600 GHz range[END_REF][START_REF] Grootendorst | Use of a handheld terahertz pulsed imaging device to differentiate benign and malignant breast tissue[END_REF] through non-destructive testing (NDT) [START_REF] Cristofani | Nondestructive testing potential evaluation of a terahertz frequencymodulated continuous-wave imager for composite materials inspection[END_REF][START_REF] Friederich | Terahertz Radome Inspection[END_REF] or even art inspection [START_REF] Dandolo | Terahertz frequency modulated continuous wave imaging advanced data processing for art painting analysis[END_REF][START_REF] Guillet | Art Painting Diagnostic Before Restoration with Terahertz and Millimeter Waves[END_REF], specific integrations are demonstrating the multiple applicative capabilities of terahertz technologies.

The wide advances in this field in the recent years enabled the development of compact and cost-effective systems, driven by the increasing demand in dedicated complementary techniques applied to specific niche problematics [START_REF] Ellrich | Terahertz Quality Inspection for Automotive and Aviation Industries[END_REF]. Among those, complex materials integrity monitoring as well as penetrating imaging and sensing for NDT tasks is proving to be an ideal application fields for frequency modulated continuous-wave (FMCW) radars while an increasing demand in stratified structures inspection for the automotive and aerospace industries is being fulfilled with the implementations of terahertz time-domain spectroscopy (TDS) techniques.

Radar imaging Capabilities for NDT

The potential of FMCW sensing comes from its intrinsic benefits such as its high level of integration and compactness, its relatively low cost as well as the harmless character of the emitted radiations and their a priori high penetration capabilities in optically opaque materials. These arguments make FMCW a first-grade choice to complete the panel of tools allowing contactless inspection. Additionally, FMCW radars intrinsically provide information on the wave propagation delay and therefore straightly enables depth sensing. Hence, direct 3D imaging can be assessed [START_REF] Perraud | Shape-from-focus for real-time terahertz 3D imaging[END_REF][10], representing, along with the previously stated specificities, a non-negligible benefit in comparison with X-Rays imaging, low frequency continuous wave systems or thermographic or ultra-sound inspection methods. The detection and localization of damages, cracks, impact, delamination or inclusions are then achievable for materials that are commonly used in the aeronautics industry, such as glass fiber-based composite, honeycomb structures, polymers or ceramics either during fabrication or maintenance processes.

FMCW radar architecture and principle

FMCW radar have found several applicative purposes, from automotive anti-collision systems as precursors, to high-end NDT tools and have been developed on various technologies, mainly Si-based fully integrated devices [START_REF] Pauli | Miniaturized Millimeter-Wave Radar Sensor for High-Accuracy Applications[END_REF] or III/V high frequency precision radars [START_REF] Wong | Computational Image Enhancement for Frequency Modulated Continuous Wave (FMCW) THz Image[END_REF]. Despite the inherent benefits and drawbacks of each technology, their implementations nevertheless rely on similar architectures, described in Figure 1 (a). The first stage of the system, working at low frequency, features a voltage control oscillator (VCO), eventually combined to a phase locked loop (PLL) for stabilization purposes, in order to generate fast saw tooth-like frequency sweep. The multipliers chain is then introduced for up-conversion to reach the desired working frequency band, typically centred around 100 GHz, 150 GHz, 300 GHz or even up to 600 GHz. Nevertheless, the higher the frequency band, the lower the emitted power due to conversion losses. This high frequency sweep, considered as a reference signal, is then emitted towards the target (plain red Figure 1 (b)). By modulating the frequency in such a way, it is equivalent as putting a unique "time stamp" on the emitted wave at every instant.

The reflected sweep, delayed by the back and forth propagation through and into the target, (dotted orange Figure 1 (b)) is then collected and guided towards the detection mixer, thanks to a directional coupler, where a beating electric field is generated via the multiplication of this return sweep with the reference electromagnetic wave or one of its lower harmonics, due to their slight frequency difference. Collected by a data acquisition unit (DAQ), this beating signal displays a main frequency 𝑓 !"#$% , proportional to the propagation delay and so the distance to the object (see Equation ( 1)). From this formulation, the longitudinal resolution 𝛿 !"# is derived (see Equation ( 2)) and solely relies on the bandwidth of the radar transceiver.

𝑑 = 𝑐 $ 𝛥 % 2 𝑛 = 𝑐 $ 𝑓 &'("! 2 𝑛 . 𝑇 #)""* 𝐵𝑊 , (1)
𝛿 !"# = 𝑐 $ 2 𝑛 𝐵𝑊 , ( 2
)
where 𝑐 $ is the speed of light in vacuum, 𝛥 & , the back and forth propagation-induced time delay, 𝑛, the refracting index of the propagation media at the working frequency, 𝐵𝑊, the frequency bandwidth, 𝑓 !"#$% , the beating frequency and 𝑇 '($$) , the sweep period. For such radar unit architectures, the imaging capabilities are not intrinsic to those implementations due to the emission of a diverging beam from the antenna. The imaging process [START_REF] Cristofani | Nondestructive testing potential evaluation of a terahertz frequencymodulated continuous-wave imager for composite materials inspection[END_REF], [START_REF] Ding | THz 3-D Image Formation Using SAR Techniques: Simulation,ProcessingandExperimentalResults[END_REF] can nevertheless be performed through synthetic aperture radar reconstruction or more classically thanks to a focused raster scanning implementation. Due to the complexity of the inspected objects for see through imaging problematics, the latter solution is chosen and a quasioptical coupling system for focalized beam shaping is implemented with 2 Teflon ® lenses, as displayed in Figure 2. Such an integration enables a simplified optimisation of the beam focus on the object of interest, beyond an eventual extra partially transparent cover, to reach the optimum returned power and resolution.

The Gaussian beam model, provides the direct relation between the beam waist diameter ω $ , hence the resolution limit, and the focal length, f', of the focusing lens, as follows:

ω $ = 4 𝜆 𝑓′ π D , (3)
with D the diameter of the lens, and 𝜆 the wavelength related to the emitting source.

Therefore, for such a transceiver implementation, the resolution relies on the frequency band for the longitudinal resolution, as well as the optical configuration through the beam shaping system for the lateral resolution. At given working frequency, a trade-off between achievable lateral resolution and optimum imaging distance needs to be considered. Moreover, the adequate frequency band selection thus represents one of the key optimization since higher frequency bands would ensure an improved lateral resolution (see Equation 3) and would allow larger achievable bandwidth, therefore, providing enhanced longitudinal resolution as well (see Equation 2). Nevertheless, a wide range of aeronautics related materials remains sufficiently transparent for millimeter wave radiation, but commonly displays a significant increase in absorption with frequency [START_REF] Perraud | Terahertz imaging and tomography as efficient instruments for testing polymer additive manufacturing objects[END_REF], therefore resulting as well in a trade-off between material penetrability and targeted resolution depending on the defect size.

3D Imaging results toward see-through applications

A wide range of samples based on glass-fiber and reinforced plastic have been inspected using a SynView imaging system in reflection or transmission mode and have been evaluated using 100 GHz and 300 GHz radar transceivers, as displayed in Figure 3 for health monitoring applications. Those images, are resulting from processed 2D radar scans, providing 3D inspection data thanks to the depth sensing capability of the FMCW method, with data retrieval in the object's best imaging plan. The high measurement rate of FMCW radar units, typically in the range of few kHz, correlated with the millimeter achievable resolution, ensures an adequacy with the imaging cadencies expected for NDT applications, with a typical measurement time of a few minutes for such sample dimensions. As expected, higher lateral and longitudinal resolutions are achieved with the 300 GHz scan while the lower measurement dynamic range of the systems as well as the reduced material penetrability represents a large drawback for such inspection tasks. Nevertheless, intentional foreign inclusions remain discernible for both measurements. On the contrary, carbon fiber-doped samples, displays high reflectivity over the terahertz range, preventing in-depth or see-through inspection (maximum of one carbon fibre layer can be screened) but allowing topographic inspection or coatings analysis. Focusing on the see-through imaging problematics nevertheless brings to light few challenges. Targeting in-situ direct inspections while avoiding cover dismantlement, the reflection geometry represents the only adequate sensing implementation. Due to lowered material absorption levels at low frequencies, the 100 GHz band radar transceiver turned out to be the most suitable in order to minimize the losses when passing back and forth through the fiber-glass panel while enabling objects sensing and guaranteeing adequate resolution of a few millimetres.

As a result, Figure 4 (a) displays the first demonstration of this system, in such a complex configuration of see-through plane cover imaging. A 3D animated representation, resulting from the 2D scan, of multiple reflective samples behind a plane cover, is available through the attached supplementary materials. The goal is to evaluate if we can properly detect and identify a metallic reflective object, such as a wrench, a tube or the sample holders. All the ½ inch metallic post are also clearly imaged and resolved with a very good contrast. To push further the investigation, other objects that displays much lower reflectivity levels, such as a dry and moist sponges (Figure 4 (a) grey dotted box) held below this metallic wrench, are inserted, aiming to check the water content detection capability in case of crack or fissure inducing any fluid leakage. Clearly noticeable, a differentiation of the moisture level is achieved between the two sponges in presence and absence of water content.

Reminiscences of the panel geometry remains nevertheless discernible. We obtained useful images through the cover over distances up to 150 mm between the imaging area and the covering panel. Pushing further the tests, towards real case implementation, we demonstrated that the insertion of extra thermal insulation layers has no effect on the imaging capabilities thanks to their full transparency in the low frequency terahertz range (confirmed by terahertz spectroscopy). Beside simple object sensing and identification, Figure 4 (c) depicts the capability of see-through radar imaging for detection of defects and alterations, such as cracks or broken components, through a cover. For this implementation, a zoom on part of a 19 mm diameter aluminium tube (delimited by grey horizontal lines) displaying two controlled breaks (c1 and c2 in the grey dotted square box).

Finally, besides the many parameters impacting the sensing capabilities of radar units, such as the investigated feature's depth and size and their optical proprieties, the visualization of defects like foreign material, object, alterations or moisture contamination is made possible through a glass-fiber plane cover. A monitoring of the spatial location of objects of interests can moreover be performed thanks to the inherent depth sensing capabilities of radar imaging techniques.

The embedded phase information provided by the FMCW sensing process turns out to be an adequate tool for 3D imaging and in-depth sensing for a large scope of applicative cases of failure. Among those, the inspection of multi-layered polymers and foams structures for thickness inspection illustrates those capabilities [START_REF] Schreiner | All-Electronic High-Resolution Terahertz Thickness Measurements[END_REF]. With similar motives [START_REF] Roman | Terahertz dielectric characterisation of polymethacrylimide rigid foam: The perfect sheer plate?[END_REF] and capabilities, but targeting scaled down structures, terahertz thickness inspection have been successfully implemented for micrometric stacked paint structures in the aeronautics and automotive industry with the use of time domain spectroscopy tools

Thickness extraction and optical paths

Although the use of FMCW radars has proven to be effective for contactless and non-destructive in-depth measurements for the aviation industry, complementary technologies whose operating frequency range remains within the terahertz band can address additional issues encountered by this industrial sector. That is notably the case of terahertz time-domain spectroscopy, primarily introduced in the field to control the thickness of aircraft coatings. The sensitivity of pulsed terahertz radiations to individual coatings forming a stratified structure was reported for different study cases, as in automotive, art science and obviously, aeronautics [START_REF] Krimi | Highly accurate thickness measurement of multi-layered automotive paints using terahertz technology[END_REF][START_REF] Picot | Terahertz Paint Thickness Measurements : from lab to automotive and aeronautics industry[END_REF][START_REF] Van Mechelen | Stratified dispersive model for material characterization using terahertz time-domain spectroscopy[END_REF][START_REF] Cassar | Characterization of Varnish Ageing and its Consequences on Terahertz Imagery: Demonstration on a Painting Presumed of the French Renaissance[END_REF]. Indeed, when reaching a dielectric interface, a portion of the incident frequency-dependent terahertz electric field is transmitted through the following layer, while another portion Inverse electromagnetic problems are commonly performed to extract, from a recorded terahertz electric field the individual layer dielectric properties and thicknesses. Usually, it consists in minimizing an objective function that denotes the discrepancies between the experimental waveform and the ones computed from a set of candidate parameters that are aimed to impact the shape of the signal as the complex refractive index and the thickness. One can, non-exhaustively, distinguish three main different ways to proceed: (i) incremental Newtonian methods; (ii) full-volume test; (iii) dynamic stochastic strategy. In the following, each method is briefly explained and their main advantages and inconveniences are underlined.

Incremental Newtonian methods

The incremental Newtonian approach [START_REF] Goodfellow | Deep learning[END_REF] consists in evaluating iteratively the curvature of the objective function based on the ratio of the gradient of this objective function over its Hessian. While such a process would make high progress initially, the speed of convergence decreases as it gets closer to the physical solution. Hence, it is suitable to set a threshold under which, the method ends.

Full-volume test

In comparison to the previous method, a full-volume test successively computes the objective function for all possible combinations of allowed-parameters. Therefore, the method is expected to converge to the physical solution. However, it is worth to note that the nominal number of calculations can be immense so that the process becomes cumbersome and would not, ultimately, be able to provide the physical parameters that compose the stack in a decent period of time.

Dynamic stochastic strategy

In comparison to the two previous methods, a stochastic process is aimed to randomly pick different combinations of parameters that are subsequently injected into the objective function. Two different approaches are possible to select the best parameters: (i) the test ends when a combination provides an objective function minimization that is greater than a predefined correlation threshold; (ii) the process terminates after a given number of shots and selects the combination that provides the smallest value of the objective function. As dealing with random selection processes, the convergence time of the first approach remains erratic and unpredictable. The later approach is accomplished within a fixed period of time but may not grasp the physical combination if the number of stochastic shots revealed to be too small. Note that in a similar fashion to the definition of convergence time for the first method, the number of required shots to find a valid combination cannot be predicted. However, in both cases the probability to pick the valid combination increases as the number of allowed parameters decreases.

Iterative tree algorithm (ITA)

The aforementioned procedures are then commonly implemented to extract the physical parameters of each individual layers that constitute a stack. However, those methods gather indistinguishably the contributions of all the possible propagation paths through the paint stack, thus requiring a large computation power. The number of calculated paths could be gigantic while a precise knowledge of the most contributing optical paths would awfully minimize the quantity of computed routes and clarify the involved propagation processes. Furthermore, the derivation of a simplified parametric transfer function would drastically reduce the required computation power to minimize the objective function of an inverse electromagnetic formulation. A new algorithm monitoring the pulse subdivisions at dielectric interfaces has been to this end, developed [START_REF] Cassar | Iterative Tree Algorithm to Evaluate Terahertz Signal Contribution of Specific Optical Paths within Multi-Layered Materials[END_REF]. The algorithm to be reported does not aim at evaluating the respective thicknesses of individual coatings but intends to establish a parametric transfer function, also referred as "reduced transfer function" in the following, which is derived from the analysis of individual optical route contribution to the responsive terahertz frequency-depend electric field. Ultimately, the so-called parametric transfer function can be used within the designed objective function to be minimized to extract, from one of the aforementioned property extraction methods, the thicknesses of coatings.

The iterative tree algorithm (ITA) classifies each portion of the signal as a function of the number of encountered dielectric interfaces instead of considering the signal as a whole. Additionally, it identifies the various optical paths involved for each contribution. The optical routes that contribute the most to the power are labelled and used to derive a parametric transfer function.

The central principle of the algorithm is to increment, at each encountered dielectric interface, a parameter k denoting pulse subdivisions (Figure 5). Each propagating radiation is thus the combination of reflection and transmission processes involving surrounding dielectric layers, and occurring at the previous algorithm step, i.e. k-1. One has then the following equivalent mathematical description of the scheme: [START_REF] Cristofani | Nondestructive testing potential evaluation of a terahertz frequencymodulated continuous-wave imager for composite materials inspection[END_REF] where and are the reflected and transmitted signals in frequency domain of the k-th iteration at each dielectric interface.

In order to pilot and to illustrate the iterative tree algorithm, an aeronautic stack sample on metallic substrate made of 4 layers of tens of microns has been extensively investigated in the THz domain and with conventional optical destructive methods for reference. The first step to achieve this investigation, was to carefully extract the dielectric properties of each individual painting material in isolation. To do so, the full-volume test reported in section 3.2 was followed. Each individual layer dielectric properties were then fitted to a single-pole Debye profile, as a consequence of a generalised thermal effect. The full stack was then reconstructed by means of equations ( 4) and its reconstruction iteratively processed is given in Figure 6.

Discrepancies between the computed waveform and the experimental one decreases when considering an increasing number of contributions until it reaches a threshold under which the error remains overall stable. The error metric, defined as the difference squared between the two waveforms, as a function of the number of iterations is given in Figure 7.

The non-monotonous nature of the dependence reflects that time features are in majority due to the combination of several physical processes. Despite additional iterations, the error function does not reach zero. The presence of a non-zero asymptote may be attributed to various roots: the non-perfect extraction of individual dielectric profiles, experimental limits and additionally due to the fact that time features corresponding to further iterations are out of the time window in which the comparison metric is defined.

One has to keep in mind that such a reconstruction actually considers the whole possible optical routes within the stack. However, one can anticipate that all of these routes are not significantly contributing to the reflected signal. Thanks to the mathematical design of the algorithm, it is possible to analyse and quantify the contribution of each optical path individually and consequently to sort them as a function of their respective contributions to the responsive waveform. Doing so, one can therefore calculate a parametric transfer function that only considers a selection of the most contributing routes. Some of these identified optical paths are shown in Figure 8. It is worth to note that the procedure initially dealt with 1,094 distinct paths between k = 1 and k = 17. However, the identification of the most contributing paths led to a shortlist of 33 routes. Globally, the time needed to reduce the number of optical paths to be considered per pulse subdivision varies from few seconds to few minutes, but remains dictated by solver configuration. Notably, the higher the iteration is, the longer the required time for optical routes minimization is.

, 1 1 1 1 1 1 1 1 ( ) ( ) ( 1 1, , 1 , 1 , 1 1, , 1 , 1 , 1 1, , 1 , 1 , 1 1
( ) i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i k k k r r r t k k k r r r t k k k t t )
- + + + + - - - - + + - - - - - - - - ì = + ï = + ï í = + ï ï = + î k r S k t S
Deeper explanation about selection rules can be found in [START_REF] Goodfellow | Deep learning[END_REF]. Ultimately, these 33 selected optical paths were used to reconstruct the previously simulated waveform (Figure 9). This new algorithm is used to understand a stratified material response by terahertz time domain characterization in reflection. The identification of the main optical paths responsible for the reflected signal led to a reduced transfer function that was used to reconstruct almost perfectly the original responsive waveform. The generated parametric transfer function can eventually be implemented in thickness determination algorithms for micrometric coating stack. This algorithm also carries interest to deepen theoretical understanding on the way a responsive terahertz electric field forms upon interacting with a multilayer structure.

Conclusions

Two real case applications using state of the art apparatus working in the so-called terahertz band were presented. Ultra-Wide Bandwidth FMCW radars are known for their contactless Non-Destructive-Testing capabilities for the aeronautics industry. Here, they proved to be an adequate, compact and sensitive imaging system for remote sensing of inaccessible targets through plane panel thanks to the object visualization and identification as well as alteration and defect detections. Further ongoing works are oriented toward the adequacy of this technique applied to other types of materials. The implementation of a more versatile FMCW radar system together with other optical NDT tools to reach a resourceful imaging system is also under development and is intended to be integrated with an automated positioning system. For the potential of terahertz time domain measurements on painting thickness, a combination of numerical strategies was described to extract the parameters that describe the best a measured signal arising from a multilayered sample. From simulated reflected waveforms, through Fresnel's equations, an iterative and genetic algorithm then selects the highest contributing optical path (depending on the material response and its thickness) amongst the stack for the generation a reduced transfer function. This new predictive tool is crucial to analyse any deviation for the awaited signal and to identify which layer does not follow the right thickness or dielectric properties.

Figure 1 :

 1 Figure 1: (a) Typical radar unit architecture (b), frequency sweep sensing scheme for FMCW measurement

Figure 2 :

 2 Figure 2: Optical configuration diagram for radar imaging in reflection mode with intermediate covering panel.

Figure 3 :

 3 Figure 3: (a) 100 GHz (b) and 300 GHz radar scan in reflection mode of a glass fiber composite test sample displaying inclusions, (c) picture of the corresponding sample.

Figure 4 :

 4 Figure 4: (a) 100 GHz image of a wrench and sponges with their support through a plane covering panel (b) Synview imaging system with implemented plane cover (side view and top view) metallic samples and sponges are visible on the low right side (c) zoom on a 100 GHz Image of an altered aluminum tube that displays 2 mm (c1) and 5 mm (c2) cracks on its right part (dashed frame) , through a plane covering panel.

  . Eventually, a time-of-flight analysis is conducted to evaluate the respective thicknesses.

Figure 5 :

 5 Figure 5: Scheme of the iterative algorithm. At each further dielectric interface, the k-parameter is incremented and the corresponding transfer function H is calculated.

Figure 6 :

 6 Figure 6: Iterative reconstruction (red) of the experimental waveform (black) for the aeronautic stack made of four layers, each of tens of microns.

Figure 7 :

 7 Figure 7: Profile of the integrated comparison metric between 22 ps and 37 ps as a function of the number of iterations.

Figure 8 :

 8 Figure8: Representation of a portion of the selected optical routes found to be necessary to reconstruct the experimental waveform.

Figure 9 :

 9 Figure 9: Reconstruction of the signal provided by the restricted selection of optical paths (red) in comparison to the previously reconstructed waveform by means of all propagation routes (black). Divergences appear after 32 ps, where additional subdivided pulses have been deliberately omitted.

Acknowledgements

The authors thank Nouvelle Région d'Aquitaine and BPI France for Athermo project funding. They also thank ATR company for providing real sample for the demonstration. The authors would like to thank OPTIKAN SAS for fruitful discussions.

Supplementary materials

An animation of the FMCW radar 3D scan for objects detections through the airplane cover is attached to the online version


R Fillet

V Nicolas

V Fierro

A Celzard
email: alain.celzard@univ-lorraine.fr

Modelling heat and mass transfer in solar evaporation systems

Keywords: Simulation, Biosourced material, Solar energy, Solar evaporation, Porous media, Heat and mass transfers

A two-dimensional axisymmetric model of a solar evaporator was developed, based on a rigid biosourced open-cell foam floating on the surface of water contained in a beaker, including a set of boundary conditions and water transport in the porous medium.

Experiments were carried out to validate the model for two cases, with and without a layer of insulating material around the beaker. Several concrete cases of experimental solar evaporation studies recently reported in the literature, with or without insulation, radiation intensification, etc., have then been studied in the light of our model. We have thus been able to compare the observed trends with those of the numerical model in order to explain the results reported in the literature and, in some cases, to find a geometric optimum. The effect of ambient conditions on the evaporation rate was investigated, and significant changes were observed by modifying solar input, air temperature and air humidity. We also showed that increasing the evaporation area and using an insulating layer between the water and the evaporative material increases the evaporation performance.

Introduction

In modern society, energy plays a major role, and its development as we know it today is only possible with access to a significant amount of resources. Nowadays, improving the effectiveness of systems and producing energy from more sustainable sources is more relevant than ever. Water evaporation, for example, is a subject of great interest [START_REF] Zhao | Materials for solar-powered water evaporation[END_REF] and a wide variety of materials have recently been developed to produce clean water from any source, even polluted water and seawater, without the need for an energy grid. Evaporation indeed occurs with or without solar flux at the water surface, but it can be enhanced by adding a floating porous medium.

Water evaporation systems (WESs) are composed of one or more materials, sometimes synthesised or manufactured from a natural resource, whose characteristics can be modified to improve evaporation efficiency [START_REF] Zhang | Direct solar steam generation system for clean water production[END_REF]. In general, studies aim to improve materials properties, e.g. to increase the absorption of sunlight, to improve the energy conversion of light into heat or to improve the transport of water to the surface. Thermal conductivity can also be modified to reduce transfer to bulk water and thus conserve energy at the surface where evaporation occurs. In addition, the resistance of materials to clogging has been studied in seawater evaporation processes, where salt can affect long-term performance [START_REF] Kuang | A High-Performance Self-Regenerating Solar Evaporator for Continuous Water Desalination[END_REF]. Changes in geometry have also been made to exploit heat exchange from the air to the surface instead of letting the transfer occur from the hot evaporation surface to the air [START_REF] Song | Cold Vapor Generation beyond the Input Solar Energy Limit[END_REF].

Besides the large variability of solar evaporation systems mentioned above, a common problem encountered when reading the literature is that the protocol for measuring experimental performance is not identical in all publications. For example, some studies report mass loss measurements under a transparent bell jar and others do not. The bell jar has the effect of increasing the temperature and humidity of the air around the material compared to an open air system. In addition, tests can be carried out in different locations around the world with different humidity and temperature conditions. In some cases, an insulating layer can be used to support the evaporation material, which would not necessarily be able to float on its own, and to provide thermal insulation with the bulk water. Examples of such a system have been published based on pomelo peal [START_REF] Zhang | Converting Pomelo Peel into Eco-Friendly and Low-Consumption Photothermic Biomass Sponge Towards Multifunctional Solar-to-Heat Conversion[END_REF] or plant stems [START_REF] Zhang | Vertically aligned Juncus effusus fibril composites for omnidirectional solar evaporation[END_REF]. All of these differences lead to more or less favourable conditions for the evaporation of water, regardless of the material.

The interest in proposing a numerical model that takes into account the mass and heat transfer that underlies water evaporation is the ability to accurately capture the behaviour of such WESs, and to explore the effect of experimental parameters such as those presented above. Whether it is the layout of the system, the geometry of the material or even the properties of the water within the material, it is possible to examine all of these impacts on the overall performance and behaviour in order to better understand how to design an efficient evaporation system. Furthermore, although numerical simulations have been used in the past for solar evaporation systems, they were limited to predicting certain aspects such as temperature distribution or water diffusion [START_REF] Wu | Highly efficient three-dimensional solar evaporator for high salinity desalination by localized crystallization[END_REF][START_REF] Zhao | Highly efficient solar vapour generation via hierarchically nanostructured gels[END_REF][START_REF] Guo | Biomass-Derived Hybrid Hydrogel Evaporators for Cost-Effective Solar Water Purification[END_REF]. As far as we know, a general model taking into account the different types of heat exchange between the system and the air, the transport of water in the material and the evaporation from the water surface is missing in the literature.

In this work, we developed a model based on experimental results obtained under strictly controlled conditions, using a mimosa tannin foam that we have synthesised and used to improve water evaporation. Two cases were experimented, with and without insulation around the vessel containing the water to be evaporated, to limit heat exchange at the boundaries. These two cases were simulated, as the physics at the boundaries are different, and were validated using the same model apart from the system geometry and a few parameters. Various cases from the open literature have also been considered, allowing the effects of various operating conditions and materials geometry to be explored, and to which our model has been applied. Although the results between the literature experiments and our simulations cannot be directly compared due to the different materials used, the main trends could be successfully compared. Therefore, our model should be particularly interesting in revealing the optimal geometry parameters, such as the dimensions of the insulation to be used beneath a material, or to study the impact of the insulation provided by synthetic foams used to allow the flotation of the material. Prior to our research, a guideline for the main experimental conditions, such as ambient temperature, relative humidity, geometry or position of the different elements, had been discussed [START_REF] Li | Measuring Conversion Efficiency of Solar Vapor Generation[END_REF]. Here, these conditions are investigated to describe their impact on the final evaporation rate and efficiency results, as our model takes into account the ambient air temperature, relative air humidity, air convection type and solar flux.

Materials and Methods

Evaporation material

Firstly, a bio-based foam was prepared, from which water evaporation experiments were conducted to validate our model. This foam was synthesised from tannin, a natural extract from the bark of the mimosa tree (Acacia Mearnsii), using a process described elsewhere [START_REF] Tondi | Tannin-based carbon foams[END_REF]. The material was obtained by physical foaming, which involves producing a gas in situ in a thermosetting resin during its polymerisation. Table 1 shows all the ingredients used in the process, indicating their relative amounts and their role in the formulation. Briefly, tannin is first dissolved in water under magnetic stirring and then mixed with furfuryl alcohol. Next, diethyl ether and a 37 wt. % aqueous solution of formaldehyde are added and the mixture is carefully homogenised. Finally, a 65 wt. % aqueous solution of para-toluenesulphonic acid (pTSA) is added under stirring, after which the mixture expands rapidly while polymerising at the same time. Once fully expanded and cured, the foam is dried at room temperature for one week; further details are given elsewhere [START_REF] Marie | First approach for modelling the physical foaming of tannin-based thermoset foams[END_REF]. The foam is then cut to the proper dimensions, as shown in Fig. 1 (about 58 mm in diameter and 15 mm thick).

Experimental protocol and measurement methods

The experiment consists in illuminating a beaker (60 mm diameter), filled with tap water and with the foam floating on the surface, with a solar light simulator (SOLAR LIGHT LS-1000). The experimental setup is shown in Fig. 1. A Sartorius ENTRIS weighing scale accurate to the milligram is used to measure the mass of the whole system, and it is thus possible to determine the amount of evaporated water by difference. T-type thermocouples are used to measure temperatures at the interface between the evaporation material and the water on the one hand, and within the evaporation material on the other hand. A thermal camera (OPTRIS PI) pointing at the surface of the evaporation material is used to measure its temperature locally. Air temperature and relative humidity are measured using a sensor placed in the room. To explore the effect of heat exchange from the sides, a second experiment was also conducted with insulating foam around the beaker. Therefore, case 1 discussed below will refer to the case without insulation and case 2 to the case with insulation. In order to compare it with the numerical model, the experiment was carried out under three different solar concentrations: 0.7, 1.0 and 1.3 sun (1 sun being equal to 1 kW m -2 ) for about one hour for each illumination. This allowed the evolution of heat and mass to be studied as a function of solar concentration. The system efficiency was calculated using Eq. 1:

Eq. 1 where (kg m -2 h -1 ) is the evaporation rate, (sun) is the solar concentration, and

(1000 W m -2 ) is the power density of the average solar light on earth. Typical results in the literature present the evaporation rate and efficiency for a given specific solar power.

Assuming a theoretical efficiency of 100 % and with 1 sun illumination contributing only to the evaporation of water, whose enthalpy of vaporisation is 2450 kJ kg -1 under STP conditions, is equal to 1.47 kg m -2 h -1 , which therefore represents the theoretical limit of evaporation. Sometimes in the literature, evaporation in the dark is subtracted from the evaporation rate with illumination to take into account only evaporation due to light. It is useful to note that Eq. 1 does not take into account the heat exchange with the air, but only the solar illumination input and the evaporation rate output. This is deliberate, as Eq. 1 is most often used in the literature, and this is how we can compare our results with what has already been published.

Physical model

A 2D axisymmetric model was written using COMSOL Multiphysics 5.6. Bulk water and evaporation material were modelled as rectangles (in blue and red, respectively, in Fig. 2).

This basic geometry was adapted to cases 1 or 2 of the study, i.e., without or with an insulating layer around the beaker, respectively. In this model, the upper surface is illuminated by a solar radiation, is cooled by convection with air, and evaporates water. On the sides and bottom, there are only convection and conduction, respectively.

General equations

In order to simulate the heat and mass transfers existing during the evaporation process, energy and mass conservation equations were used to obtain the temperature, T, (Eq. 2) and the volume fraction of liquid in the pores, (Eq. 3).

The conservation of energy equation describes conductive heat transfer and temperature changes over time. The energy associated with mass transfer is small and due to the very slow kinetics of this phenomenon, it has therefore been neglected in this study. This conservation of energy equation is applied to calculate the temperature in the beaker and in the insulator (corresponding to case 2).

Eq. 2

where (kg m -³) is the density, (J kg -1 K -1 ) is the heat capacity, and (W m -1 K -1 ) is the thermal conductivity of component j (water w, evaporation material m or insulating layer i).

The mass transfer is described by a phenomenon of water diffusion in the porous material and is not solved in water or in the insulated medium. A transient equation is used in order to allow the material to dry out and thus reduce surface evaporation when it is no longer saturated. The equation is obtained by keeping the water liquid inside the material (Eq. 3).

Eq. 3 where (dimensionless) is the water volume fraction in the evaporation material, (dimensionless) is its skeleton volume fraction, (kg m -³) is the density of liquid water, (kg m -³) is the skeleton density of the evaporation material, and (m² s -1 ) is the diffusion coefficient of water in the evaporation material.

Boundary conditions

In this section, the boundary conditions used are specified. It should be noted that thermal continuity conditions are used for the heat equation for all internal boundaries.

Interface between the air and the material

The heat exchange at the surface (Eq. 4) takes into account convection, incoming light radiation, as well as losses by emission and by water evaporation (Fig. 2). The latter involves the surface evaporation, (kg m -2 h -1 ), multiplied by the energy consumed to raise the temperature of water and the enthalpy of vaporisation of water, (J kg -1 ).

Eq. 4

with (W m -1 K -1 ) the thermal conductivity of the evaporation material, (W m -2 K -1 ) the convection coefficient at the surface, (J K -1 kg -1 ) the heat capacity of water, (K) the initial temperature of water, (1000 W m -2 ) the standard solar surface power, (sun) the solar concentration (used to achieve different values of solar power), (dimensionless) the emissivity of the material and (5.67 × 10 -8 W m -2 K -4 ) the Stefan-Boltzmann constant, and the unit vector normal to the boundary.

To calculate the convective heat transfer coefficient, the following equation was used:

Eq. 5

where L c is the characteristic length of the surface (which corresponds to the diameter for the top surface and the height for the side), and Nu is the Nusselt number (plane surface assumption), which can be calculated as follows:

Eq. 6

In Eq. 6, Re is the Reynolds number and Pr the Prandtl number, such that: Eq. 7

in which is the thermal diffusivity and the dynamic viscosity of the fluid, here the air. The Reynolds number is calculated using the following equation:

Eq. 8

with u the velocity of air (m.s -1 ) and is kinematic viscosity (m 2 .s -1 ).

Since is fitted for the case without any forced convection, the Nusselt number for the case of natural convection, , is determined using Eq. 5. Then, a Nusselt number , combination of both forced and natural convection, can be calculated [START_REF] Cengel | Heat Transfer a practical approach[END_REF]:

Eq. 9

At the upper surface, the water loss in the diffusivity equation is equal to the evaporation rate:

Eq. 10

Finally, since water evaporation comes from the difference between the mass concentration of the surface and the air, then is calculated as follows [START_REF] Nicolas | Modeling bread baking with focus on overall deformation and local porosity evolution[END_REF]:

Eq. 11

where (18 × 10 -3 kg mol -1 ) is the molar mass of water, (8.314 J mol -1 K -1 ) is the universal gas constant, (Pa) and (Pa) are the water vapour pressures at the local temperature (K) and at the air temperature (K), respectively, is the water activity, (%) is the relative humidity and (m s -1 ) is the convective mass exchange coefficient, calculated using the Colburn analogy [START_REF] Cengel | Heat Transfer a practical approach[END_REF]:

Eq. 12 and are calculated from the Clapeyron formula as follows:

Eq. 13

where is the water vapour pressure at the temperature of the environment, is the standard pressure of 101 325 Pa and is the boiling point of water (373.15 K).

Interface between the air and the side of the system

Heat exchange between the air and the wall of the beaker involves convection and radiation losses:

Eq. 14 where (W m -1 K -1 ) has the same meaning as in Eq. 1, (W m -2 K -1 ) is the convection coefficient at the side, and (dimensionless) is the emissivity of the side.

For forced convection at the side of the system, the following Nusselt expression is used assuming a vertical cylinder:

Eq. 15

Interface between the weighing scale and the bottom of the beaker

At the interface between the weighing scale and the bottom of the beaker, only a convective heat flow was assumed, formalised by: Eq. 16

where (W m -2 K -1 ) is the convection coefficient at the bottom.

Interface between the water and the evaporation material

At the interface between the water and the evaporation material, as the latter is considered saturated, the liquid fraction is fixed at a value equal to the porosity, , and thermal continuity is assumed, i.e.:

Eq. 17

Simulation parameters

Table 2 lists the input parameters used in the simulation, and their symbols, definitions, units and values considered are also given.

,

Heat transfer coefficients at the surface, the side and the bottom W m -2 K -1 5, 5, 1

Emissivities of the material and of the side of the beaker -0.8, 0.8

Water activity

Fig. 3 shows the water activity, a w , as a function of water uptake (wt.%). a w was obtained from the experimentally obtained moisture sorption isotherm, taken from a previous study [START_REF] Delgado-Sánchez | Hydrophobisation of tannin-based foams by covalent grafting of silanes[END_REF], and corresponds to the partial vapour pressure of water divided by the standard vapour pressure of water. From the knowledge of the moisture content inside the material, it was possible to obtain the partial vapour pressure of water with the water activity, and thus to calculate the evaporation rate using Eq. 11. Since the evaporation rate is related to the amount of water at the surface, the WES performance is limited if there is not enough water at the surface. At 31% water uptake, the water activity is 1, and remains constant for higher water uptakes.

Fig. 3.

Water activity in the tannin foam as a function of water uptake (after data from [START_REF] Delgado-Sánchez | Hydrophobisation of tannin-based foams by covalent grafting of silanes[END_REF]).

Thermal conductivity

Since the porous evaporation material is partially filled with water, its average thermal conductivity (W m -1 K -1 ) must be calculated taking into account the thermal conductivities of the material skeleton, water and air:

Eq. 18

where is the volume fraction of solid in the porous evaporation material, is the corresponding volume fraction of gas such that , and , and (W m -1 K -1 ) are the thermal conductivities of water, material skeleton and air.

Heat capacity

The mass heat capacity of the material, , varies much with the presence of water but the effect of air on can be neglected. Therefore, the volumetric heat capacity was calculated as follows:

Eq. 19

where (kg m -3 ) and (J kg -1 K -1 ) are the density and the mass heat capacity of the simulated material filled with water, respectively, and (J kg -1 K -1 ) is the mass heat capacity of the dry, porous, evaporation material.

Vaporisation enthalpy

The enthalpy of vaporisation of water, (kJ kg -1 ), changes slightly with the temperature according to [START_REF] Henderson-Sellers | A new formula for latent heat of vaporization of water as a function of temperature[END_REF]:

Eq. 20 16

Numerical model

Calculated variables T and were obtained using the commercial software COMSOL Multiphysics 5.6. First, a preliminary parametric study on the triangular mesh was carried out in the first study case (i.e., without insulation) for three meshes comprising 2569 elements (mesh #1), 2873 elements (mesh #2), and 5553 elements (mesh #3). Since the impact of mesh size on the results were extremely limited, mesh #2 with 2873 elements was used. Fig. 4a shows a mesh consisting of 2873 elements for the first study case and 5431 for the second (Fig. 4b), both refined at the interfaces between water and evaporation material, as well as the interface between air and the evaporation material, with 100 elements. The rest of the mesh is created automatically by the software, depending on the physics involved. Using a desktop computer with a quad-core Xeon E5 coupled with 64 GB of RAM, the model for a system without insulation is calculated in 1 minute and 5 seconds, while the model for a system with insulation around the beaker is calculated in around 1 minute and 40 seconds. The simulation required 90 calculation steps, with a step time self-determined by the software. with insulation (case 2).

Results and discussion

In order to validate the numerical model, the results of the calculations were compared with the results of the experiments obtained with the floating tannin-based foam, for the two cases, i.e., with and without insulation foam. Then, different cases from the open literature were examined with the proposed model, and the model was modified if necessary, taking into account the different parameters considered in these former studies. This includes ambient air temperature, relative humidity, solar power density, convection coefficient, and system geometry. Overall, the objective is not to reproduce the same results from the literature but to compare trends.

Validation of the model

Evaporation experiments were conducted without and with insulation (i.e., corresponding to case 1 and 2, respectively) around the beaker containing the floating tannin foam, under different solar inputs. To validate the simulation, the external conditions were imported as input ( and ) and the duration was set to that of the experiment.

Non-insulated configuration

Fig. 5a shows the vertical section of the configuration without insulation, where 1, 2 and 3 are the temperature measurement points, and Fig. 5b and Fig. 5c show the experimental and predicted results. The experimental evaporation rates are equal to 0.51, 0.70 and 1.01 kg m 2 h -1 with efficiencies of 49, 47 and 52 % for 0.7, 1 and 1.3 sun, respectively. In comparison, the calculated evaporation rates are equal to 0.48, 0.74 and 1.03 kg m 2 h -1 with efficiencies of 46, 50 and 53 % for 0.7, 1 and 1.3 sun, respectively. The difference between the simulated mass of water and the experimental measurements is no higher than 5% (see Fig. 5c). The three different solar inputs are clearly visible on the temperature curves, showing three successive jumps. The agreement between measured and calculated temperatures is better in the evaporation material (point 2) and at the interface between the latter and the water (point 3)

than at the surface (point 1) (see Fig. 5b). This discrepancy at point 1 is attributed to the setting of the thermal camera (position, emissivity, measurement area) and the estimation of fitted parameters in the model such as the convection coefficient. The average differences between the simulated and experimental temperatures in the first case study are indeed, from top to bottom, 1.01, 0.69 and 0.44 °C, which shows that the model reproduces the experiment fairly closely. sun, respectively. In comparison, the calculated evaporation rates are equal to 0.54, 0.88 and 1.25 kg m 2 h -1 with efficiencies of 52, 59 and 64 % for 0.7, 1 and 1.3 sun, respectively. The average difference between experimental and simulated masses of water is 8%, so it is higher than in the non-insulated configuration. The differences for the three temperatures were 2.71, 1.31 and 0.64 °C for points 1, 2 and 3, respectively (Fig. 5d). As in the previous case (Fig. 5b), and for the same reasons, the difference between experimental and simulated surface temperature is higher (point 1 in Fig. 5d) than for the other two temperatures, in the bulk (point 2) and at the interface (point 3), which were recorded with thermocouples. However, the evolution of the water mass over time is fairly well predicted by the model for the three illuminations tested, as well as the transient thermal evolution that occurs when the input sunlight changes.

Insulated configuration

The numerical model has therefore been validated and we can conclude that it can be used to simulate water evaporation in the following. This validation was done on the basis of the results obtained with tannin foam, but this material did not show particularly interesting performances compared to other materials reported in the literature, and for which the evaporation rate were close to 1.3 kg m -2 h -1 [START_REF] Fillet | A review of natural materials for solar evaporation[END_REF]. We now extend the application of this model to these more efficient materials in terms of solar evaporation.

Study of different cases using the model

Among the large number of evaporation results for different materials and systems in the literature, some phenomena and techniques are discussed repeatedly. The aim of this section is to test the behaviour of our model when applied to the cases seen in the literature. Those cases can be a good source of inspiration for the design of an efficient evaporation system.

For all studies considered, unless otherwise specified, the temperature of the air was set at 20

°C with a relative humidity (H R ) of 50 %. The initial temperature of all components was also set at 20 °C.

Impact of ambient conditions

The previous section 4.1.2 has shown the significant impact of the beaker insulation on the system performance, which indicates a dependence of the performance on the ambient conditions. In order to better understand their impact, several parametric studies on the exchanges between the air and the material (i.e., air temperature, relative humidity, convection and solar input) were carried out using the numerical model. The evaporation performance was calculated from the final results, which correspond to a steady state.

Air temperature

The first parameter that has an impact on the evaporation rate is the ambient air temperature, which can be easily measured and controlled. In the simulation, we set the relative humidity at 50 % with natural convection, and a parametric study was conducted on the air temperature.

Fig. 6a shows the simulated evaporation rates for different values of air temperature, and the corresponding evaporation efficiency. With increasing temperature, the evaporation rate increases significantly, leading to changes of about 0.15 kg m -2 h -1 and 10 % efficiency for a temperature change of 10 °C, which may be a reasonable change over a whole day.

Evaporation in the dark is also likely to change, as it is also shown in the same figure for three different temperatures (taken from [START_REF] Song | Cold Vapor Generation beyond the Input Solar Energy Limit[END_REF]). °C (from [START_REF] Song | Cold Vapor Generation beyond the Input Solar Energy Limit[END_REF]).

Relative humidity

The relative humidity (H R ) is considered in Eq. 11 and it also has an impact on the evaporation rate at 20 °C, which was studied by varying the H R from 10 to 90 %, with a 10 % step. Fig. 6b shows the evaporation rate and efficiency as a function of H R . Lower H R increases the evaporation rate and efficiency, for instance the evaporation rate increases from 0.60 kg m -2 h -1 , with an efficiency of 40 % at 90% H R , to 0.72 kg m -2 h -1 , with an efficiency of 49% at 10 % H R . Although H R is a factor that has a significant impact on evaporation rate and efficiency, it is easier to control and can simply be recorded as an ambient parameter. The effect of H R at 23.6 °C on the evaporation rate in the dark is also plotted for three H R values for reference (from [START_REF] Song | Cold Vapor Generation beyond the Input Solar Energy Limit[END_REF]).

Heat exchange with air by convection

Insulating the system increases evaporation performance by increasing the temperature and thus the evaporation rate [START_REF] Wang | Reversing heat conduction loss: Extracting energy from bulk water to enhance solar steam generation[END_REF], as shown when validating our model in section 4.1.2. Although convection increases the heat exchange of the system with the air [START_REF] Ali | Challenges and opportunities in functional carbon nanotubes for membrane-based water treatment and desalination[END_REF], which means that there is more heat loss (Eq. 4, Eq. 14 and Eq. 16), convection also increases the evaporation rate (Eq. 11 and Eq. 12) with or without solar illumination. Here, three types of forced convection, with and without insulation, are tested by adjusting the convection coefficients to three air flow rates. These three flow rates are 0.15, 0.30 and 0.6 m s -1 which corresponds to a heat exchange coefficient at the evaporation surface h c of 6.8, 9 and 12.4 W.m -2 .K -1 . Fig. 7 shows the surface evaporation rates in the dark and under 1 sun illumination with their corresponding efficiencies as a function of the air flow rate, with and without insulation around the beaker. Forced convection improves water evaporation performance in the dark and in the light, which means that although heat loss increases with convection, the evaporation rate benefits more. Insulation around the beaker also improves evaporation in darkness and light. Despite slightly reducing the slope of the trend line (compare Fig. 7b and Fig. 7a), beaker insulation gives the same results by increasing convection. As shown in Fig. 1c and Fig. 7d, the insulation reduces the amount of heat exchange between the material and the bulk water as well as between the material and its side, which is divided by 3 in the insulated case. The global heat exchange between the air and the evaporation surface is lower in the insulated case than in the non-insulated case. Considering that the solar input is the same, this is consistent with the evaporation results. The more water evaporates, the more heat is transferred to the air, and the less heat reaches the material surface. Overall, insulation can only reduce the heat loss from the sides of the beaker, , which does not have as much impact on evaporation as convection at the surface of the material. Nevertheless, convection may improve evaporation in a large open area when a much larger evaporation surface is used, so it is an interesting factor to consider when evaluating evaporation performance.

Solar power and efficiency

One of the most mentioned phenomena in the literature is the increase in evaporation efficiency with increasing solar energy [START_REF] Sheng | Bamboo decorated with plasmonic nanoparticles for efficient solar steam generation[END_REF][START_REF] Chen | Highly Flexible and Efficient Solar Steam Generation Device[END_REF]. Many papers show results with illuminations above 1 sun as solar energy input [START_REF] Jiang | Bilayered Biofoam for Highly Efficient Solar Steam Generation[END_REF][START_REF] Zhu | Plasmonic Wood for High-Efficiency Solar Steam Generation[END_REF][START_REF] Zhu | Tree-Inspired Design for High-Efficiency Water Extraction[END_REF][START_REF] Li | Broadbandabsorbing WO3-x nanorod-decorated wood evaporator for highly efficient solar-driven interfacial steam generation[END_REF][START_REF] Ghafurian | Wood surface treatment techniques for enhanced solar steam generation[END_REF][START_REF] Yu | Enhancing efficiency of carbonized wood based solar steam generator for wastewater treatment by optimizing the thickness[END_REF][START_REF] Lewis | Powering the planet: Chemical challenges in solar energy utilization[END_REF]. Fig. 8 shows evaporation rate and efficiency results from the literature as a function of solar power concentration. The more power there is on the WES, the higher is the efficiency.

Fig. 8 Evaporation rates (E.R.) and evaporation efficiencies for different solar concentrations, as reported in the literature: (a) and (c) for a wood-based material [START_REF] Sheng | Bamboo decorated with plasmonic nanoparticles for efficient solar steam generation[END_REF] (reproduced from Ref.

[20] Copyright (2017), with permission from Wiley); and (b) and (d) for an improved bamboo-based material [START_REF] Chen | Highly Flexible and Efficient Solar Steam Generation Device[END_REF] (reproduced from Ref. [START_REF] Chen | Highly Flexible and Efficient Solar Steam Generation Device[END_REF] Copyright (2020), with permission from Elsevier).

Based on these articles, Fig. 9a shows the simulated evaporation rate and efficiency for our system. The efficiency increases from 49 to 67 % when the illumination increases from 1 to 7

sun and so the evaporation rate increases from 0.7 to 6.9 kg m -2 h -1 . It is important to note that, according to Eq. 1, if the efficiency had remained constant, the evaporation rate at 7 suns would have been 4.9 kg m -2 h -1 (only 7 times the evaporation at 1 sun). With a 15 % increase in efficiency from 1 to 7 sun for our system, Fig. 8c showing a 12 % increase and Fig. 8d a 14 % increase for both systems from the literature, the increase in efficiency with the solar power concentration is consistent. Natural convection alone cannot be the only cause of the increase in efficiency and, although the enthalpy of vaporisation also changes with the temperature according to Eq. 20, it produces only a slight increase in the evaporation rate. The main cause of the increase in efficiency is the surface temperature, which increases steadily with increasing illumination.

As is a function of temperature (Eq. 13), the latter plays a major role in the evaporation rate (Eq. 11). Fig. 9b shows that the evaporation rate increases exponentially with surface temperature while it increased linearly with illumination (Fig. 9a).

Increase of the exchange area with air

Several studies have shown that increasing the evaporation area effectively increases the performance of the system [START_REF] Song | Cold Vapor Generation beyond the Input Solar Energy Limit[END_REF][START_REF] Wu | Highly efficient three-dimensional solar evaporator for high salinity desalination by localized crystallization[END_REF][START_REF] Wu | Solar-driven evaporators for water treatment: challenges and opportunities[END_REF]. One of its consequences is a decrease in surface temperature, due to the fact that sunlight is then scattered over a larger area and the exchange with air is higher, thus reducing the local heating and thus the temperature. This technique has already been explored in the literature [START_REF] Song | Cold Vapor Generation beyond the Input Solar Energy Limit[END_REF], as shown in Fig. 10, with a carbon-coated paper placed on a polystyrene foam. By increasing the receiving area (Fig. 10d and Fig. 10e) the evaporation rate also increases (Fig. 10c). Copyright (2018).

Based on this work, we have modified the geometry of our model to represent the increase in surface area in the simulation. A conical geometry, not identical but close to the geometry illustrated in Fig. 10a, was chosen for its simplicity to be captured by the axisymmetric model. With such a geometry, we can easily conduct a parametric study on the height of the cone. For that purpose, a variable N has been created, which is written as follows:

Eq. 21 [START_REF] Wu | Solar-driven evaporators for water treatment: challenges and opportunities[END_REF] where h (m) and r (m) are the height and the radius of the base of the cone, respectively. This allows the cone height to be easily fixed in the simulation and to produce a parametric study with integer values. The solar power, initially 1 kW m -2 at the base of the cone, is decreased when the tip of the cone is stretched as this increases the area. A coefficient (dimensionless) is defined to account for the corresponding scattering of the solar power. It represents the increase in surface area when changing from a flat surface to a conical surface and is calculated as follows:

Eq. 22 where (m 2 ) is the area of the base of the cone and (m 2 ) is the lateral area of the cone. Therefore, the power concentration at the surface of the cone reads:

Eq. 23 where (kW m -2 ) is the surface solar power arriving on the material, and (kW m -2 ) is the surface solar power arriving from the sun. A with the decrease in power concentration, the evaporation rate increased, as the evaporation rate is still determined by the flat surface that receives the sunlight ( ). Fig. 11 shows the conical geometry used in the simulation with the parameters introduced previously, for three values of N. With the increase of the parameter N, k decreases, which reduces the solar power density. Fig. 12a shows the simulation results of evaporation rates and efficiencies by varying N from N = 0 (flat surface) to N = 8, applying the model presented here. The temperatures are also shown for each simulated cone. It can be observed that the evaporation rate and the efficiency increase linearly with N while the surface temperature decreases. This observation is contrary to the previous study, which showed that the evaporation flux increases with the surface temperature. It can also be observed that the evaporation rate decreases with increasing the area of the cone (Fig. 12b). The calculation of the performance of the system, concerning evaporation, is not related to the exchange area but to the area at the base of the structure. By calculating the flux related to the base, we observe an increase in the evaporated flux as a function of N and therefore of the exchange area (Fig. 12b). It is also observed that the ratio 1/k between the exchange area and the area at the base increases linearly with N (Fig. 12c), which shows the direct relationship between the evaporation area and the evaporated flux related to the base of the cone. Overall, the increase in evaporation area is beneficial, since it doubles the efficiency from N = 0 to 8, passing from 44 to 88%. It should be noted that these efficiencies were calculated using Eq. 1, which does not take into account the heat exchange with the air, but only the solar illumination input and the evaporation rate output.

However, such an improvement could have limitations, for example, a slender geometry with N = 8 could be difficult to produce. Another possible geometry could be a recessed spherical area, which might be easier to manufacture [START_REF] Li | Simple Hierarchical Interface Design Strategy for Accelerating Solar Evaporation[END_REF]. With such a geometry, the surface allows the same kind of reduction of heat loss with air, but also a reduction of radiation loss, since the light beam is partially reflected inside the sphere. Such a system could achieve 1.65 kg m -2 h -1 for a solar illumination of 1 sun [START_REF] Li | Simple Hierarchical Interface Design Strategy for Accelerating Solar Evaporation[END_REF].

Insulation under the material

Due to the natural shape of some materials, it is not possible to use them alone for evaporation and therefore an additional material must be used to maintain the evaporative surface above the water [START_REF] Zhang | Converting Pomelo Peel into Eco-Friendly and Low-Consumption Photothermic Biomass Sponge Towards Multifunctional Solar-to-Heat Conversion[END_REF][START_REF] Xu | Mushrooms as Efficient Solar Steam-Generation Devices[END_REF][START_REF] Lin | Poly(Nphenylglycine)-Based Bioinspired System for Stably and Efficiently Enhancing Solar Evaporation[END_REF][START_REF] Liao | Lotus leaf as solar water evaporation devices[END_REF][START_REF] Zhu | Biomass Carbon Materials for Efficient Solar Steam Generation Prepared from Carbonized Enteromorpha Prolifera[END_REF][START_REF] Wang | Functionalized biomass-derived composites for solar vapor generation[END_REF]. This can be a standard foam for thermal insulation, which has the ability to float on water but also allows less heat exchange with water due to its very low thermal conductivity. This was done by Xu et al. [START_REF] Xu | Mushrooms as Efficient Solar Steam-Generation Devices[END_REF], whose system is shown in Fig. 13a, where a mushroom was placed on a polystyrene foam support. The path taken by the water to evaporate is in the centre of the insulation foam, where a hole has been made to allow the mushroom to touch the water. Similarly, Fig. 13b shows another type of WES that has paper wrapped around an insulating foam [START_REF] Li | Graphene oxide-based efficient and scalable solar desalination under one sun with a confined 2D water path[END_REF] to allow water to be transported through the sides of the system to an evaporation surface at the top. Inspired by these two uses of insulating foam, we studied two configurations, one where the insulator is located in the core of the evaporation material, with water transport taking place at the sides, and the other where the insulator is located around the material, with water transport taking place through a hole in the foam. These configurations are shown in Fig. 14a and Fig. 14b, respectively. Thus, the evaporation material is positioned on the insulating/floating foam, and water transport is allowed either through a central hole or through the periphery. To study the impact of the foam, parametric studies were conducted on the effects of varying the radius (R i ) and the height (h i ) of the insulating foam. and represent the ratios of the radius and heights of the foam to those of the evaporation material (see Eq. 24 and Eq. 25 and Fig. 15).

Eq. 24

Eq. 25

where R t and h t are the radius and the height of the evaporation material, respectively. The higher and , the larger the foam and the more space it takes under the material.

with thus correspond to the case without insulation layer. The variations of N r for the two configurations are represented in Fig. 15.

Adding an insulation layer under the evaporation material should theoretically improve performance as more heat is retained on the surface. However, the addition of insulation limits the transport of water to the surface, so a balance must be struck between heat and water management. The simulation shows that increasing the size of the insulation, both in radius and in height, increases evaporation to a certain point, where it reaches its maximum, and then water transport limits the performance. The diffusion coefficient of water in the evaporation material plays an important role in the dimensions of the insulation material that can be used. show the case of insulation in the centre. In this case, as shown in Fig. 15h, increasing the size of the insulating foam in both directions increases the evaporation, and as in the first case, up to a certain maximum after which it no longer increases. However, in this case, the evaporation rate starts to decrease at higher values of and . This is due to the peripheral transport of water, which spreads over a larger area compared to the hole in the centre.

Therefore, to obtain the same kind of proportion between the volume of the evaporation material and the volume of the insulator, the radius of the insulating foam must be higher.

Nevertheless, this geometry, with a thin layer around the insulator, is limited by the possibilities of preparation. A sheet of paper, like the one presented above, is feasible but a foam-like product like our tannin foam would be impossible to cut so thin. Overall, the peripheral insulation with N r = 0.9 and N h = 0.4 achieves an evaporation rate of 0.87 kg m -2 h - 1 , while the central insulation with N r = 0.95 and N h = 0.9 gives 0.85 kg m -2 h -1 .

These observations may be clues to making evaporation systems more efficient, since the thermal conductivity of evaporation materials is strongly influenced by the amount of water inside their porous structure. Having an additional insulating layer to prevent heat exchange with water results in a higher evaporation rate, when water supply to the surface is not a problem. However, this can also be an obstacle when comparing materials. Again, it is essential to consider the whole system, with all the parameters and appropriate geometry, to study its evaporation efficiency.

Conclusion

A two-dimensional axisymmetric model, validated for temperature and water mass evolutions by experiments on a solar water evaporation system based on tannin foam, has made it possible to study the impact of different operating parameters on the performance of this type of system. As present, the tannin foam is of no interest in this type of system. This is due to the fact that the open porosity of the foam is saturated with water, which considerably limits its insulating properties.

We have shown that air temperature, relative air humidity and the type of air convection around the system all have an impact on the evaporation rate. When the air temperature increases from 20 to 50 ° C, the evaporation rate is doubled, while a large change in relative humidity has very little impact on performance. In addition, a system where the beaker is insulated increases the performance of our material by 15 % on average, and increasing the convection at the evaporation surface significantly increases the evaporated flux. This study on the experimental conditions shows the need to approach the experimental conditions recommended by Li et al. [START_REF] Li | Measuring Conversion Efficiency of Solar Vapor Generation[END_REF] with the air temperature of 25 °C and the relative air humidity of less than 60 %, all in a breeze-free environment. Additionally, the cone study shows that increasing the exchange area allows the evaporation rate of the system to be increased, which is a considerable potential improvement for future development.

Finally, the addition of a thermal insulation layer to the evaporation material, with two possible positions for the water transport channel, either at the centre or at the periphery of the axisymmetric system, proved effective, providing good thermal insulation to prevent energy from going to the water while not impeding water transport. With thermal insulation at the periphery, N r of 0.9 and N h of 0.4 are the best parameters, while N r of 0.95 and N h of 0.9 is the optimal combination when the insulation is in the centre of the system.

Fig. 1 .

 1 Fig. 1. Experimental setup used to measure the evaporation rate of water with a dark porous

Fig. 2 .

 2 Fig. 2. Types of heat exchange at the boundaries of the system: (a) case 1: without insulation;

Fig. 4 .

 4 Fig. 4. Mesh used for the water evaporation system: (a) without insulation (case 1); and (b)

  The sunlight power was set to three different values during the same test: 0.7 sun, then 1 sun and finally 1.3 sun, each for a duration of 3300 seconds for the non-insulated configuration and 3500 seconds for the insulated configuration. Changing the illumination in the middle of the test allows the evolution over time of parameters such as temperature and mass loss to be validated when the power input changes. The results of the simulation were compared to the average of three experiments and are shown in Fig. 5.

Fig. 5 .

 5 Fig. 5. Revolution cross-section of the model water evaporation system in the case where this

Fig. 5d ,

 5d Fig. 5d, e and f are equivalent to Fig. 5a, b and c, respectively, the only difference being

Fig. 6 .

 6 Fig. 6. Simulated evaporation rate as a function of: (a) air temperature at 50% H R , with the

Fig. 7 .

 7 Fig. 7. Simulated evaporation rate: (a) without, and (b) with insulation for several air

Fig. 9

 9 Fig.9(a) Simulated evaporation rate and simulated efficiency for different solar

Fig. 10 (

 10 Fig. 10 (a) Geometries used in the literature for increasing the exchange area; (b) surface

Fig. 11

 11 Fig. 11 Conical geometry used in the simulation, here shown for three parameters N: 0 (flat),

Fig. 12 .

 12 Fig. 12. (a) Simulated evaporation rate (left axis) and efficiency (right axis) calculated using

Fig. 13 (

 13 Fig. 13 (a) Mushroom-based water evaporation system on a floating polyurethane foam

Fig. 14 .

 14 Fig. 14. Simulated temperature of the system with the insulation layer located: (a) in the core,

Fig. 15 .

 15 Fig. 15. Top: Cross-section of revolution of the evaporation system shown in Fig. 14b with

Fig. 15a ,

 15a Fig. 15a, b and e show the simulation performed when insulation at the periphery was

Table 1 .

 1 Standard formulation of tannin-based foam.

	Ingredients	Role in the formulation	Mass (g)
	Mimosa tannin extract	Base of the resin	60
	pTSA (65 wt. % in water)	Acidic catalyst	22
	Furfuryl alcohol	Co-monomer	21
	Formaldehyde (37 wt. % in water)	Crosslinker	14.8
	Water	Solvent	12
	Diethyl ether	Blowing agent	6

Table 2 .

 2 Input parameters used in the simulation.

	Symbols	Definitions	Units	Values
	, , ,	Thermal conductivities of water, skeleton of the	W m -1 K -1	0.6, 0.32,
			evaporation material, air, insulating material		0.0262, 0.043
			Porosity of the evaporation material	-	0.955
	,	,	Heat capacities of water, evaporation material	J K -1 kg -1	4180, 1800, 871
			and insulating material
	D		Coefficient of diffusion of water inside the	m 2 s -1	10 -8
			material
	, ,		Density of water, skeleton density of the	kg m -3	1000, 1550, 30
			evaporation material, and bulk density of the
			insulating material

Acknowledgements

This study was supported by the TALiSMAN project, funded by ERDF (2019-000214).