Fakhreddine Ababsa
email: fakhreddine.ababsa@ensam.eujunhui.he@ensea.fr

Junhui He

Jean-Remy Chardonnet
email: remy.chardonnet@en-sam.eul

Lispen

Free Hand-Based 3D Interaction in Optical See-Through Augmented Reality Using Leap Motion

Keywords: augmented reality, natural interaction, calibration, HoloLens, leap motion

Introduction

In augmented reality environments, the natural hand interaction between a virtual object and the user is a major issue to manipulate a rendered object in a convenient way. Microsoft's HoloLens (Microsoft 2018) is an innovative augmented reality (AR) device that has provided an impressive experience for the user. However, the gesture interactions offered to the user are very limited. HoloLens currently recognizes two core component gestures: Air tap and Bloom. To solve this issue, we propose to integrate a Leap Motion Controller (LMC) within the HoloLens device (Figure 1). We thus used 3D hand and finger tracking provided by the LMC [START_REF] Lu | Dynamic hand gesture recognition with leap motion controller[END_REF] to propose new free hand-based interaction more natural and intuitive. We implemented three fully 3D techniques for selection, translation and rotation manipulation. In this work, we first investigated how to combine the two devices to get them working together in real time, and then we evaluated the proposed 3D hand interactions.

Communication between LMC and HoloLens

LMC is attached to a separate computer to communicate with HoloLens and User Datagram Protocol (UDP) provides a real-time data transmission between the two devices. The HoloLens runs as a server, and the computer runs as a client. We use Leap Motion SDK to get hand positions at each frame. As long as a new frame received, data will be sent to the server (HoloLens). On the client side, the message reception/sending processes are running in different threads. While sending data, the client can listen to the channel at the same time. For the server, we invoke a unityevent each time to react the sending from the client. This solution makes it possible to synchronize the two devices and thus to limit the latency.

Calibration Procedure

The main idea of the calibration procedure is to get the coordinates of the same point in the two devices (LMC and HoloLens) coordinate systems. After careful consideration, we choose fingertips as the detecting points because the fingertips is easy to access. Moreover, fingers can move flexibly and can be detected easily by LMC. The entire calibration process is as follows:

1. Detect the position of fingertip and get the coordinate from LMC. 4. HoloLens records the coordinates of blue ball and LMC records the fingertip simultaneously. Sending the two coordinate to the computer. 5. After collecting a set of points, analyse the data and getting R and t using ICP algorithm [START_REF] Zhang | Iterative point matching for registration of free-form curves and surfaces[END_REF] 6.

Use the computed transformation to change the coordinates of LMC points in HoloLens coordinate frame and display by HoloLens. Compare the error between the virtual points and real ones. The calibration results are very good as shown in figure (2.e). We obtain an average error equal to 12 mm.

Free-Hand Interaction

We used the 3D coordinates of the fingertips and the centre point of the left hand in order to define three intuitive hand-based interaction techniques in HoloLens environment: a) Select by touch: when the index finger collides with the virtual cube in the HoloLens coordinates frame, the cube is selected, and its colour changes to green. b) Translate: To move the selected object the thumb and index finger of the hand must touch each other. To drop the object it is just necessary to separate the two fingers. c) Rotate: The selected virtual object is rotated according to the 3D rotation of the left hand. For that, we use the coordinates of the five detected fingertips and the centre of the hand to estimate the plane corresponding to the hand. Finally, the normal to this plane is computed, it defines the 3D rotation applied to the virtual cube when the hand turns.

Results

In order to evaluate the performance of our free-hand interaction approach, we proposed several scenarios: select different cubes with different sizes at different depths (figure 3

Conclusion

Our results indicate that combining LMC with Hololens produce viable techniques for intuitive 3D interaction, which tack into account the 3D localisation of the virtual objects, improving the sensation of presence in the augmented reality environment. Future works will continue this study and develop complex interactions that use both hands.

Figure 1 .

 1 Figure 1. LMC attached on the top of HoloLens

  2. Attach virtual plane with the real plane. (Desktop or wall) (Figure 2.c) 3. Coincide the fingertip with the marked point (Blue ball) on the plane (Figure 2.d).

  -a), Translate different cubes with different sizes (figure 3-b), rotate a 3D plane (figure 3-c). The obtained results are interesting in terms of affordance and simplicity of the interactions.


Laser-processed three dimensional graphitic electrodes for diamond radiation detectors

Benoît Caylar, Michal Pomorski, Philippe Bergonzo CEA-LIST, Diamond Sensors Laboratory, Gif-sur-Yvette, F-91191, France

We have used an original approach for diamond detectors where three dimensional buried graphitic electrodes are processed in the bulk of a diamond substrate via laserinduced graphitization. Prototype made of polycrystalline chemical vapor deposition diamond has been fabricated using a nanosecond UV laser. Its charge collection efficiency was evaluated using α-particles emitted by a 241-Americium source. An improved charge collection efficiency was measured proving that laser micromachining of diamond is a valid option for the future fabrication of three dimensional diamond detectors.

Due to its remarkable physical properties: wide band gap, high charge carrier mobility, and high radiation hardness, diamond exhibits superior interests for the fabrication of radiation detection devices 1,2,3 . Conventional diamond detectors such as beam position monitors, beam loss monitors or trackers rely on planar electrodes 4,[START_REF] Ilgner | IEEE Nuclear Science Symposium Conference Record[END_REF][START_REF] Wallny | [END_REF] , a configuration where the distance between collecting electrodes, i.e. the charge carrier drift distance, is defined by the sample thickness. But although diamond devices show extremely high radiation hardness, a degradation of the charge collection efficiency (CCE) could also be observed under prolonged irradiations of the detector, due to radiation damage and related phenomena such as polarization 7 . An alternative approach has been proposed for silicon by Parker et al 8 , where the collecting electrodes are buried within the bulk of the detector and separated by a few dozens of microns. In such a configuration the distance between electrodes, i.e. the charge carrier drift path, is not defined by the sample thickness and can be further reduced. Thus the trapping probability is lowered and higher CCE can be achieved. This geometry is called 3D geometry 8 , and a radiation hardness improvement has already been observed in Si from this approach 9,10 .

Since diamond radiation hardness is far superior to silicon for high energy particles 7,11 , one can expect for 3D structured diamond detectors to withstand higher fluences than their silicon analogues and be able to operate after unprecedented radiation levels. Thus future 3D diamond detectors might be interesting candidates for tracking applications in future colliders such as Super Large Hadrons Collider.

Here we demonstrate the possibility to use laser graphitization of diamond to produce 3D electrodes.

A nanosecond Nitrogen UV laser (MNL 100, Lasertechnik Berlin) operating at λ = 337 nm, was used to induce diamond graphitization 12,13 . The process consists in focusing a laser beam beyond the sample thickness and moving it so that the focal point moves backwards through the diamond thickness. By adjusting the stage movement, columns of precise length can be drawn through the sample (Fig. 1(b).). Laser graphitization of diamond with nanosecond pulses (τ=3ns) is a complex phenomenon: during laser irradiation, two processes are competing: conversion of diamond into graphite and graphite sublimation. The graphite generated during the beginning of a laser pulse has a higher absorption than the diamond and a lower thermal diffusion. Thus, the energy contained in the end of the pulse is mainly deposited in the newly formed thin graphite layer, leading to vaporization of this layer at high fluences. In order, to avoid sublimation of graphite a fluence of F=4,72 J.cm -2 was used.

According to the litterature 14 , this fluence level avoid graphite sublimation. Then the stage was moved with a 20µm/s velocity, which was high enough to avoid multiple shots on the same graphite part, but also low enough to grow continuous graphitic wires.

On Fig. 1(b) one can notice large cracks surrounding the electrodes. Indeed when part of the diamond lattice was turned into graphite its density became lower resulting in a volume increase and internal strain formation. Diamond graphitization is a process that takes few picoseconds 15 . If the laser pulse length is equal or shorter than few picoseconds, it has been shown 15,16,17 that when it ends, a graphite nucleus of minimal dimension is formed without damage in the crystal lattice due to after-pulse lattice heating. Here nanosecond pulses were used. As a consequence after the formation of the initial graphite nucleus a large amount of energy was still deposited leading to the graphite nucleus overgrowth. It resulted in the creation of highly strained areas where macroscopic cracks occurred because the local strain overcame the diamond tensile strength (2.8 GPa 18 ). Using the setup shown in Fig. 1(a), graphitic channels (500µm long, 30µm diameter) were fabricated in a 10x10x0.7mm 3 poly-crystalline Chemical Vapor Deposition (pc-CVD) diamond purchased by Element 6 Ltd. The diamond had two active areas: the first one with graphitic buried electrodes, and the second one with conventional planar electrodes (Fig. 2.).

The device with planar structure was fabricated as a reference to evaluate the CCE on the same diamond sample and quantify any improvement of the detector performances. The buried structure enables to assess that graphitic columns are electrically active and suitable as electrodes for radiation detector applications. Dark currents of both structures were measured (Fig. 3) to assess their suitability for radiation detection application. The hybrid device exhibits a leakage current two orders of magnitude above that of the planar device. In fact, this directly results from the shorter distance between graphitic electrodes and the aluminum back contact (200µm) leading to locally enhanced electric fields in the sample. According to electric field simulations achieved at +500V (Fig. 2.), buried electrodes generate a tip termination effect that locally enhances the electric field up to 4V/µm. However, despite these high electric fields, the leakage current remained below 50pA thus making this configuration suitable for radiation detection. Devices were subsequently tested using α-particles through a collimated 241-Americium source (5,486 MeV α-particles). CCE was measured using a fast charge sensitive amplifier [START_REF] Ciobanu | IEEE Nuclear Science Symposium Conference Record[END_REF] . The signal was recorded using a Lecroy Wavepro 725Zi digital oscilloscope.

Calibration of the acquisition chain was performed using the same source and a 500µm thick CVD single crystal (Element 6 Ltd.) biased at +500V exhibiting a CCE close to 100% as we already showed in a previous communication [START_REF] Pomorski | Electronic Properties of Single Crystal CVD Diamond and Its Suitability for Particle Detection in Hadron Physics Experiments[END_REF] . Fig. 4 shows spectra of both devices in top and back injection configurations. As no resolved peak can be observed, average efficiencies were calculated as the weighted means of each spectrum. In both configurations, the device with planar electrodes exhibited a CCE of 15% and values up to 50% could be noticed locally most probably correlated with large grains with lower grain boundaries density. CCE performances of the hybrid device change depending on the side where α-particles were injected. When α-particles were injected from the top side, CCE of only 6% was observed compared to 15% in planar configuration. This behavior was expected from the electric field distribution shown in Fig. 2. Two different regimes must be considered in this case : firstelectrons are created at the vicinity of the top electrode where the local electric field is very low (0.15V/µm with respect to 0.75V/µm in planar configuration). To contribute to the detection signal, they have to drift in a low gradient electric field along the graphitic columns, during several hundreds of µm before they only enter the high field region for at most 200µm. During the long drift path at low field, the drift distance being low in this kind of pc-CVD material, several electrons will not manage to reach the high electric field region thus nor contribute to the signal, thus explaining the lower detection efficiencies observed when particles are injected from the top side.

On the other hand, when α-particles were injected from the back side, electrons were created in a high electric field region where their drift path is short (only 200µm), and thus more likely to efficiently contribute to the signal. As a result, the hybrid device exhibited an average efficiency of 20% with maximum values reaching 75%. In this configuration, the hybrid device CCE was 1.5 times larger than the planar one. Data extracted from spectra are consistent with the electric field simulations and show a slight but real improvement in the charge collection properties of the hybrid device compared to planar. However the drift path was reduced by a factor 3.5, so a better improvement of the CCE was expected. Damages and strain created during the electrodes fabrication might be responsible of this by increasing the trapping probability. This constitutes the major drawback of the technique. But it has been shown that using a picosecond [START_REF] Shimizu | [END_REF] or femtosecond 22 laser can lead to the production of graphitic wires with minimal surrounding defects.

In conclusion, we reported on the fabrication of 3D graphitic electrodes buried in diamond using laser micro-machining. Our results demonstrate that these electrodes are electrically active and allow charge collection. As a consequence, they appear suitable for future radiation detector applications. However, the electrodes fabrication process still need to be optimized to achieve better CCE. For this, the use of other lasers sources is under investigation.

Nevertheless, this work constitutes a major progress toward the fabrication of a 3D diamond detector which could be a significant improvement for radiation hardness issues in high energy physics applications.

FIG. 1 .

 1 FIG. 1. (Color Online) (a) Laser setup schematic. (b) Optical microscope picture of three graphitic columns with various depths.

FIG. 2 .

 2 FIG. 2. (Color online) Electric field simulations [COMSOL Multiphysics] in the hybrid and planar geometries. Electric field color map is in logarithmic scale. (a) Hybrid geometry: an array of 5 by 5 channels (30µm diameter, 500µm depth) with 150µm pitch. The channels are connected all together with an aluminum pad on the top surface. On the back surface, an aluminum pad has been evaporated too. (b) Planar geometry: two aluminum contacts evaporated on the top and back surfaces of the pc-CVD diamond sample.

FIG. 3 .

 3 FIG. 3. Dark current measurements within both active areas. Planar detector's characteristic is plotted in filled grey triangle. Hybrid detector's I(V) curve is plotted in filled black square.

FIG. 4 .

 4 FIG. 4. CCE comparison of both devices. α-particles impinging the device from the top (a) or the back (b) sides.

ACKNOWLEDGEMENTS

The authors would like to thank Alex ander Oh from The University of Manchester for providing us the diamond sample used in this work.

REFERENCES:


Suelen Gasparin

Denys Dutykh

Nathan Mendes

A spectral method for solving heat and moisture transfer through consolidated porous media

Keywords: Spectral methods, Chebyshev polynomials, Tau-Galerkin method, numerical simulation, heat and moisture transfer, reduced-order modeling primary), 35K05, 80A20, 65M32 (secondary) 02.70.Bf (secondary)

This work presents an efficient numerical method based on spectral expansions for simulation of heat and moisture diffusive transfers through multilayered porous materials. Traditionally, by using the finite-difference approach, the problem is discretized in time and space domains (Method of lines) to obtain a large system of coupled Ordinary Differential Equations (ODEs), which is computationally expensive. To avoid such a cost, this paper proposes a reduced-order method that is faster and accurate, using a much smaller system of ODEs. To demonstrate the benefits of this approach, tree case studies are presented. The first one considers nonlinear heat and moisture transfer through one material layer. The second case -highly nonlinear -imposes a high moisture content gradient -simulating a rain like condition -over a two-layered domain, while the last one compares the numerical prediction against experimental data for validation purposes. Results show how the nonlinearities and the interface between materials are easily and naturally treated with the spectral reduced-order method. Concerning the reliability part, predictions show a good agreement with experimental results, which confirm robustness, calculation efficiency and high accuracy of the proposed approach for predicting the coupled heat and moisture transfer through porous materials.

Introduction

Energy consumption of conditioned spaces is strongly related to heat transfer through porous structures, which are dependent of external (weather) and internal conditions. Moreover, moisture migration and accumulation within the porous structures considerably affect the transient hygrothermal performance of porous elements, especially in buildings. Thus, coupled hygrothermal simulations are important to estimate envelope energy performance and risks associated to the presence of moisture such as material degradation, mold growth and related health aspects of occupants.

Initially, the major hypothesis was of no coupling between governing equations of heat and moisture transfers. Although, in porous materials, moisture transfer and accumulation have a direct impact on the heat transfer, especially when there is phase change [5, [START_REF] Mulay | Simultaneous transport of heat and moisture in a partially saturated porous media[END_REF]43]. As presented by Deru [13], to precisely determine the heat losses, simultaneous calculations with moisture content are required, as they are closely interdependent. Temperature and moisture contents are highly interconnected phenomena and, hence, must be simulated together [START_REF] Berger | Factors governing the development of moisture disorders for integration into building performance simulation[END_REF]33]. For instance, moisture can affect the effective thermal conductivity by a factor of ten.

However, there are some challenges on how to characterize mathematically those physical processes, due to the complexity of the physical phenomena and to the porous structure. Studies of heat and moisture transfer have been published since 1950's by Philip and De Vries [START_REF] Philip | Moisture movement in porous materials under temperature gradients[END_REF] and Luikov [29]. They represented the coupled processes of heat and moisture transfer by a system of two nonlinear second-order partial differential equations which uses as driving potentials the temperature and the moisture content gradients. The system is nonlinear mainly since of the phase change term in the energy conservation equation and also due to the fact the transport coefficients are highly moisture content and/or temperature dependent. Although, the paramount importance of accurately predict moisture content and temperature fields in several fields of science, detailed simulation of the heat and moisture transfer has only been performed with the improvement of the computer systems. The main reason for that is due to the complexity of the problem and to the high computational cost to solve the coupled heat and moisture transfer equations, which is at least two orders of magnitude greater than that needed to solve only the heat conduction equation [13].

The numerical models used to predict heat and moisture transfer have to deal with multiscale effects in both space and time domains, with different characteristic times and lengths. For example, simulation of building behavior is generally analyzed on the scale of one year (or more). However, the phenomena and particularly the boundary conditions evolve on a time scale of minutes or even seconds. The geometric configurations of the buildings require three-dimensional modeling of lengthy elements such as facades and ground. Furthermore, when dealing with heat and moisture, the nonlinear behavior of the materials should be taken into account. The combination of all those factors implies the use of more robust and efficient numerical methods since no analytical solution does exist for those problems.

For modeling purposes, the problem previously described is generally solved by the finitedifference methods [START_REF] Fan | An improved model of heat and moisture transfer with phase change and mobile condensates in fibrous insulation and comparison with experimental results[END_REF]19], the finite-volume methods [14, [START_REF] Manz | Experimental and numerical study of a mechanically ventilated glass double façade with integrated shading device[END_REF][START_REF] Mendes | A method for predicting heat and moisture transfer through multilayered walls based on temperature and moisture content gradients[END_REF] and the finite-element methods [27, [START_REF] Rouchier | Influence of concrete fracture on the rain infiltration and thermal performance of building facades[END_REF][START_REF] Thomas | A fully nonlinear analysis of heat and mass transfer problems in porous bodies[END_REF], which are well-established methods in the field of heat and mass transfer. The basics of a numerical approach is the idea of making an approximation of the solution, which takes a solution for a finite number of degrees of freedom (DOFs). The greater the number of discrete points, the closer to the exact solution will be the approximated solution [31].

However, when the solution of one problem requires long simulation periods (years), considering an entire building, with a very fine time step and mesh refinement, computation becomes too time-consuming. To solve this problem, model reduction techniques can be used as an alternative to approach the solution of the problem and reduce the computational cost.

The intent to construct reduced-order models (ROMs) is to provide an accurate description of the physical phenomena by decreasing the number of degrees of freedom, while retaining the model's fidelity, at a computational cost much lower than the large original model [START_REF] Reddy | Constrained reduced-order models based on proper orthogonal decomposition[END_REF]. In recent years, reduced-order modeling techniques have proven to be powerful tools for solving various problems. Important efforts have been dedicated to developing reduced-order models that can provide accurate predictions while dramatically reducing computational time, for a wide range of applications, covering different fields such as fluid mechanics, heat transfer, structural dynamics, among others [START_REF] Bai | Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems[END_REF][START_REF] Herzet | Model reduction from partial observations[END_REF][START_REF] Lucia | Reduced-order modeling: new approaches for computational physics[END_REF]. Reduced-order models -such as POD (Proper Orthogonal Decomposition), MBR (Modal Basis Reduction) and PGD (Proper Generalized Decomposition) -have shown a relevant reduction of the computational cost and have been successfully employed by the building physics community [3]. In those works, they have applied reduced-order models to build accurate solutions with less computational effort than the complete original model. Reduced-order models can be classified as a priori or a posteriori methods. The a posteriori approaches need a preliminary computed (or even experimental) solution data of the large original problem to build the reduced one. Whereas the a priori ones do not need preliminary information on the studied problem. The reduced-order model is unknown a priori and is directly built. A careful attention must be paid regarding the definition of ROMs since sometimes it is related to the degradation of the physical model [46], which is not the case of the present work.

Another promising approach to solve the coupled transfer problem is the spectral method, which is a robust and highly accurate method that has been applied to solve partial differential equations since the 70's, but lost its spot due to the difficulties to treat nonlinearities, complex geometries, irregular domains, and non-periodic boundary conditions. However, the Spectral methods have overcome some of the mentioned difficulties and now they are successfully applied in studies of wave propagation, meteorology, computational fluid dynamics, quantum mechanics and other fields [START_REF] Canuto | Spectral Methods: Fundamentals in Single Domains[END_REF]. The main attractiveness of this method is the superior rate of convergence and the low dissipation and dispersion errors, making its use also attractive to the industry. Nonetheless, spectral methods still have some constraints. For example, complex geometries are one of their main drawbacks as they work better when the geometry of the problem is fairly smooth and regular [START_REF] Gasparin | UMR 5127[END_REF], while finite-element methods are particularly well suited to problems in very complex geometries. Moreover, spectral methods can offer higher accuracy in geometries like boxes and spheres, which can be combined into more complex shapes [START_REF] Canuto | Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics[END_REF]17]. In applications where geometry-related disadvantages are not present, the classic finite-element, finite-volume and finite-difference methods do not come close in terms of efficiency.

Some works related to transport phenomena can be found in literature involving diffusive [START_REF] Guo | The Application of the Chebyshev-Spectral Method in Transport Phenomena[END_REF][START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF], convective [START_REF] Chen | Spectral collocation method for natural convection in a square porous cavity with local thermal equilibrium and non-equilibrium models[END_REF]41] and radiative [START_REF] Chen | Chebyshev collocation spectral method for solving radiative transfer with the modified discrete ordinates formulations[END_REF][START_REF] Li | Iterative and direct Chebyshev collocation spectral methods for one-dimensional radiative heat transfer[END_REF][START_REF] Ma | Thermal radiation heat transfer in one-and two-dimensional enclosures using the spectral collocation method with full spectrum kdistribution model[END_REF] heat transfer. The spectral techniques applied in these works are diversified, adopted according to the geometry, boundary conditions and field of application. In recent works, researchers have implemented spectral methods for solving heat and moisture transfer in food engineering [38] and on fluid flow [36]. Recently, in [20], the authors have studied the moisture transfer in porous building materials considering layered domains, and in [START_REF] Gasparin | Advanced Reduced-Order Models for Moisture Diffusion in Porous Media[END_REF], they have compared the Spectral method to others ROMs, applied to parametric problems of the building physics field.

Therefore, the scope of this work is to continue the investigations presented in [20] and [START_REF] Gasparin | Advanced Reduced-Order Models for Moisture Diffusion in Porous Media[END_REF], extending it to the coupled heat and mass transfer. Here, the Spectral method is used to compute one-dimensional heat and moisture diffusion transfer in porous materials, which is validated against experimental data from the literature. The problems treated here involve irregular domains, high nonlinear transport coefficients and non-periodic boundary conditions. The objective is to significantly reduce the computational cost while maintaining high fidelity solutions. This technique assumes separated tensorial representation of the solution by a finite sum of function products. It fixes a set of spatial basis functions to be the Chebyshev polynomials and then, a system of ordinary differential equations is built to compute the temporal coefficients of the solution using the Tau-Galerkin method.

The efficiency of the Spectral approach will be analyzed for simple and multilayered domains with highly nonlinear properties and with sharp boundary conditions and profiles of solutions. For this purpose, the manuscript is organized as follows. First, the description of the physical phenomena is presented (Section 2). Then, the Spectral technique is described (Section 3). In the sequence, the proposed method is applied to two different cases: (i) considering heat and moisture transfer through a single layer (Section 4.1) and (ii) focusing on the heat and moisture transfer through a multilayered domain (Section 4.2). Finally, the method and the model are compared with experimental data from literature in Section 5, considering a single material with real boundary conditions. The main conclusions of the study are outlined in Section 6.

Mathematical model

The physical problem considers heat and moisture transfer through a porous material defined in the one-dimensional spatial domain Ω x = [ 0, L ] and in the time horizon Ω t = [ 0, τ ] . The moisture transfer occurs due to capillary migration and vapour diffusion.

The heat transfer is governed by diffusion and latent mechanisms. The physical model of the problem can be formulated as [3]:

1 P s ∂w ∂φ ∂P v ∂t - ∂ ∂x k l R v T ρ l P v + δ v ∂P v ∂x = 0 , (2.1a)
ρ 0 c 0 + w c w ∂T ∂t - ∂ ∂x λ ∇T + L v δ v ∂P v ∂x = 0 , (2.1b)
where w [kg/m 3 ] is the material volumetric moisture content, φ [-], the relative humidity, δ v [s] and k l [s], the vapour and liquid permeabilities, P v [Pa], the vapour pressure, P s [Pa], the saturation pressure, T [K], the temperature, R v [J/(kg • K)] , the water vapour gas constant, c 0 [J/(kg • K)], the material specific heat, ρ 0 [kg/m 3 ], the material specific mass, ρ l [kg/m 3 ] , the water specific mass, c w [J/(kg • K)] , the water specific heat, λ [W/(m • K)], the thermal conductivity, and, L v [J/kg] , the latent heat of vaporization. Table 1 presents the values of the water properties considered in this work.

The relation between the moisture content w and the relative humidity φ is given by the sorption isotherm and the relation between the vapour pressure P v and the relative humidity φ is given by φ = P v /P s(T ) . In addition, the following assumptions are adopted in this study: (i) no hysteresis effect; (ii) temperature within the range [0, 40] • C ; (iii) no temperature dependency on the mass balance equation and (iv) properties are dependent only on the vapour pressure field.

Thus, considering the following notation: the energy storage coefficient, system (2.1) can be rewritten in one-dimensional form as:

k M def := k l ρ l R v T P v + δ v :
c M (P v ) ∂P v ∂t - ∂ ∂x k M (P v ) ∂P v ∂x = 0 , (2.2a) c T (P v ) ∂T ∂t - ∂ ∂x k T (P v ) ∂T ∂x + k T M (P v ) ∂P v ∂x = 0 . (2.2b)
Finally, the problem of interest is a coupled system of two nonlinear parabolic partial differential equations, with vapour pressure P v and temperature T gradients as driving potentials.

Boundary conditions. The moisture exchange between the environment and the surface is driven by vapour exchange, including evaporation and condensation, and by driving rain: The heat balance at the boundary includes the convective exchange, the latent heat transfer due to vapour exchange, and the sensible heat transfer due to precipitation, which is expressed as:

n • k M ∂P v ∂x = h M P v -P v, ∞ (t) -g ∞ (t) , (2.3
n • k T ∂T ∂x + k T M ∂P v ∂x = h T T -T ∞ (t) + L v h M P v -P v, ∞ (t) -H l g ∞ (t) , (2.4)
where T ∞ [K] is the temperature of the air that varies over time and h T [W/(m Initial conditions. The initial conditions can either have a uniform distribution or a profile more appropriated to the boundary conditions to reduce a warm-up simulation period, which can be very significant depending on the material hygrothermal properties and on the thickness of the building component.

P v (x , t = 0) = P v, 0 (x) , T (x , t = 0) = T 0 (x) .
Interface. The configuration assumed at the interface between materials follows the hydraulic continuity [START_REF] De Freitas | Moisture migration in building walls -Analysis of the interface phenomena[END_REF], which considers interpenetration of both porous structure layers. Consider two different materials, both of them are homogeneous and isotropic, and the coupled heat and moisture transfer are simulated, through a perfectly airtight structure.

The hydraulic continuity assumes that there is a continuous moisture distribution of vapour content and temperature:

P v, 1 (x int , t) = P v, 2 (x int , t) , (2.5a
) [Pa] a continuous moisture flow and a continuous heat flux across the interface verify:

T 1 (x int , t) = T 2 (x int , t) , (2.5b)
k M, 1 ∂P v, 1 ∂x x int = k M, 2 ∂P v, 2 ∂x x int , (2.6a)
k T, 1 ∂T 1 ∂x + k T M, 1 ∂P v, 1 ∂x x int = k T, 2 ∂T 2 ∂x + k T M, 2 ∂P v, 2 ∂x x int , (2.6b)
where x int ∈ Ω x represents the location of the interface between materials. Therefore, we can split the spatial domain in two parts, Ω x, 1 = [ 0 , x int ] and Ω x, 2 = ( x int , L ] , which represent the spatial domain of material 1 and material 2 .

Fluxes and flows. One of the interesting outputs in the building physics framework is the heat flux, divided into sensible q s and latent q l heat fluxes [W/m 2 ], which are defined as:

q s (t) def := -k T ∂T ∂x x 0 and q l (t) def := -k T M ∂P v ∂x x 0 .
The moisture flow g [kg/(m 2 • s)] is similarly computed:

g (t) def := -k M ∂P v ∂x x 0 , where x 0 ∈ [ 0 , L ] .

Dimensionless formulation

Before solving directly the problem, it is of capital importance to get a dimensionless formulation of the problem under consideration [25]. In this way, we define the following dimensionless quantities:

u def := T T ref , v def := P v P v, ref , x ⋆ def := x L ref , t ⋆ def := t t ref , c ⋆ M def := c M • L 2 ref k M, ref • t ref , c ⋆ T def := c T • L 2 ref k T, ref • t ref , k ⋆ M def := k M k M, ref , k ⋆ T def := k T k T, ref , k ⋆ T M def := k T M • P v, ref k T, ref • T ref , Bi T def := h T • L ref k T, ref , Bi M def := h M • L ref k M, ref , Bi T M def := h M • L v • L ref • P v, ref k T, ref • T ref , g ⋆ ∞ def := g ∞ • L ref k M, ref • P v, ref , H ⋆ l def := H l • k M, ref • P v, ref k T, ref • T ref .
where the subscript ref represents a reference value, chosen according to the application problem and the superscript ⋆ represents a dimensionless quantity of the same variable. Therefore, the governing system (2.2) can be written in a dimensionless form as:

c ⋆ M ∂v ∂t ⋆ - ∂ ∂x ⋆ k ⋆ M ∂v ∂x ⋆ = 0 , (2.7a) c ⋆ T ∂u ∂t ⋆ - ∂ ∂x ⋆ k ⋆ T ∂u ∂x ⋆ + k ⋆ T M ∂v ∂x ⋆ = 0 . (2.7b)
The dimensionless formulation of the boundary conditions are:

n • k ⋆ M ∂v ∂x ⋆ = Bi M v -v ∞ (t ⋆ ) -g ⋆ ∞ (t ⋆ ) , n • k ⋆ T ∂u ∂x ⋆ + k ⋆ T M ∂v ∂x ⋆ = Bi T u -u ∞ (t ⋆ ) + Bi T M v -v ∞ (t ⋆ ) -H ⋆ l g ⋆ ∞ (t ⋆ ) ,
and of the initial conditions are:

u (x ⋆ , t ⋆ = 0) = u 0 (x ⋆ ) , v (x ⋆ , t ⋆ = 0) = v 0 (x ⋆ ) .
Interface conditions in the dimensionless form are written as:

v 1 (x ⋆ int , t ⋆ ) = v 2 (x ⋆ int , t ⋆ ) , (2.8a
)

u 1 (x ⋆ int , t ⋆ ) = u 2 (x ⋆ int , t ⋆ ) , (2.8b
)

k ⋆ M, 1 ∂v 1 ∂x ⋆ x ⋆ int = k ⋆ M, 2 ∂v 2 ∂x ⋆ x ⋆ int , (2.8c
)

k ⋆ T, 1 ∂u 1 ∂x ⋆ + k ⋆ T M, 1 ∂v 1 ∂x ⋆ x ⋆ int = k ⋆ T, 2 ∂u 2 ∂x ⋆ + k ⋆ T M, 2 ∂v 2 ∂x ⋆ x ⋆ int .
(2.8d)

In the following, we drop ⋆ for the sake of clarity.

Spectral reduced-order model

Spectral methods consider a sum of polynomials that suit for the whole domain, providing a high approximation of the solution. The smoother a function is, the faster the convergence of its spectral series [START_REF] Gasparin | UMR 5127[END_REF]. For considerably smooth problems, the error decreases exponentially, making the solution with the same order of accuracy of other methods but with a much lower number of degrees of freedom. As a result, this method has a low memory usage, allowing to store and operate a lower number of variables [48]. The Spectral methods used in this work are the Chebyshev polynomials on the basis function and the Tau-Galerkin method to compute the temporal coefficients.

Method description

Problem (2.7) has an important difficulty in dealing with the nonlinearities of the storage c and diffusion k coefficients, all of them depending on the moisture content field. These coefficients are usually given by empirical functions from experimental data. For this reason, Eq. (2.7a) and (2.7b) are recalled with a simplified notation:

c M (v) ∂v ∂t - ∂ ∂x k M (v) ∂v ∂x = 0 , (3.1a) c T (v) ∂u ∂t - ∂ ∂x k T (v) ∂u ∂x + k T M (v) ∂v ∂x = 0 . (3.1b)
In addition, boundary conditions are simplified and also written with a shorter notation, just for the sake of explaining the method in a pedagogical way/manner:

n • k M (v) ∂v ∂x = Bi M, L v -v ∞, L (t) , (3.2a) n • k T (v) ∂u ∂x + k T M (v) ∂v ∂x = Bi T, R u -u ∞, R (t) (3.2b) + Bi T M, R v -v ∞, R (t) . (3.2c)
A special attention must be given to the spatial domain because the Chebyshev Spectral method is traditionally presented on the canonical interval -1 , 1 . Thus, if the dimensionless interval is not within -1 , 1 , a change of variables (domain transformation) must be performed for the computational domain.

In order to apply better the spectral method, Eqs. (3.1a) and (3.1b) are written in the non-conservative form as:

∂v ∂t -ν (v) ∂ 2 v ∂x 2 -λ (v) ∂v ∂x = 0 , (3.3a) ∂u ∂t -α (v) ∂ 2 u ∂x 2 -β (v) ∂u ∂x -γ (v) ∂ 2 v ∂x 2 -δ (v) ∂v ∂x = 0 , (3.3b)
where,

ν (v) def := k M (v) c M (v) , λ (v)
def := 1 c M (v) • ∂ k M (v) ∂x , α (v) def := k T (v) c T (v) , β (v)
def := 1 c T (v) • ∂ k T (v) ∂x , γ (v) def := k T M (v) c T (v) , δ (v)
def := 1 c T (v) • ∂ k T M (v) ∂x .
The unknowns u ( x, t ) and v ( x, t ) from Eq. (3.3) are accurately represented as a finite sum [START_REF] Mendes | Numerical methods for diffusion phenomena in building physics[END_REF]Chap. 6]:

v ( x, t ) ≈ v n ( x, t ) = n i = 0 a i (t) T i (x) , i = 0, 1, 2, . . . , n , (3.4a) u ( x, t ) ≈ u n ( x, t ) = n i = 0 b i (t) T i (x) , i = 0, 1, 2, . . . , n . (3.4b)
Here, {T i (x)} n i = 0 are the Chebyshev polynomials, {a i (t)} n i = 0 are the corresponding time-dependent spectral coefficients and n represents the number of degrees of freedom of the solution component. Eqs. (3.4a) and (3.4b) can be seen as a series truncation from N = n + 1 modes. The Chebyshev polynomials are chosen as the basis functions since they are optimal in L ∞ approximation norm [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]. Therefore, the expression of the derivatives in the Chebyshev basis are:

∂v n ∂x = n i = 0 a i (t) ∂T i ∂x (x) = n i = 0 ã i (t) T i (x) , (3.5a) ∂ 2 v n ∂x 2 = n i = 0 a i (t) ∂ 2 T i ∂x 2 (x) = n i = 0 ã i (t) T i (x) , (3.5b
)

∂v n ∂t = n i = 0 ȧ i (t) T i (x) , (3.5c
ã i = 2 c i n-1 p = i + 1 p + i odd p a p , i = 0, 1, . . . , n -1, ã n ≡ 0 , ã i = 1 c i n-2 p = i + 2 p + i even p p 2 -i 2 a p , i = 0, 1, . . . , n -2, ã n-1 ≡ ã n ≡ 0 ,
with,

c i = 2 , if i = 0 , 1 , if i > 0 .
The derivatives of u (x, t) are written in a similar way by just replacing v by u and a by b . Thus, the derivatives defined by Eqs. (3.5a), (3.5b) and (3.5c) are replaced into Eqs. (3.3a) and (3.3b) to provide the residuals:

R 1 = n i = 0 ȧ i (t) T i (x) -ν n i = 0 a i (t) T i (x) n i = 0 ã i (t) T i (x) -λ n i = 0 a i (t) T i (x) n i = 0 ã i (t) T i (x) , (3.6a) R 2 = n i = 0 ḃ i (t) T i (x) -α n i = 0 a i (t) T i (x) n i = 0 b i (t) T i (x) -β n i = 0 a i (t) T i (x) n i = 0 b i (t) T i (x) -γ n i = 0 a i (t) T i (x) n i = 0 ã i (t) T i (x) -δ n i = 0 a i (t) T i (x) n i = 0 ã i (t) T i (x) , (3.6b)
which are considered a misfit of the approximate solution. The purpose is to minimize the residual so the solution satisfies the governing equations. To this end, the residual is minimized via the Tau-Galerkin method, which requires that Eqs. (3.6a) and (3.6b) be orthogonal to the Chebyshev basis functions R , T j = 0 :

R , T j = ˆ1 -1 R ( x , t ) T j ( x ) √ 1 -x 2 dx = 0 , j = 0, 1, 2, . . . , n -2 . (3.7)
As a result, the project residuals are:

M • ȧ i (t) = G i, j • ã i (t) + Λ i, j • ã i (t) , (3.8a) M • ḃ i (t) = M i, j • b i (t) + N i, j • b i (t) + F i, j • ã i (t) + J i, j • ã i (t) , (3.8b)
where, M is a diagonal and the singular matrix (rank ( M ) = N -2) which contains the coefficients of the Chebyshev weighted orthogonal system. The matrix M has the following form:

M =           π π 2 0 . . . π 2 0 0 0          
, and, the matrices with indices (i , j) written as:

G i, j = ˆ1 -1 ν n i=0 a i (t) T i (x) T i (x) T j (x) √ 1 -x 2 dx , Λ i, j = ˆ1 -1 λ n i=0 a i (t) T i (x) T i (x) T j (x) √ 1 -x 2 dx , M i, j = ˆ1 -1 α n i=0 a i (t) T i (x) T i (x) T j (x) √ 1 -x 2 dx , N i, j = ˆ1 -1 β n i=0 a i (t) T i (x) T i (x) T j (x) √ 1 -x 2 dx , F i, j = ˆ1 -1 γ n i=0 a i (t) T i (x) T i (x) T j (x) √ 1 -x 2 dx , J i, j = ˆ1 -1 δ n i=0 a i (t) T i (x) T i (x) T j (x) √ 1 -x 2 dx .
with indices i , j are the ones defined in Eqs. (3.4) and (3.7). By using the Chebyshev-Gauß quadrature, the integrals are also approximated by a finite sum:

G i, j ≈ π m m k = 1 ν k T i (x k ) T j (x k ) , Λ i, j ≈ π m m k = 1 λ k T i (x k ) T j (x k ) , M i, j ≈ π m m k = 1 α k T i (x k ) T j (x k ) , N i, j ≈ π m m k = 1 β k T i (x k ) T j (x k ) , F i, j ≈ π m m k = 1 γ k T i (x k ) T j (x k ) , J i, j ≈ π m m k = 1 δ k T i (x k ) T j (x k ) ,
where,

ν k def := ν n i = 0 a i (t) T i (x k ) , λ k def := λ n i = 0 a i (t) T i (x k ) , α k def := α n i = 0 a i (t) T i (x k ) , β k def := β n i = 0 a i (t) T i (x k ) , γ k def := γ n i = 0 a i (t) T i (x k ) , δ k def := δ n i = 0 a i (t) T i (x k ) ,
and x k are the Chebyshev nodes:

x k = cos 2 k -1 2 m π , k = 1, 2, . . . , m .
The value of m is approximately the same as the number of modes, as discussed in [20].

To complete the problem, the boundary conditions are also written in the form of the residuals. First, the boundary conditions associated to the moisture transport -Eq. (3.2a):

ω 1 = k M n i = 0 a i (t) T i (-1) n i = 0 ã i (t) T i (-1) -Bi M, L n i = 0 a i (t) T i (-1) -v ∞, L , ω 2 = k M n i = 0 a i (t) T i (1) n i = 0 ã i (t) T i (1) + Bi M, R n i = 0 a i (t) T i (1) -v ∞, R .
Then, the boundary conditions for Eq. (3.2c), regarding the heat transport:

κ 1 = k T n i = 0 a i (t) T i (-1) n i = 0 b i (t) T i (-1) + k T M n i = 0 a i (t) T i (-1) n i = 0 ã i (t) T i (-1) -Bi T, L n i = 0 b i (t) T i (-1) -u ∞, L -Bi T M, L n i = 0 a i (t) T i (-1) -v ∞, L , κ 2 = k T n i = 0 a i (t) T i (1) n i = 0 b i (t) T i (1) + k T M n i = 0 a i (t) T i (1) n i = 0 ã i (t) T i (1) + Bi T, R n i = 0 b i (t) T i (1) -u ∞, R + Bi T M, R n i = 0 a i (t) T i (1) -v ∞, R ,
with T i ( -1) = ( -1) i and T i ( 1 ) ≡ 1 (see [39] for more details).

In this way, it is possible to compose the system of ODEs to be solved, plus the four additional algebraic expressions regarding the boundary conditions. Finally, the system of differential-algebraic equations (DAEs) has the following form:

M 0 0 M ȧ n ḃ n = A 0 B C • a n (t) b n (t) + b 1 (t) b 2 (t) , (3.9)
where b 1 (t) and b 2 (t) are vectors containing the boundary conditions, previously defined by ω 1 , ω 2 , κ 1 and κ 2 :

b 1 (t) =           0 0 . . . 0 ω 1 ω 2           and b 2 (t) =           0 0 . . . 0 κ 1 κ 2           .
Matrix A is written from the right member of Eq. (3.8a), and, matrices B and C are written from the right member of Eq. (3.8b). Initial values of the coefficients {a i (t = 0)} and {b i (t = 0)} are calculated by the orthogonal projection of the initial condition [START_REF] Canuto | Spectral Methods: Fundamentals in Single Domains[END_REF]:

a 0, i ≡ a i (0) = 2 π c i ˆ1 -1 v 0 (x) T i (x) √ 1 -x 2 dx , i = 0, 1, . . . , n , (3.10a) b 0, i ≡ b i (0) = 2 π c i ˆ1 -1 u 0 (x) T i (x) √ 1 -x 2 dx , i = 0, 1, . . . , n , (3.10b)
where, v 0 (x) and u 0 (x), are the dimensionless initial condition. Therefore, the reduced system of ODEs composed from Eqs. (3.9) and (3.10) can be solved. Different approaches can be used to solve the system of ODEs (3.9). The most straightforward solution is to apply a numerical integration scheme, with moderate accuracy. So, with an embedded error control and not so stringent tolerances, it can be done very efficiently. In this work, the Matlab ™ environment was used to perform simulations, and the solvers ODE15s or ODE23t were used to solve the differential-algebraic system of equations (DAEs). The output are the vectors of spectral coefficients {a i (t)} n i = 0 and {b i (t)} n i = 0 . Then, it enables to reconstruct the solution thanks to spectral representations.

Extension to multilayered domains

We present the idea for two materials, but it can be generalized to any finite number of subdomains straightforwardly. Consider that the original spatial domain

Ω x = [ 0, L ] is decomposed into two subdomains Ω x, 1 = [ 0, x int ]
and Ω x, 2 = ( x int , L ] , which represent each material surface and x int represents the location of the interface, as previously defined. From this, the unknowns v (x, t) and u (x, t) are then written as: for t 0 . Thus, v 1 , v 2 , u 1 and u 2 are written respectively as:

v (x, t) = v 1 (x, t) , x ∈ Ω x, 1 , v 2 (x, t), x ∈ Ω x, 2 , u (x, t) = u 1 (x, t), x ∈ Ω x, 1 , u 2 (x, t), x ∈ Ω x, 2 ,
v 1 (x, t) = n i = 0 a i, 1 (t) T i (x) , v 2 (x, t) = n i = 0 a i, 2 (t) T i (x) , u 1 (x, t) = n i = 0 b i, 1 (t) T i (x) , u 2 (x, t) = n i = 0 b i, 2 (t) T i (x) .
Note that the Chebyshev polynomials T (x) are always the same and x must always be in the closed set [ -1, 1 ] . To assure this, the subdomains Ω x, 1 and Ω x, 2 are linearly transformed to the spectral domains Ω x, 1 and Ω x, 2 , as illustrated in Figure 1, so they can fit within the interval of interest, which should be seen as a horizontal coordinate in each material. In addition, by separating the real domain in two spectral domains, the smoothness of the solution is assured. The condition at the interface between the two materials states the continuity of the fields and the flows as in Eq. (2.8). It implies that the derivatives of the fields u and v are not continuous at the interface. This important remark has to be taken into account in the construction of the spectral reduced-order model. Indeed, the domain is decomposed in subdomains to maintain a smooth solution and particularly a continuous derivative on each sub-domain. In this way, the model order reduction is optimal and ensures the error of the Spectral-ROM to decrease exponentially. It is possible to build the reduced order model considering the whole domain (without decomposition). However, the convergence is undermined since the solution and its derivatives are not smooth at the interface between two materials. Many more modes would be necessary to reach the same accuracy, as detailed in Theorem 1 of [START_REF] Trefethen | Spectral methods in MatLab[END_REF]Chap. 4].

By considering the two materials, the DAE system represented in Eq. (3.9) becomes:

     M 0 M M 0 M      •         ȧ i, 1 ḃ i, 1 ȧ i, 2 ḃ i, 2         =          A 1 0 0 B 1 C 1 A 2 0 0 B 2 C 2          •         a n, 1 (t) b n, 1 (t) a n, 2 (t) b n, 2 (t)         +         b 1 (t) b 2 (t) b 3 (t) b 4 (t)         .
The interface conditions -Eq. (2.8) -are then written in the spectral form as:

ϑ 1 = n i = 0 a i, 1 (t) - n i = 0 a i, 2 (t) (-1) i , ϑ 2 = k M, 1 n i = 0 a i, 1 (t) n i = 0 ã i, 1 (t) -k M, 2 n i = 0 a i, 2 (t) (-1) i n i = 0 ã i, 2 (t) (-1) i , ϑ 3 = n i = 0 b i, 1 (t) - n i = 0 b i, 2 (t) (-1) i , ϑ 4 = k T, 1 n i = 0 a i, 1 (t) n i = 0 b i, 1 (t) + k T M, 1 n i = 0 a i, 1 (t) n i = 0 ã i, 1 (t) -k T, 2 n i = 0 a i, 2 (t) (-1) i n i = 0 b i, 2 (t) (-1) i -k T M, 2 n i = 0 a i, 2 (t) (-1) i n i = 0 ã i, 2 (t) (-1) i ,
which are included in vectors b 1 and b 2 and set equal to zero. In the same way, the boundary conditions are written in the spectral form as:

σ 1 = k M, 1 n i = 0 a i (t) T i (-1) n i = 0 ã i (t) T i (-1) -Bi M, L n i = 0 a i (t) T i (-1) -v ∞, L , σ 2 = k M, 2 n i = 0 a i (t) T i (1) n i = 0 ã i (t) T i (1) + Bi M, R n i = 0 a i (t) T i (1) -v ∞, R , σ 3 = k T, 1 n i = 0 a i (t) T i (-1) n i = 0 b i (t) T i (-1) + k T M, 1 n i = 0 a i (t) T i (-1) n i = 0 ã i (t) T i (-1) -Bi T, L n i = 0 b i (t) T i (-1) -u ∞, L -Bi T M, L n i = 0 a i (t) T i (-1) -v ∞, L , σ 4 = k T, 2 n i = 0 a i (t) T i (1) n i = 0 b i (t) T i (1) + k T M, 2 n i = 0 a i (t) T i (1) n i = 0 ã i (t) T i (1) + Bi T, R n i = 0 b i (t) T i (1) -u ∞, R + Bi T M, R n i = 0 a i (t) T i (1) -v ∞, R ,
b 1 =           0 0 . . . 0 ϑ 1 ϑ 2           , b 2 =           0 0 . . . 0 ϑ 3 ϑ 4           , b 3 =           0 0 . . . 0 σ 1 σ 2           and b 4 =           0 0 . . . 0 σ 3 σ 4           .
With all elements listed before, it is possible to set the system to be solved. The system of ODEs for solving the coupled heat and moisture transfer with two layers has the size of 4 • N with four additional algebraic expressions for the boundary and interface conditions. The initial condition is also given by Eq. (3.10) and the DAE system is solved by the Matlab ™ solver ODE15s. In this work, the approach was presented for a two-layered wall for the sake of clarity, knowing it can be easily extended to any number of layers and any other boundary conditions at their interfaces.

Numerical benchmark

To analyze the accuracy of the proposed method, the error between the solution Y num , obtained by the Spectral or IMEX methods, and the reference solutions Y ref , are computed as functions of x using the following Euclidean norms:

ε Y 2 ( x ) def := 1 N t N t j = 1 Y num j ( x , t j ) -Y ref j ( x ,
( x ) : ε Y ∞ def := sup x ∈ 0 , L ε Y 2 ( x ) .

Single-layered domain

Simulations of one-dimensional coupled heat and moisture transport are carried out with the spectral method, the IMplicit-EXplicit (IMEX) method and a reference solution to verify the applicability of the spectral method. The IMEX scheme approximates the continuous operator to order O (∆x 2 + ∆t) . The advantage of this semi-implicit scheme over the fully implicit one is to avoid sub-iterations in the solution procedure and, at the same time, being stable and consistent.

This case study considers moisture-dependent coefficients c M , k M , c T , k T and k M as illustrated in Figure 17 of Appendix B. Their variations are similar to the load bearing material from Hagentoft [START_REF] Hagentoft | Assessment Method of Numerical Prediction Models for Combined Heat, Air and Moisture Transfer in Building Components: Benchmarks for One-dimensional Cases[END_REF], which are presented in Table 4 of Appendix A. The material of 0.1 m of length has uniform vapour pressure and temperature initial conditions, in which P v, 0 = 1.16 • 10 3 Pa and T 0 = 293 K. The ambient temperature and relative The Spectral method is composed by N = 10 modes with m = 15 . The ODE15s was used to solve System (3.9), with a tolerance set to 10 -5 . This solver is adaptive in time, although it can provide the integration in time for the instants required. For this case, the Spectral method was compared to the IMEX and to a reference solution computed using the Chebfun Matlab ™ toolbox [15]. All solutions have been computed with ∆x ⋆ = 10 -2 . For the time discretization, the IMEX solution was computed with ∆t ⋆ = 10 -2 and the Spectral with ∆t ⋆ = 10 -1 , to have the same order of error of the solution.

Vapour pressure variations at the boundaries are shown in Figure 3(a). In the first four days, it is possible to observe the influence of the initial condition in which the vapour pressure rises significantly, meaning that when the simulation started the material was not in balance with the ambient environment. The vapour pressure at x = 0 m and x = 0.1 m oscillates according to the boundary conditions while also retaining the moisture. As the convective vapour transfer coefficient is higher at the left side, the material will exchange more with the external environment making variations more visible.

Figure 3(b) shows the temperature variations at the boundaries of the material. The temperature slowly oscillates according to the temperature ambient conditions and also with the vapour pressure variations. As one could expect, vapor pressure and temperature values inversely oscillate. Results of the error ε 2 as a function of x are shown in Figure 4(a). The errors associated to the IMEX scheme and to the Spectral method are of order of O (10 -5 ), for temperature, and of order of O (10 -4 ) for moisture. The computer run time used to perform this case and the error ε ∞ are given in Table 2, which shows the spectral method is more efficient, even compared to the IMEX, which is an improved method based on the Euler scheme as it does not need sub-iterations. The number of degrees of freedom (DOF) of the Spectral solution is considerably lower if compared to the IMEX solution, which impacts directly on the computational time, making the Spectral approach 7 times more efficient for this case.

Figure 4(b) displays the last coefficient a n of the Spectral-ROM solution for the temperature u and vapour pressure v solutions. The last spectral coefficient is the smallest one and it gives the order of approximation of the residual [START_REF] Gasparin | UMR 5127[END_REF]Page 51]. It works as an upper bound on the error. For example, the error of the u solution in Figure 4(a) is of order of O (10 -5 ), which has the same magnitude of the last coefficient, as can be seen in Figure 4(b).

To illustrate the rate of convergence of the methods, Figure 5(a) presents the error ε ∞ as a function of the number of modes N for the Spectral solution, and, as a function of the number of spatial nodes N x for the IMEX solution. As the number of modes increases, the solution of the Spectral method converges exponentially and stabilizes with approximatively 13 modes, reaching an accuracy of order of O (10 -5 ) , which is equivalent to the tolerance set on the solver ode15s. On the other hand, the IMEX method converges slower than the Spectral method. It needs at least 100 spatial nodes to reach the same accuracy of the Spectral method: the CPU time of each simulation was measured and it is presented in Figure 5(b). The computational time increases faster for the Spectral method. However, the spectral solution does not need many modes to converge to an acceptable accuracy (N ≃ 10). The Spectral method has demonstrated good agreement to represent the physical model and the fidelity of the model is not deteriorated with the use of this approach.

Multi-layered domain

This case study considers a porous wall formed by 2 layers: 8-cm load bearing material and 2-cm finishing material. Figure 6 shows its schematic representation. The first layer has a faster liquid transfer, while the second layer acts as a hygroscopic finish. The properties of these materials are given in Tables 4 and5. Considering the temperature range of interest in building applications, temperature dependence was neglected when compared to their dependence on moisture content, with the transport coefficients calculated as a function of the moisture content. Initial conditions are considered to be uniform over the spatial domain, with the initial temperature of T 0 = 293.15 K and the initial vapour pressure of P v, 0 = 1.16 • 10 3 Pa, regarding to a relative humidity of 50 % . The boundary conditions oscillate sinusoidally during 7 days of simulation, which are represented in Figures 2(a Simulations with the Spectral-ROM were performed using the ode15s, with a tolerance set to tol = 10 -5 , with N = 8 modes and m = 13 . The time is incremented with a discretization of ∆t ⋆ = 10 -1 and the Spectral solution is project to N x = 101 spatial nodes. The IMEX solution is computed for a discretization parameter of ∆t ⋆ = 10 -2 and ∆x ⋆ = 10 -2 .

Figure 8(a) presents three profiles of vapour pressure field, at the time instants t = {20, 30, 40} h , which are included on the interval of the first incoming raining flow. Variations of vapour pressure are more significant on the first layer because this material is more hygroscopic than the second one. At the interface, one can notice a discontinuity on the derivative. However, the continuity of the field and the one of the flow are assured. Additionally, Figure 8(b) presents the temperature profiles for the same time instants. Variations on the boundaries occur mainly due to the variations of P v, ∞ in both fields, and not only due to the rain flow, which has an opposite effect on the temperature.

The evolution of the temperature and the vapour pressure at the boundary surfaces (x = 0 m and x = 0.1 m) is shown in Figures 9(b) and 9(a), respectively. The vapour pressure varies according to the sinusoidal fluctuations of the boundary conditions until the rain hits the surface. In the first peak, it is not possible to observe high changes but in the second peak, as the material accumulates moisture, the vapour pressure is suddenly augmented and it diffuses through both layers. In fact, only when the vapour pressure is higher inside both material that is possible to observe the impact of the rain flux because of the material properties, more precisely due to the sorption isotherm. As can be observed in Figure 9(b), the temperature at the boundaries of the composite wall varies according to the gradients of vapour pressure rather than the temperature gradients.

By examining the profiles and evolution of temperature and vapour pressure profiles, it is possible to notice that the solutions of the Spectral and IMEX methods are both in a good agreement with the reference solution given by Chebfun. The distribution of the error ε 2 as a function of x are shown in Figure 10. The L ∞ error of the spectral solution is ε v ∞ = 2.9 • 10 -3 for the vapour pressure and ε u ∞ = 1.29 • 10 -4 for the temperature. These values depend on the order of the time discretization, on the chosen tolerance of the solver and on the number of modes. For the IMEX solution, the L ∞ error is ε v ∞ = 2.6 • 10 -3 for the vapour pressure and ε u ∞ = 1.1 • 10 -4 for the temperature. To obtain the same accuracy of the solutions, the Spectral-ROM was 7 times more efficient than the IMEX approach, as it can be seen in Table 3. One should recall that the IMEX has the same accuracy of the Euler explicit and is more efficient, in terms of computer run time, than the largely used Crank-Nicolson scheme. Just for curiosity, the computational time of the Chebfun simulations have been added. As one can notice they are much higher than the values obtained by the other methods. In fact, Chebfun is made for numerical computing of a wide range of problems. As it does not know about the specific problem, it does not take advantage of the problem structure. Anyway, it is a good tool for comparing purposes.

Figure 11(a) presents the last spectral coefficients of u and v solutions for material 1, while Figure 11 upper limit, as previous mentioned. In material 1, spectral coefficients needed more modes to have an acceptable solution, while in the other material, an accurate solution can be built with a lower number of modes. This difference occurs due to the properties of material 1 that are more nonlinear than those of material 2. One may notice in Figure 11(b) a peak after day 5, which disappears shortly after. This may occur due to aliasing errors. However, the magnitude of this error is very small if compared with the magnitude of the solution and it does not occur on the other spectral coefficients. As the problem demands only a few modes (around 10), this kind of error does not really affect the final solution.

In addition, Figure 12 displays the Fourier power spectrum function of the signal frequency per unit of time, on the left and right boundaries, for u and v solutions. Oscillations occurring for u , from 2 • 10 -1 Hz, are attributed to aliasing errors. However, the power of this frequency is very low if compared with the highest peak, corresponding to a difference of 6 orders of magnitude.

The results for the sensible, latent and total heat fluxes at the left boundary x = 0 m are given graphically in Figure 13(a), the boundary that receives the rain flow. The sensible and latent heat fluxes have high values but with opposite signs. Although, as they do not have the same value, they do not cancel each other, as shown for the total heat flux, which is the sum of the latent and sensible heat fluxes. Furthermore, Figure 13(b) presents the total moisture flow at the same boundary. The rain flow is observed in the two maximal

Validation of the model

In this section, the physical model and the spectral solution are compared with experimental data gathered from the French project HYGRO-BAT [START_REF] Hygro-Bat | Vers une méthode de conception HYGRO-thermique des BA-Timents performants[END_REF]. The measurements were performed at the French laboratory LOCIE (Laboratory of Optimisation of the Conception and Engineering of the Environment) [START_REF] Bouddour | Heat and mass transfer in wet porous media in presence of evaporation-condensation[END_REF]. One-dimensional coupled heat and moisture transfer through a single-layered wall is monitored by sensors placed inside of the material and on its surfaces. In this case, they have not considered the liquid transfer in the moisture balance equation -Eq. (2.1a). Surface sensors provide boundary conditions for the coupled simulation, while the other sensors provide reference measurements for the model validation.

The relative error ǫ is computed to compare simulations with the experimental data and it is defined as:

ǫ ( t ) def := Y num k ( x k , t ) -Y meas k ( x k , t ) 2 Y meas k ( x k , t ) ,
where Y num k is the computed solution and Y meas k is the measured data.

Experimental setup

The wall considered in this study is composed of a 16cm layer of wood fibre material, which is subjected to variations in terms of temperature and relative humidity for a 14-days period. The material properties are given in Table 6. The reference data for the evaluation is provided by temperature and humidity sensors SHT75 Sensirion, located at x = {4, 8, 12} cm within the wall.

The sensors have a measurement uncertainty of ± 0.3 • C for temperature and of ± 0.018 for relative humidity. Furthermore, the uncertainty regarding the position of the sensors is of ± 1 cm for the sensors located at x = {4, 12} cm and of ± 0.5 cm for the other sensor. The uncertainty on the positions are different because the sensor at x = {4, 12} cm have been settled by perforating a whole in the material layer.

The measured temperature and relative humidity at the interior and exterior boundaries are given in Figures 14(a) and 14(b). The gray color around the curves represent the uncertainties related to the measurements. At the interior boundary, the temperature is set to approximately 24 • C and the relative humidity set to 40% in the first week and to 70% in the second week. The exterior temperature and relative humidity values are given by their measurement at the boundary, which is filtered by a 2 cm of a coating layer, excluding solar radiation and driven rain phenomenon. Thus, both boundaries are expressed as Dirichlet-type conditions for the model validation.

Simulation

Simulations are performed for a single-layer wall of a material that separates two well defined environments. The initial condition for temperature and vapour pressure fields are given by an interpolation of the measurements at x = {0, 4, 8, 12, 16} cm at the first time instant.

The Spectral method is composed of N = 8 modes with m = 13 . The ODE15s was used to solve numerically System (3.9), with a tolerance set to 10 -5 . The solution has been computed with a time step of ∆t ⋆ = 0.1 , the equivalent of 6 min.

Simulations are compared with the dynamic temperature and vapour pressure. Figure 15 presents the measured data and the simulation results in each one location of the sensors, at x = {4, 8, 12} cm.

The error between the predicted solution and the experimental observations is between the uncertainties of the sensors during almost all simulation period. For the temperature evolution, the highest discrepancy is when the moisture at the left boundary changes from 40% to 70% of relative humidity. From the 7 th day, the measured temperature rises more than the calculated one as shown in Figure 15(a), which are higher than the uncertainties of the sensors. The difference between the simulated and measured temperature reaches a maximum of 2 • C in this period. The temperature simulated in at x = {8, 12} cm are presented in Figures 15(c) and 15(e). The predicted values are closer to the measurements than in x = 4 cm, and they follow considerably well variations of the outside boundary. It seems that the total diffusion coefficient of the temperature used for the simulation is higher than the real one.

The influence of the step on the relative humidity can strongly be observed at x = 4 cm in Figure 15(b) and less far from this boundary. Simulations are able to follow the variations of the measured vapour pressure better than the temperature. In Figures 15(d) and 15(f), the difference between simulations and measurements become more important, with the incoming moisture flow, reaching a maximum difference of 70 Pa and 80 Pa, respectively. The absolute difference is higher for x = 4 cm as it is closer to the left boundary and consequently to the incoming flow. The discrepancies, come from the physical model, which does not consider liquid transport neither hysteresis. Another explanation is the fact that the properties also were estimated with and admissible error, which influences the predicted solution. This difference on the vapour pressure simulations impacts on the predictions of temperature.

The relative error ǫ of the computed solutions are present in Figure 16. Solutions simulated with the Spectral method showed a good agreement with the reference data. The maximum relative error for the temperature solution is of 0.7% and for the vapour pressure solution is of 5.6%. The average error on the dynamic profiles is 0.2% for temperature and 1.75% for vapour pressure, which are close to the values obtained by Rouchier et al. [START_REF] Bouddour | Heat and mass transfer in wet porous media in presence of evaporation-condensation[END_REF].

The spectral method can provide a very accurate solution of the physical model. In this case, the last spectral coefficient for the temperature is of order of b n ≃ O (10 -5 ) and for the vapour pressure of order of a n ≃ O (10 -3 ) , for N = 8 . The difference observed in Figure 15(a) can be caused by differences in the mathematical model and by the uncertainties of the estimated properties that are not taken in consideration [START_REF] Busser | Comparison of model numerical predictions of heat and moisture transfer in porous media with experimental observations at material and wall scales: An analysis of recent trends[END_REF]. Despite of this disparity, simulated results can be considered satisfactory to predict the heat and moisture transport.

As the boundary conditions are of the Dirichlet-type, the IMEX approach is more efficient in this case, taking only 4 times more CPU time than the Spectral reduced-order model approach. For a single simulation, it may not appear important, however, when the direct model has to be simulated thousands of times as in [START_REF] Bouddour | Heat and mass transfer in wet porous media in presence of evaporation-condensation[END_REF], the gains become much more expressive.

Conclusions

In the present work, we showed that the unsteady heat and moisture transfers problem is solved efficiently with the spectral reduced-order model. This approach has been evaluated on two numerical unidimensional case studies of heat and moisture transfer in porous media. Each case aimed at exciting the nonlinear properties of the material to induce sharp profiles of temperature and vapour pressure. The first case considered a single material layer with sinusoidal boundary conditions. The second case took into account the rain effect at one boundary through a multi-layered material. For all cases, the spectral method has shown a high accuracy and perfect agreement with the respective reference solutions. The maximum global error was of the order of O (10 -3 ) . The advantage of the proposed method is the low computational burden 7 times more efficient than the IMEX approach. In sensitivity analysis and inverse problems, when one has to simulate the system of differential equation many times, the spectral approach becomes very attractive.

Besides the numerical benchmark, one case with experimental data was performed to confirm the physical model and the spectral-based solution method. Results of the simulation showed a good agreement with the experimental data, which confirms all the process of simulation. Results can be improved by adding the hysteresis and the liquid transport effects. Although, this information was not available for this material and this experiment.

Further research should be dedicated to multiple space dimensions and to the application considering both the diffusive and advective transfer, which complexity increases considerably due to the nonlinearities of the problem. For the dimensionless properties of the material, they can be written as: C.2. Case from Section 4.2

Nomenclature

c ⋆ M (v) =
The dimensionless temperature and vapour pressure at the boundaries are written as in the previous case and also the materials properties of material 1. Thus, for the second material, properties are written as:

  the total moisture transfer coefficient under vapour pressure gradient,k T M def := L v δ v :the heat coefficient due to a vapour pressure gradient, k T def := λ : the heat transfer coefficient under temperature gradient, 0 c 0 + w c w :

Figure 1 .

 1 Figure 1. Schematic representation of the domain division with the real domain (a) transformed linearly to obtain the spectral domain (b).

which are included in b 3 and b 4

 4 and set equal to zero. Vectors b 1 , b 2 , b 3 and b 4 are N × 1 column vectors with the form:

Figure 2 .

 2 Figure 2. Ambient temperature (a) and relative humidity (b).

Figure 3 .

 3 Figure 3. Vapour pressure evolution at the boundaries (a) and vapour pressure profiles inside the material (b).

Figure 4 .

 4 Figure 4. Error of both solutions (a) and last spectral coefficient (b).

Figure 5 .

 5 Figure 5. Error ε ∞ (a) and the equivalent CPU time (b) as a function of the number of modes N for the Spectral solution and as a function of the number of spatial nodes N x for the IMEX solution.

Figure 6 .Figure 7 .

 67 Figure 6. Schematic representation of the two-layers wall.

  ) and 2(b). The convective mass and heat transfer coefficients are set toh M = 2 • 10 -7 s/m, h M, R = 3 • 10 -8 s/m, h T, L = 25 W/(m 2 • K) and h T, L = 8 W/(m 2 • K).The liquid water flow (rain) has two peaks, one at 42 h and the other one at 126 h, reaching a maximum value of 1.53 • 10 -4 kg/(m 2 • s) as shown in Figure 7(a), which generates a sensible heat flux of 11 W/m 2 and 14 W/m 2 as displayed in Figure 7(b).

Figure 8 .

 8 Figure 8. Vapour pressure and (a) temperature (b) profiles at t = {20, 30, 40} h.

Figure 9 .

 9 Figure 9. Evolution of vapour pressure (a) and temperature (b) at the boundaries of the composite wall.

Figure 10 .

 10 Figure 10. Error ε 2 given in function of x .

Figure 11 .

 11 Figure 11. Spectral coefficients of material 1 (a) and of material 2 (b).

Figure 12 .

 12 Figure 12. Fourier power spectrum of the Spectral solution computed in the right and left boundaries of the domain.

Figure 13 .

 13 Figure 13. Latent, sensible and total heat fluxes (a) and moisture flow (b) at the left boundary (x = 0) .

Figure 14 .

 14 Figure 14. Ambient temperature (a) and relative humidity (b) of the experimental investigation.

Figure 15 .

 15 Figure 15. Measured and predicted temperature and vapour pressure values at x = {4, 8, 12} cm.

Figure 16 .

 16 Figure 16. Relative errors regarding temperature (a) and vapour pressure (b) values at different positions within the porous material.

- 4 .

 4 15 • 10 6 • v 3 + 7.76 • 10 6 • v 2 -5.44 • 10 5 • v + 3.53 • 10 6 v 4 + 9.18 • 10 5 • v 3 -2.47 • 10 6 • v 2 + 1.68 • 10 6 • v + 5493 , k ⋆ M (v) = -4.23 • v 8 -4.901 • v 7 + 120.6 • v 6 -340.7 • v 5 + 417.8 • v 4 -255.9 • v 3 + 77.15 • v 2 -10.33 • v + 1.57 , c ⋆ T (v) = 9.327 • exp(-2.4 • 10 -4 • v) + 1.457 • 10 -14 • exp(15.58 • v) , k ⋆ T (v) = 0.9996 • exp(-8.813 • 10 -4 • v) + 5.65 • 10 -15 • exp(15.58 • v) , k ⋆ T M (v) = 0.1276 • exp(-1.651 • 10 -4 • v) ,withk M, ref = 5.4712 • 10 -9 s and k T, ref = 0.5021 W/(m • K) .

2

 2 (v) = 1.221 • v -0.878 , k ⋆ M, 2 (v) = -1.084 • 10 -4 • v 15.44 + 11.34 , c ⋆ T, 2 (v) = 4.52 • exp(0.1058 • v) + 1.79 • 10 -11 • exp(12.81• v) , k ⋆ T, 2 (v) = 0.8686 • exp(0.1414 • v) + 9.498 • 10 -7 • exp(7.968 • v) , k ⋆ T M, 2 (v) = -1.884 • 10 -11 • exp(14.25 • v) + 1.216 • exp(-0.0284 • v) ,with k M, ref and k T, ref equal to the previous case. In addition, the rain flow is expressed as:g ⋆ ∞ (t ⋆ ) = 2.4 • sin(π t ⋆ /84) 70 .C.3. Case from Section 5.2The boundary conditions are gathered from the experimental data and just admensionalized. The dimensionless form of the initial condition are:u 0 (x ⋆ ) = -0.08806 • (x ⋆ ) 4 + 0.1688 • (x ⋆ ) 3 -0.1143 • (x ⋆ ) 2 -0.01621 • x ⋆ + 1.015 , v 0 (x ⋆ ) = -0.408 • (x ⋆ ) 4 + 1.188 • (x ⋆ ) 3 -1.053 • (x ⋆ ) 2 + 0.08969 • x ⋆ + 1.092 .The dimensionless properties of the material can be written as:c ⋆ M (v) = -0.663 • v + 37.52 , k ⋆ M (v) = 0.007289 • v + 0.9854 , c ⋆ T (v) = 0.08587 • v + 16.53 , k ⋆ T (v) = 0.0005546 • v + 0.9989 , k ⋆ T M (v) =3.465 • 10 -5 • v + 0.004684 , with k M, ref = 3.34 • 10 -11 s and k T, ref = 6.98 • 10 -2 W/(m • K) .

Figure 17 .

 17 Figure 17. Coefficients c M (a), c T (b), k M (c), k T (d) and k T M (e) of the load bearing material, of the finishing material and of the wood fibre.

  )where P v, ∞ [Pa] stands for the vapour pressure far from the surface and, h M [s/m], is the convective moisture transfer coefficient. If the bounding surface is in contact with the outside air, g ∞ [kg/(m 2 • s)] is the liquid flow from wind driven rain. The normal n assumes + 1 or -1 at the left or right boundary sides.

Table 1 .

 1 Hygrothermal properties of water.

	Property	Value		Unit
	Heat capacity, c w Latent heat of evaporation, L v 2.5 • 10 6 4180 Water gas constant, R v 462 Density, ρ l 1000		[J/(kg • K)] [J/kg] [J/(kg • K)] [kg/m 3 ]
	Saturation pressure, P s (T )	997.3 •	120.6 T -159.5	8.275

Table 2 .

 2 Some features of the one-layer heat and moisture transfer case.

		IMEX	Spectral N = 10
	∆t ⋆ ε u ∞ ε v ∞ DOF	1.00 • 10 -2 1.48 • 10 -5 1.43 • 10 -4 200	1.00 • 10 -1 3.30 • 10 -5 2.31 • 10 -4 20
	CPU time (s)	6.46	0.94
	CPU time (%)	100	15

Table 3 .

 3 Computer run time required for the numerical schemes to perform simulations.

	Case study Spectral IMEX Chebfun Ratio (IMEX/Spectral)
	1 layer	0.9 s	6.4 s	46.6 s	7.1
	2 layers	1.8 s	12.5 s	227.9 s	6.9

Acknowledgements

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -in the framework of the International Cooperation Program CAPES/COFECUB (Grant #774/13). The Authors acknowledge also the support from CNRS/INSIS and Cellule Énergie under the grant MN4BAT-2017. Finally, Professor Mendes thanks the Laboratory LAMA UMR 5127 for the warm hospitality during his visits in 2018, which were supported by the project MN4BAT under the AAP Recherche 2018 programme of the University Savoie Mont Blanc.

A. Material properties Table 4. Hygrothermal properties of the load bearing material [START_REF] Hagentoft | Assessment Method of Numerical Prediction Models for Combined Heat, Air and Moisture Transfer in Building Components: Benchmarks for One-dimensional Cases[END_REF].

Property

Value Unit Volumetric heat capacity

Sorption isotherm w (φ) = 47.1 1 + -1692.94 • ln(φ)

Table 5. Hygrothermal properties of the finishing material [START_REF] Hagentoft | Assessment Method of Numerical Prediction Models for Combined Heat, Air and Moisture Transfer in Building Components: Benchmarks for One-dimensional Cases[END_REF].

Table 6. Hygrothermal properties of the wood fibre [START_REF] Bouddour | Heat and mass transfer in wet porous media in presence of evaporation-condensation[END_REF].

B. Coefficients of the different materials

See Figure 17.

C. Dimensionless values

Some reference values were used on all the simulation cases. The reference time is t ref = 1 h , the equivalent to 3600 s. The reference temperature was T ref = 293.15 K and the reference of the vapour presure was P v, ref = 1166.9 Pa . The reference length is the total length of the spatial domain L ref = L m , so it is possible to have a dimensionless domain between x ⋆ = [ 0, 1 ] .

C.1. Case from Section 4.1

The temperature boundary conditions are expressed as:

and, of the vapour pressure boundary conditions as: