Julien Berger

Denys Dutykh

Nathan Mendes

Bolatbek Rysbaiuly
email: b.rysbaiuly@mail.ru

A new model for simulating heat, air and moisture transport in porous building materials

Keywords: heat and mass transfer, porous material, Benchmarking with experimental data, Scharfetter-Gummel numerical scheme, Du Fort-Frankel numerical scheme MSC: [2010] 35R30 (primary), 35K05, 80A20, 65M32 (secondary) PACS: [2010] 44.05.+e (primary), 44.10.+i, 02.60.Cb, 02.70.Bf (secondary) mass transfer, porous material, Benchmarking with experimental data, Scharfetter-Gummel numerical scheme, Du Fort-Frankel numerical scheme

This work presents a detailed mathematical model combined with an innovative efficient numerical model to predict heat, air and moisture transfer through porous building materials. The model considers the transient effects of air transport and its impact on the heat and moisture transfer. The achievement of the mathematical model is detailed in the continuity of Luikov's work. A system composed of two advection-diffusion differential equations plus one exclusively diffusion equation is derived. The main issue to take into account the transient air transfer arises in the very small characteristic time of the transfer, implying very fine discretisation. To circumvent these difficulties, the numerical model is based on the Du Fort-Frankel explicit and unconditionally stable scheme for the exclusively diffusion equation. It is combined with a two-step Runge-Kutta scheme in time with the Scharfetter-Gummel numerical scheme in space for the coupled advection-diffusion equations. At the end, the numerical model enables to relax the stability condition, and, therefore, to save important computational efforts. A validation case is considered to evaluate the efficiency of the model for a nonlinear problem. Results highlight a very accurate solution computed about 16 times faster than standard approaches. After this numerical validation, the reliability of the mathematical model is evaluated by comparing the numerical predictions to experimental observations. The latter is measured within a multi-layered wall submitted to a sudden increase of vapor pressure on the inner side and driven climate boundary conditions on the outer side. A very satisfactory agreement is noted between the numerical predictions and experimental observations indicating an overall good reliability of the proposed model.

Introduction

Moisture is a key factor on durability and performance of buildings. An excessive level compromises the construction quality, impacts the indoor air quality and the thermal comfort, as well as the building energy efficiency [START_REF] Berger | Factors governing the development of moisture disorders for integration into building performance simulation[END_REF]. As a consequence, a number of models for predicting the impact of moisture on building energy efficiency are proposed in the literature. A primary overview may be consulted in [START_REF] Mendes | Numerical methods for diffusion phenomena in building physics[END_REF]. Among the physical phenomena, the air transfer through the porous building media has a crucial impact on the amount of moisture. Diverse studies enhance these effects using both experimental and numerical results [START_REF] Desta | Experimental data set for validation of heat, air and moisture transport models of building envelopes[END_REF][START_REF] Kalamees | Moisture Convection Performance of External Walls and Roofs[END_REF][START_REF] Wang | Effect of air leakage on the hygrothermal performance of highly insulated wood frame walls: Comparison of air leakage modelling methods[END_REF].

Several numerical models are proposed in the literature for the prediction of physical phenomena of coupled heat, air and moisture transport through porous building materials. Their physical representations are based on the mass conservation laws for the dry air, vapor and liquid water, as well as the energy conservation law as detailed in the early work of Luikov [31]. As the continuity of his work, the numerical models proposed in the literature can be divided into two main groups. The first group considers the three evolution differential equations to compute the temperature, the mass content and the air pressure in the porous media. In [11], a model is proposed for the simulation of transfer through hollow porous blocks. It is based on an implicit finite-difference numerical scheme. More recently, in [START_REF] Abahri | Sensitivity analyses of convective and diffusive driving potentials on combined heat air and mass transfer in hygroscopic materials[END_REF][START_REF] Mnasri | Numerical analysis of heat, air, and moisture transfers in a wooden building material[END_REF], the commercial COMSOL TM software is used to propose a numerical model for such physical problems. As mentioned by the authors, the scheme is based on an explicit in time finite element approach. The main drawback of these numerical models is their computational cost. The implicit approach requires costly sub-iterations at each time step to handle severe nonlinearities of the problem. The explicit scheme requires very fine time steps to satisfy the so-called Courant-Friedrichs-Lewy (CFL) stability conditions. Indeed, the characteristic time of air transfer is very small compared to the ones for heat and mass transfer.

To handle this computational issue, the second group of models does not consider the transient phenomena for the air transport through the porous matrix. In other words, the evolution differential equation for air transfer is transformed into a simple steady boundary value problem which is solved at prescribed time instances. It enables somehow to relax the stability conditions and, therefore, to save computational efforts comparing to the models from the first group. Some examples can be found in [START_REF] Tariku | Transient model for coupled heat, air and moisture transfer through multilayered porous media[END_REF] or [START_REF] Langmans | A quasi-steady state implementation of air convection in a transient heat and moisture building component model[END_REF] for one-dimensional transfer. The former references uses an explicit scheme provided by the commercial COMSOL TM software. The latter is based on an implicit scheme based on the generic ODE solver from the SUN-DIALS solver package [23]. More recently, a numerical model is proposed in [START_REF] Belleudy | Experimental and numerical investigations of the effects of air leakage on temperature and moisture fields in porous insulation[END_REF][START_REF] Belleudy | A 2D model for coupled heat, air, and moisture transfer through porous media in contact with air channels[END_REF] for the simulation of two-dimensional transfer in building structures. It is also based on an explicit scheme from COMSOL TM software. The assumption of neglecting the transient term is justified by the low velocities occurring through the porous matrix. As mentioned by authors, the numerical predictions are reliable only in the context of simulations with standard hourly climate driven boundary conditions. It should be noted that this condition is not often satisfied, particularly for the air pressure surrounding building walls as detailed in [27]. Moreover, even if the restriction on the time discretisation is relaxed, the numerical models are based on standard approaches and still have a high computational cost.

Therefore, the goal is to propose an efficient numerical model considering the three transient equations to improve the reliability of its predictions. It requires to be accurate with a reduced computational cost. To address this issue, this paper proposes to use the Scharfetter-Gummel scheme combined with a two-step Runge-Kutta approach. The Scharfetter-Gummel numerical scheme was proposed in 1969 in [START_REF] Scharfetter | Large-signal analysis of a silicon Read diode oscillator[END_REF] with very recent theoretical results in [START_REF] Gosse | Viscous equations treated with L-splines and Steklov-Poincare operator in two dimensions[END_REF][START_REF] Gosse | L-Splines and Viscosity Limits for Well-Balanced Schemes Acting on Linear Parabolic Equations[END_REF]. In the context of building porous media, it is successfully applied in [5] to water transport and then in [START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF] to combined heat and moisture transfer. The contributions of the present paper is two fold. First, the model proposed in [5] is extended by including the air transport equation. Then, the extension of the Scharfetter-Gummel approach to a system of three coupled advection-diffusion equations is proposed. The combination with a two-step Runge-Kutta scheme is investigated in order to relax the stability restrictions on the choice of the time discretisation.

The paper is organized as follows. Section 2 presents the demonstration of the mathematical model to describe the physical phenomena and its dimensionless formulation. Then, Section 3 presents the numerical method to solve the system of three differential advection-diffusion equations. In Section 4, one case studie is considered to validate the numerical model. The purpose is to quantify the accuracy and efficiency in terms of computational time and relaxation of the stability condition. For each case, a reference solution is proposed, computed by a numerical pseudo-spectral approach. In the last Section, the reliability of the numerical predictions is evaluated by comparing them to experimental observations. A wall composed of two layers of wood fiberboard is submitted to a controlled environment on the inner side and to climate driven boundary condition on the outer side. Three points of measurements of temperature and vapor pressure within the wall are used for comparison purposes.

Formulation of the physical phenomena

First, the mathematical model including the governing equations to describe the physical phenomena is presented.

Porous medium basics

The term moisture is used to designate the water vapor, denoted by index 1, and the liquid water, denoted by index 2 . The gas mixture (moist air) is composed of water vapor and dry air indexed by 3 . The porous matrix of the material is symbolized by index 0 . The volumetric concentration w i kg/m 3 of specie(s) i , which varies in space and time, is defined by

w i def := lim V → 0 m i ( V V .
We denot by m i ( V kg the mass of the specie(s) in the representative volume V m 3 . The representative volume of material V is composed of the volume of the solid porous matrix V 0 , the volume of water in its liquid form V 2 and the volume of the voids formed by the pores V p . Another important property of a porous material is the so-called sorption isotherm curve, which provides a relation between the moisture content w 12 , in liquid and vapor phases, with the relative humidity:

w 12 def := m 12 V = f φ , (2.1)
where f can be a fitted function from experimental data. By introducing this property, an instantaneous thermodynamic equilibrium is assumed in the material. As reported in [START_REF] Soudani | Assessment of the validity of some common assumptions in hygrothermal modeling of earth based materials[END_REF], most of the models assume the water content provided by the sorption curve corresponds approximately to the measure of the liquid water content, since the mass of the vapor phase is negligible compared to the one of liquid, i.e., m 1 ≪ m 2 . In the present work, this simplification is not considered. The measure of the sorption curve includes both vapor and liquid water masses. The saturation rate σcorresponds to the volume of liquid water out of the volume of pores. The estimation of the saturation is a difficult task since the weighting of a material only provides the amount of both liquid and vapor within the porous structure. When the assumption m 12 ≈ m 2 is considered, the saturation rate is evaluated directly by the sorption curve. Presently, the saturation rate is given by:

σ = 1 ρ 2 R 1 T -P 1 R 1 T Π w 12 -P 1 .
Assuming the perfect gas law, the volumetric vapor and dry air content are given by:

w 1 = P 1 R 1 T Π 1 -σ , w 3 = P 3 R 3 T Π 1 -σ .
Last, the liquid volumetric mass content is given by:

w 2 = w 12 -w 1 = w 12 - P 1 R 1 T Π 1 -σ . (2.2)
Interested readers may consult [31] for further information on these definitions. Using Whitaker volume averaging method to link the microscopic and macroscopic approaches [44, 45], the conservation law for each specie can be written according to the following differential equation:

∂w i ∂t = -∇ • j c , i + I i , i ∈ 1 , 2 , 3 , (2.3)
where I i kg/(m 3 . s) is the volumetric term source or sink of specie i . Since chemical reactions have not been taking into account, I 3 = 0 . Moreover, it is assumed that the temperature is always higher than the freezing point so no ice water may appear in the material. Thus, the following relation is stated:

I 1 + I 2 = 0 .
The flow of mass is denoted by j c , i kg/(m 2 . s) . Due to the airflow occurring through the porous medium, the mechanism of mass transfer is driven by both diffusion j d and advection j a flows:

j c = j d + j a .
It is important to mention that the gravity effects are neglected since the transport phenomena will be investigated on a horizontal plane perpendicular to the gravity forces.

Moisture mass balance

By summing Eq. (2.3) (i ← 1) and Eq. (2.3) (i ← 2), the moisture transfer equation is obtained:

∂ ∂t w 1 + w 2 = -∇ • j c , 1 + j c , 2 .
The term w 1 + w 2 is equal to the moisture content w 12 . Using the expression of the sorption isotherm curve Eq. (2.1), the partial differential of the moisture content becomes:

∂w 12 ∂t = ∂w 12 ∂φ 1 P s ∂P 1 ∂t ,
where P s is the saturation pressure, depending on temperature, and ∂ ∂φ w 1 + w 2 is the derivative of the sorption curve of the material. In this way, the moisture transfer equation becomes:

∂w 12 ∂φ 1 P s ∂P 1 ∂t = -∇ • j c , 1 + j c , 2 .

Moist air mass balance

Since the air content source term vanishes and adding Eq. (2.3) (i ← 1) to Eq. (2.3) (i ← 3), one obtains:

∂ ∂t w 1 + w 3 = -∇ • j c , 1 + j c , 3 + I 1 , (2.4)
where w 1 + w 3 accounts for the amount of mass of dry air and vapor in the whole mixture and can be expressed through the total pressure P :

w 1 + w 3 = P R 13 T Π ( 1 -σ ) , (2.5)
where P is the total pressure of the mixture of vapor and dry air. The term R 13 is the gas constant for the mixture of dry air and vapor. It is computed using the following expression:

R 13 = R 3 1 -P 1 P 3 1 -R 3 R 1 . (2.6)
Assuming the porosity to be invariant, the following equation is obtained by computing the differential of Eq. (2.5):

∂ w 1 + w 3 = Π ( 1 -σ ) R 13 T ∂P + Π P R 13 T ∂( 1 -σ ) - Π ( 1 -σ ) P R 13 T 2 ∂T . (2.7)
The variation of the gas mixture mass (composed of vapor and dry air) depends on the variation of the total pressure P , of the volume gas 1 -σ and of the temperature T . Thus, Eq.(2.4) becomes:

Π ( 1 -σ ) R 13 T ∂P ∂t = -∇ • j c , 1 + j c , 3 + I 1 + Π P R 13 T ∂σ ∂t + Π ( 1 -σ ) P R 13 T 2 ∂T ∂t .
The variation of the pressure of the mixture depends on the divergence of the fluxes, on the vapor source due to vaporization/condensation and on the variation of the gas volume and the temperature.

Energy balance

The energy balance equation is derived from the first law of thermodynamics, which states the internal energy variation in time is due to the balance of energy across the volume control, i.e., that the volumetric concentration of the enthalpy h J/kg equals the divergence of the conduction and enthalpy flux:

∂ ∂t h 0 ρ 0 + 3 i=1 h i w i = -∇ • j q + 3 i=1 h i j c , i , (2.8)
where ρ 0 kg/m 3 is the dry material density and c i J/(kg.K) the specific heat of each species. The expression of the volumetric water content w i is detailed in Section 2.1. By assuming a constant volume, Eq (2.8) becomes:

c 0 ρ 0 + 3 i=1 c i w i ∂T ∂t + 3 i=1 h i ∂w i ∂t = -∇ • j q - 3 i=1 h i ∇ • j c , i - 3 i=1 ∇ h i • j c , i .
(2.9)

Since I 3 = 0, by summing Eq. (2.3) (i ← 1), multiplied by h 1 , Eq. (2.3) (i ← 2), multiplied by h 2 and Eq. (2.3) (i ← 3), multiplied by h 3 , we obtain that:

3 i=1 h i ∂w i ∂t = 3 i=1 h i ∇ • j c , i + h 1 I 1 + h 2 I 2 .
Thus, Eq. (2.9) becomes:

c 0 ρ 0 + 3 i=1 c i w i ∂T ∂t = -∇ • j q -h 1 I 1 -h 2 I 2 - 3 i=1 ∇ c i T • j c , i .
We denote by r 12 def := h 1 -h 2 J/kg as the latent heat of vaporization. According to Kelvin's law, it is defined as [START_REF] Soudani | Assessment of the validity of some common assumptions in hygrothermal modeling of earth based materials[END_REF]:

r 12 = r • 12 + c 1 -c 2 • T -T ref -R 1 T ln φ .
(2.10)

Finally, the energy conservation equation yields:

c 0 ρ 0 + 3 i=1 c i w i ∂T ∂t = -∇ • j q -r 12 I 1 - 3 i=1 ∇ c i T • j c , i .
The temperature variations in the material depend on the variation of the thermal flux, on the energy due to vaporization and the enthalpy transport of each species by advection.

Volumetric vapor source

In order to write the system of governing differential equations, it is important to express the volumetric vapor source I 1 . For that, several possibilities are suggested in literature. In its original work [31], Luikov neglected the time variation of the vapor mass and expressed the source term as

I 1 = ∇ • j c , 1 .
(2.11)

Similar approach has been adopted in [START_REF] Kuenzel | Simultaneous heat and moisture transport in building components: one-and two-dimensional calculation using simple parameters[END_REF][START_REF] Tariku | Transient model for coupled heat, air and moisture transfer through multilayered porous media[END_REF]. This assumption was discussed in [START_REF] Soudani | Assessment of the validity of some common assumptions in hygrothermal modeling of earth based materials[END_REF], for the case of rammed earth material, by comparing the results using expression (2.11) and the following one:

I 1 = - ∂w 2 ∂t + ∇ • j c , 2 .
The latter expression assumes that the moisture concentration w 2 is obtained by the sorption isotherm w 12 ≈ w 2 . This assumption is not considered in this work, a third expression of the source is proposed:

I 1 = ∂w 1 ∂t + ∇ • j c , 1 ,
where the vapor concentration w 1 is expressed as:

w 1 = Π 1 -σ R 1 T .
The differential of w 1 gives:

∂ w 1 = Π ( 1 -σ ) R 1 T ∂P 1 + Π P 1 R 1 T ∂( 1 -σ ) - Π ( 1 -σ ) P 1 R 1 T 2 ∂T .
(2.12)

Then, we obtain the final expression of the volumetric source term I 1 :

I 1 = Π ( 1 -σ ) R 1 T ∂P 1 ∂t - Π P 1 R 1 T ∂σ ∂t - Π ( 1 -σ ) P 1 R 1 T 2 ∂T ∂t + ∇ • j c , 1 ,
which depends on the vapor pressure, the amount of liquid in the pores, the temperature and the variations of the vapor fluxes.

Expression of the fluxes

A filtration of the dry air occurs through the porous material. The vapor velocity and the air velocity are assumed to be the same, v 3 = v 2 = v , and can be expressed by the Darcy law, remembering that the gravity effects are neglected:

v = - k 13 µ ∇P , (2.13)
where k 13 m 2 is the intrinsic permeability of the material and µ is the dynamic viscosity of the fluid. The diffusion of vapor and air is small compared to the filtration transfer. Thus, the total flux of the vapor and dry air is written as:

j c , 1 + j c , 3 = w 1 + w 3 Π 1 -σ v ,
which can be rewritten using Eq. (2.13) for the velocity:

j c , 1 + j c , 3 = - w 1 + w 3 Π 1 -σ k 13 µ ∇P .
The vapor flow is driven by diffusion and advection:

j c , 1 = -k 1 ∇P 1 + w 1 Π 1 -σ v = -k 1 ∇P 1 + P 1 R 1 T v ,
where k 1 is the vapor permeability of the material and R 1 J/(kg.K) the vapor gas constant. The air filtration is sufficiently slow to have any influence on the liquid water. Thus, the total moisture flow is written as:

j c , 1 + j c , 2 = -k 2 ∇P 2 -k 1 ∇P 1 + P 1 R 1 T v .
It has been shown that the diffusion of moisture can be written using the vapor pressure P 1 , introducing the global moisture permeability k m [START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF]. Thus, the total moisture flux yields:

j c , 1 + j c , 2 = -k m ∇P 1 + P 1 R 1 T v .
The heat flux is expressed as:

j q = -k q ∇T + w 1 c 1 + w 3 c 3 T Π 1 -σ v .
It is driven by conduction and advection of the moist air.

Governing Equations

For the sake of clarity, we introduce the following advection coefficients:

a q def := w 1 c 1 + w 3 c 3 Π 1 -σ v , a v def := v R 1 T ,
along with the storage coefficients:

c a def := Π ( 1 -σ ) R 13 T , c as def := Π T P R 13 - P 1 R 1 , c at def := Π ( 1 -σ ) T 2 P R 13 - P 1 R 1 , c av def := Π ( 1 -σ ) R 1 T , c qv def := Π ( 1 -σ ) R 1 T , c qs def := Π P 1 R 1 T , c m def := ∂w 12 ∂φ 1 P s , c q def := c 0 ρ 0 + 3 i=1 c i w i - Π P 1 ( 1 -σ ) R 1 T 2 r 12 ,
and the permeability coefficients:

k a def := w 1 + w 3 Π 1 -σ k 13 µ , k v def := k 1 .
Finally, the set of partial differential governing equations for the physical problem of heat, air and moisture transfer through porous building material becomes:

c m ∂P 1 ∂t = ∇ • k m ∇P 1 -a v P 1 ,
(2.14a)

c q ∂T ∂t = ∇ • k q ∇T -a q T + r 12 ∇ • k v ∇P 1 -a v P 1 (2.14b) - 3 i=1 ∇ c i T • j c , i -r 12 c qv ∂P 1 ∂t + r 12 c qs ∂σ ∂t , c a ∂P ∂t = ∇ • k a ∇P -∇ • k v ∇P 1 -a v P 1 (2.14c) + c av ∂P 1 ∂t + c at ∂T ∂t + c as ∂σ ∂t .

Boundary, interface and initial conditions

For the set of governing equations (2.14a), (2.14b) and (2.14c), boundary, interface and initial conditions are provided for a one-dimensional problem. At the interface between the material and the ambient air, the liquid flow is imposed:

j c , 2 = g m ,∞ ,
where g m ,∞ kg/(m 2 .s) is the liquid water flux due to wind driven rain that crosses the porous surface. For the vapor phase, the total vapor flow is proportional to the surface and ambient conditions:

j c , 1 = -α m P 1 -P 1 , ∞ • n ,
where α m s/m is the surface vapor transfer coefficient and n is the surface normal unitary vector. Therefore, for the moisture transfer equation, the following boundary conditions are stated:

j c , 1 + j c , 2 = -α m P 1 -P 1 , ∞ • n + g m , ∞ .
Similarly, for the heat transfer, the flow is proportional to the surface and ambient conditions. An additional heat flux g q , ∞ is introduced to take into account the radiative and convective heat transfer:

j c , q = -α q T -T ∞ • n + g q , ∞ .
Therefore, the boundary condition for the heat transfer equation can be written as:

k q ∇T -a q T + r 12 ∇ • k v ∇P 1 -a v P 1 = -α q T -T ∞ • n + g q , ∞ -r 12 α m P 1 -P 1 , ∞ • n + r 12 g m , ∞ .
Finally, the air pressure at the interface between the material and the ambient air is imposed:

P = P ∞ .
At the interface x i between two materials A and B , a perfect contact is assumed so that air, moisture and heat fluxes are assumed to be continuous across the interfaces [START_REF] De Freitas | Moisture migration in building walls -Analysis of the interface phenomena[END_REF][START_REF] Guimaraes | Interface influence on moisture transport in buildings[END_REF]:

j c , A ( x i ) = j c , B ( x i ) .
As initial conditions, the fields may depend on the space coordinate x :

T = T i ( x ) , P = P i ( x ) , P 1 = P 1 , i ( x ) .
It is important to note that at t = 0 , the boundary and initial conditions must be consistent, i.e. the initial condition must verify the boundary conditions. The numerical model to solve the mathematical problem is presented in Section 3 based on a dimensionless formulation.

Elaborating an efficient numerical model

Since the physical phenomena have been described, the second part in the elaboration of a numerical model consists in detailing the numerical method to solve the physical problem. For this, it is of major importance to define the strategy of building a numerical model that reduces the computational effort and maximize the accuracy of the solution. First, since we have a nonlinear problem, explicit scheme is preferred than implicit approaches to avoid costly subiterations to treat the nonlinearities at each time step. Then, an important observation can be made on the physical problem. In most cases, the air transfer in the porous material is faster than moisture and heat ones. The following inequalities can be set for the Fourier numbers:

Fo a ≫ max Fo v , Fo m , Fo q .
Thus, when using an Euler explicit scheme, since we have a system of equations coupled through the advection coefficient a ⋆ v , the stability condition will be imposed by the air transfer equation. To circumvent this restriction, the diffusive air transfer equation is solved using the Du Fort-Frankel scheme. It provides an efficient highly stable stable scheme. The moisture and heat equations are advection-diffusion types. The Scharfetter-Gummel approach has shown great efficiency in preliminary studies for a single equation [5] and a system of two coupled equations [START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF]. Therefore it will be used for the spatial discretisation of the moisture and heat equations. Since this scheme has a stability condition, a two-step Runge-Kutta will be used for the time discretisation to extend the stability region of the numerical scheme [8].

A uniform discretisation is considered for space and time intervals. The discretisation parameters are denoted using ∆x for the space and ∆t for the time. The spatial cell

C = x j -1 2 , x j + 1 2 is illustrated in Figure 1(a). The discrete values of function u ( x , t )
are denoted by u n j def := u( x j , t n ) with j ∈ 1 , . . . , N and n ∈ 1 , . . . , N t . For the sake of clarity to explain the numerical method, the system Eqs. (2.14) will be written using a simpler notation and considering a linear one-dimensional problem for x ∈ 0 , 1 and t > 0 :

∂u ∂t = ∂f ∂x , (3.1a
)

∂v ∂t = ∂g ∂x + S , (3.1b
)

∂w ∂t = ∂h ∂x + Σ , (3.1c)
where

f = k 1 ∂u ∂x -a 1 u , g = k 22 ∂v ∂x -a 22 v + k 21 ∂u ∂x -a 21 u , S = -c 21 ∂u ∂t + c 2s ∂σ ∂t , h = k 33 ∂w ∂x -k 31 ∂u ∂x + a 31 u , Σ = c 31 ∂u ∂t + c 32 ∂v ∂t + c 3s ∂σ ∂t .
The term

3 i=1 ∇ c i T • j c , i in Eq. (2.14b
) is first omitted for the description of the numerical method. Its inclusion is detailed in Section 3.3.

Spatial discretisation

3.1.1 Equation (3.1a), moisture mass balance

The discretisation of Eq. (3.1a) gives the following semi-discrete difference relation:

du j dt = 1 ∆x f j+ 1 2 -f j-1 2 .
The Scharfetter-Gummel scheme was first proposed in [START_REF] Scharfetter | Large-signal analysis of a silicon Read diode oscillator[END_REF]. It assumes that the numerical flux is constant on the dual cell C ⋆ =

x j , x j+1 , which is illustrated in Figure 1(a). Additional theoretical results have been proposed recently in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws: Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF][START_REF] Gosse | Viscous equations treated with L-splines and Steklov-Poincare operator in two dimensions[END_REF]. The description of the scheme for advection-diffusion transport in porous material are provided in [5,[START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF]. Thus, for the sake of compactness, interested readers may refer to the above mentioned references. The semi-discrete difference relation for Eq. (3.1a) is directly given by:

du j dt = d 1 ∆x B Θ u j+1 -B Θ + B -Θ u j + B -Θ u j-1 . (3.2)
As mentioned in [START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF], one important advantage of the Scharfetter-Gummel scheme is the possibility of computing the flux f j+ 1 2 as wells as the exact solution of u :

u ( x ) = - 1 a 1 f j+ 1 2 + C j e a 1 x k 1 , ∀ x ∈ C ⋆ , (3.3)
where

C j def := u j -u j+1 1 -e Θ e - a 1 x j k 1
is a defined constant.

Equation (3.1b), energy balance

For the spatial discretisation of Eq. (3.1b), the semi-discrete difference relation yields:

dv j dt = 1 ∆x g j+ 1 2 -g j-1 2 + S j .
The source term S j is evaluated using Eq. (3.2):

S j = -c 21 du j dt .
By analogy, using the Scharfetter-Gummel approach, the numerical flux g j+ 1 2 is assumed constant on the dual cell C ⋆ . The exact solution of u ( x ) is already known for x ∈ C ⋆ so that the numerical flux g j+ 1 2 in the dual cell can be computed. It is important to remark that the solution considers a fully coupled approach between Equations (3.1a) and (3.1b). For the sake of notation compactness, the computation is provided in the supplementary Maple TM file. It should be noted that the exact solution v ( x ) ∈ C ⋆ can also be computed.

Equation (3.1c), moist air mass balance

The following semi-discrete difference relation stands for Eq. (3.1c):

dw j dt = 1 ∆x h j+ 1 2 -h j-1 2 + Σ j ,
where the numerical flux h j+ 1 2 is given by:

h j+ 1 2 = k 33 ∂w ∂x j+ 1 2 -k 31 ∂u ∂x j+ 1 2 + a 31 u j+ 1 2 . (3.4)
The exact solution of u ( x ) is already known so that it is possible to compute the second and the third right-hand side terms of Eq. (3.4). The hypothesis of the Scharfetter-Gummel approach, assuming the flux to be constant through the dual cell C ⋆ , is not considered for the diffusion term. This term is approximated using central finite difference scheme. At the end, the semi-discrete difference relation for Eq. (3.1c) is:

dw j dt = k 33 ∆x 2 w j+1 -2 w j + w j-1 (3.5) - k 31 a 1 ∆x k 1 C j e a 1 x j+ 1 2 k 1 -C j-1 e a 1 x j-1 2 k 1 (3.6) - a 31 a 1 f j+ 1 2 -f j-1 2 + a 31 C j e a 1 x j+ 1 2 k 1 -C j-1 e a 1 x j-1 2 k 1 + Σ j , (3.7)
where the source term Σ j = c 31 du j dt is evaluated using Eq. (3.2). For the treatment of the boundary conditions, interested reader is invited to consult [START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF] for the details.

Temporal discretisation

The purpose is to relax the stability restriction of the system of differential equations (3.1). The main restriction comes from Eq. (3.1c) and will be softened by using the Du Fort-Frankel scheme presented in Section 3.2.2. Nevertheless, it is also of major importance to relax the stability condition coming from Eqs. (3.1a) and (3.1b). For this purpose, the two-step Runge-Kutta methods can be used since it enables to provide with fewer stages the same order of accuracy than standard one-step Runge-Kutta approaches [START_REF] Jackiewicz | Two-Step Runge-Kutta Methods[END_REF]. Moreover, it is possible to extend the regions of stability of the scheme. In addition, the numerical scheme is written in an explicit way, avoiding costly sub-iterations to treat the nonlinearities observed when using implicit approaches. Thus, the two-step Runge-Kutta scheme provides a competitive option compared to the classical explicit Euler, one-step Runge-Kutta or implicit methods.

∂u ∂t = F ( u ) , F ( u ) = ∂f ∂x , ∂v ∂t = G( v ) , G( v ) = ∂g ∂x + S .
The two-step Runge-Kutta(TSRK) methods give the following fully discrete scheme [8]:

u n+1 = θ u n-1 + 1 -θ u n + ∆t s k=1 ν k F U n-1 k + µ k F U n k , U n k = λ k u n-1 + 1 -λ k u n + ∆t s p=1 a k p F U n-1 p + b k p F U n p
and similarly for v :

v n+1 = θ v n-1 + 1 -θ v n + ∆t s k=1 ν k G V n-1 k + µ k G V n k , V n k = λ k v n-1 + 1 -λ k v n + ∆t s p=1 a k p G V n-1 p + b k p G V n p ,
where θ , ν k , µ k , λ k , a kp and b kp are numerical coefficients provided in [8] depending on the value of the number of stage s . According to the vocabulary in [8], the approach has two-steps since it computes the fields u n+1 , v n+1 using the previous values u n , v n and u n-1 , v n-1 at t n and t n-1 , respectively.

Equation (3.1c), moist air mass balance

For the description, Eq. (3.5), resulting from the semi-discrete spatial discretisation of Eq. (3.1c) with the method of horizontal lines [28], is written as:

dw j dt = k 33 ∆x 2 w j+1 -2 w j + w j-1 + H j , (3.8)
where

H j = k 31 a 1 ∆x k 1 C j e a 1 x j+ 1 2 k 1 -C j-1 e a 1 x j-1 2 k 1 - a 31 a 1 f j+ 1 2 -f j-1 2 + a 31 C j e a 1 x j+ 1 2 k 1 -C j-1 e a 1 x j-1 2 k 1 + Σ j .
Using the Du Fort-Frankel approach [START_REF] Dufort | Stability Conditions in the Numerical Treatment of Parabolic Differential Equations[END_REF], the numerical scheme given for Eq. (3.8) is expressed as:

w n+1 j -w n-1 j 2 ∆t = k 33 ∆x 2 w n j+1 -w n-1 j + w n+1 j + w n j-1 + H n j , (3.9)
where the term 2 w n j has been replaced by w n+1 j

+ w n-1 j
. The stencil is illustrated in Figure 1(b). Re-arranging Eq. (3.9) to compute w n+1 j , it gives the following discrete system:

w n+1 j = 1 -δ 1 + δ w n-1 j + δ 1 + δ u n j+1 + w n j-1 + 2 ∆t 1 + δ H n j ,
where

δ def := 2 ∆t ∆x 2 k 33 .
According to the standard von Neumann stability analysis, the Du Fort-Frankel scheme is unconditionally stable [START_REF] Gasparin | Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials[END_REF][START_REF] Taylor | The stability of the Du Fort-Frankel method for the diffusion equation with boundary conditions involving space derivatives[END_REF]. Further details and examples of applications of this approach are presented in [START_REF] Gasparin | An improved explicit scheme for wholebuilding hygrothermal simulation[END_REF][START_REF] Gasparin | Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials[END_REF].

Extension to nonlinear coefficients

The physical problem (2.14) involves coefficients c , k and a varying with the fields u , v and/or w . To extend the numerical method for nonlinear coefficients, the assumption of frozen coefficients on the dual cell C ⋆ is adopted. Thus, the flux f j + 1 2 is written as:

f j+ 1 2 = k 1 , j+ 1 2 ∂u ∂x -a 1 , j+ 1 2 u ,
where the coefficients are computed as:

k 1 , j+ 1 2 = k 1 u j+ 1 2 , v j+ 1 2 , w j+ 1 2 , a 1 , j+ 1 2 = a 1 u j+ 1 2 , v j+ 1 2 , w j+ 1 2 , with u j+ 1 2 = 1 2 u j+1 + u j , v j+ 1 2 = 1 2 v j+1 + v j , w j+ 1 2 = 1 2 w j+1 + w j .
Then, the numerical flux f j± 1 2 is computed using the approach described in Section 3.1.1. By analogy, the approach is extended for the computation of the fluxes g j± 1

β i j c , i ∂v ∂x t = t n = β i j n c , i ∂v ∂x .
Then, the flux g in Equation (3.1b) is written as:

g = k 22 ∂v ∂x -a 22 + 3 i=1 β i j n c , i v + k 21 ∂u ∂x -a 21 u ,
and the whole description of the Scharfetter-Gummel approach for the spatial discretisation can be applied.

Important features of the numerical scheme

A synthesis of the schemes used for building the numerical model is provided in Table 1. Some important features of the numerical scheme should be highlighted. First, the Scharfetter-Gummel scheme is well balanced and asymptotically preserving as detailed in [START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF]. When the advection coefficient is much greater than the diffusion one, the expression of the numerical flux tends to the so-called upwind scheme. Inversely, when the diffusion coefficient is more important, the flux is approximated by central finite difference. Therefore, the flux is correct independently from the grid parameters. This feature applies also to for the numerical flux g from Eq. (3.1b).

Another important point is that we have an explicit expression for each governing equations of system (3.1). Thus, when dealing with nonlinear coefficients, the computational efforts are preserved. No sub-iterations are required at each time step, contrarily to implicit numerical schemes.

It is of major importance to comment the stability condition of the system (3.1). First of all, with an explicit Euler scheme and central-finite differences, the stability condition of the system (3.1) is:

∆t ∆x 2 2 min 1 k , l 2 min 1 j N 1 k kl , j . (3.10)
Using the Du Fort-Frankel approach, we have an unconditionally stable explicit scheme to solve Eq. (3.1c). Thus, the discretisation parameters for this equation need only to be chosen in terms of the characteristic time of the physical phenomena. The scheme ensure to compute a bounded solution. However, for Eqs. (3.1a) and (3.1b), the stability condition of the Scharfetter-Gummel scheme combined with Euler explicit approach is given [START_REF] Gosse | L-Splines and Viscosity Limits for Well-Balanced Schemes Acting on Linear Parabolic Equations[END_REF]:

∆t max 1 k , l 2 max 1 j N k kl , j max 1 j N a kl , j k kl , j tanh a kl , j ∆x 2 k kl , j -1 ∆x . (3.11)
The stability condition (3.11) is nonlinear in ∆x . It can be observed that for large space discretisation ∆x , the conditions is equivalent to ∆t C 1 ∆x [6], C 1 being a defined constant. This condition is much less restrictive than the standard approach, which is given by the standard Courant-Friedrichs-Lewy (CFL) conditions (3.10): ∆t C 2 ∆x 2 for parabolic equations, with C 2 as a defined constant. Thus, it is not necessary to use fine spatial discretisation while using the Scharfetter-Gummel scheme for Eqs. (3.1a) and (3.1b). Moreover, since the Scharfetter-Gummel scheme provides an exact solution of the fields u and v in the dual cell C ⋆ . Thus, no interpolation is required, ensuring an accurate computation of u and v in any point of the spatial domain. It should be noted that for the field w, since the Scharfetter-Gummel approach is not used, an interpolation using w j and w j+1 is required to compute w( x ) , ∀ x ∈ x j , x j+1 .

As described in Section (3.2.1), a two-step Runge-Kutta scheme is used for the temporal discretisation of Eqs. (3.1a) and (3.1b). This scheme enables to extend the stability region compared to the classic Euler approach. It can be expected a less restrictive stability condition than the one given by Eq. (3.11) is obtained.

Summarizing, the System (3.1) is composed of three coupled parabolic equations. One diffusion driven equation solved using the Du Fort-Frankel approach providing an unconditionally stable explicit scheme. Two advection-diffusion equations solved by means of the Scharfetter-Gummel approach, by implementing an explicit well balanced and asymptotically preserved numerical scheme. The stability condition ∆t ∼ ∆t , for large spatial discretisation, is relaxed using an efficient two-step Runge-Kutta approach for the temporal discretisation of these two equations. It should be noted that the Scharfetter-Gummel scheme cannot be used for the diffusion-type Equation (3.1c).

Efficiency and reliability of a numerical model

To compare the efficiency of the numerical model several criteria are considered. First, the error with a reference solution u ref ( x , t ) is evaluated. According to the definition provided in Section 2, it enables to evaluate the numerical approximations introduced by the numerical model compared to the mathematical model. The error between the solution, obtained by the numerical model, and the reference one is computed as a function of x using:

ε 2 ( x ) def := 1 N t N t n = 1 u ( x , t n ) -u ref ( x , t n ) 2 ,
where N t is the number of temporal steps. The global uniform error is measured using the norm of the space

L ∞ L 2 0 , τ , 0 , L
is given by the maximum value of ε 2 ( x ) :

ε ∞ def := sup x ∈ 0 , L ε 2 ( x ) ,
where τ is the time horizon of simulation. Another criterion is the significative correct digits of the solution, computed according to [START_REF] Söderlind | Evaluating numerical ODE/DAE methods, algorithms and software[END_REF]:

scd ( u ) def := -log 10 u ( x , τ ) -u ref ( x , τ ) u ref ( x , τ ) ∞ .
The last criterion is the computational (CPU) time required by the numerical model to compute the solution. It is measured using Matlab TM platform on a computer using Intel i7 CPU and 32 GB of RAM.

For the purposes of comparing the numerical predictions of a model with experimental observations to evaluate the reliability of the model, the relative error for temperature and vapor pressure are defined as:

ε T r ( x 0 , t ) def := T num ( x 0 , t ) -T obs ( x 0 , t ) max t T obs ( x 0 , t ) -min t T obs ( x 0 , t ) , ε P 1 r ( x 0 , t ) def := P num 1 ( x 0 , t ) -P obs 1 ( x 0 , t ) P obs 1 ( x 0 , t ) ,
where super script num. stands for the output field obtained by the numerical model and obs. for the experimental observation of the field obtained at the sensor location x 0 . The relative error is computed in a such way to avoid scaling problems due to the temperature scale.

Last, some physical hypotheses assumed in the mathematical model can be discussed after obtaining the numerical predictions. For this, the relative L 2 error is computed according to:

ε 2 ,r q ( t ) def := N x j = 1 q app ( x j , t ) -q ref ( x j , t ) 2 N x j = 1 q ref ( x j , t ) 2 .
Since the hypothesis can be computed on any physical quantity q (such as the enthalpy or volumetric moisture source among others), it is important to specify it in the definition of the error ε 2 ,r . The quantity q app designates the quantity approximated by the hypothesis and q ref the reference one.

Validation of the numerical model

After presenting the numerical model, a case study with nonlinear coefficients and Robin-type boundary conditions is considered to validate its implementation. The reference solution is computed using a numerical pseudo-spectral approach obtained with the Matlab TM open source toolbox Chebfun [12].

Description of the case

The dimensionless problem is considered for the validation:

c ⋆ m ∂u ∂t ⋆ = Fo m ∂ ∂x ⋆ k ⋆ m ∂u ∂x ⋆ -Pe m a ⋆ m u , [4pt]c ⋆ q ∂v ∂t ⋆ = Fo q ∂ ∂x ⋆ k ⋆ q ∂v ∂x ⋆ -Pe q a ⋆ q v + Fo q γ r ⋆ ∂ ∂x ⋆ k ⋆ v ∂u ∂x ⋆ -Pe v a ⋆ v u -Ko qv c ⋆ qv r ⋆ ∂u ∂t ⋆ + Ko qs c ⋆ qs r ⋆ ∂σ ∂t ⋆ -Fo q 3 i=1 β i ∂v ∂x j c , i , c ⋆ a ∂w ∂t ⋆ = Fo a ∂ ∂x ⋆ k ⋆ a ∂w ∂x ⋆ -Fo a δ ∂ ∂x ⋆ k ⋆ v ∂u ∂x ⋆ -Pe v a ⋆ v u + Ko av c ⋆ av ∂u ∂t ⋆ + Ko as c ⋆ as ∂σ ∂t ⋆ + Ko at c ⋆ at ∂v ∂t ⋆
, with the following boundary conditions at the air-material interface:

k ⋆ m ∂u ∂x ⋆ -Pe m a ⋆ m u = Bi m η u ∞ -u + η g ⋆ m , ∞ , k ⋆ q ∂v ∂x ⋆ -Pe q a ⋆ q v + γ r ⋆ k ⋆ v ∂u ∂x ⋆ -Pe v a ⋆ v u = Bi q η v ∞ -v + γ Bi v r ⋆ η u ∞ -u + η g ⋆ q , ∞ , w = w ∞ ,
where η corresponds to the projection of the surface normal unity vector n on the axis x . The simulation is performed for x ⋆ ∈ 0 , 1 and t ⋆ ∈ 0 , 24 . At the initial state, all fields are set to zero. The boundary conditions are defined as:

u ∞ , L = 0.6 sin π t 2 + sin 2 π 24 t , u ∞ , R = 0.9 sin 2 π 6 t 2 , v ∞ , L = 1.2 sin 2 π 5 t 2 + sin 2 π 24 t , v ∞ , R = 0.5 sin 2 π 3 t 2 , w( 0 , t ) = 1.3 sin 2 π 6 t 2 , w( 1 , t ) = 0.7 sin 2 π 4 t 2 + sin 2 π 24 t ,
with the following numerical values for the Biot numbers:

for x ⋆ = 0 , Bi m = 200 , Bi q = 30 , Bi v = 10 4 , for x ⋆ = 1 Bi m = 125 , Bi q = 22.5 , Bi v = 5 • 10 3 .
The numerical values of other dimensionless numbers are:

Fo m = 4 • 10 -3 , Fo q = 4 • 10 -2 , Fo a = 4 • 10 -1 , Pe m = 2 • 10 -3 , Pe v = 3 • 10 -3 , Pe q = 5 • 10 -3 , Ko qv = 1 • 10 -2 , Ko av = 2 • 10 -2 , γ = 1.5 • 10 -3 , δ = 3 • 10 -4 ,
with the following functions for the parameters c ⋆ , k ⋆ and a ⋆ :

k ⋆ q = 1 + 0.05 • u + 0.03 • v + 0.01 • w , a ⋆ q = 1 + 0.01 • v , c ⋆ q = 1 + 0.6 • u 2 + 0.1 • v + 0.3 • w , a ⋆ v = 1 + 0.01 • u , k ⋆ v = 1 + 0.1 • u + 0.01 • u 2 + 0.02 • v + 0.02 • w , c ⋆ qv = c ⋆ av = 1 + 0.06 • u , c ⋆ m = 1 + 0.05 • u 2 + 0.03 • v 2 + 0.01 • w 2 , a ⋆ m = 1 + 0.03 • u , k ⋆ m = 1 + 0.6 • u + 0.1 • v + 0.3 • w , k ⋆ a = 1 + 0.01 • u + 0.02 • u 2 + 0.07 • v 2 + 0.04 • w 2 , c ⋆ a = 1 + 0.04 • u + 0.03 • v + 0.01 • w ,
Other parameters are equal to zero and r ⋆ is set to one.

The efficiency of the numerical model proposed in Section 3 is evaluated according to the criteria defined in Section 3.5. The model is tested for different values of s ∈ 2 , 3 for the two-step Runge-Kutta method, which coefficients are given in the following tables [8]: The proposed numerical model will be compared to the reference solution. The abbreviation NM 3 corresponds to the proposed one for different values of parameter s : NM 3b stands for s = 2 and NM 3c for s = 3 . The following discretisation parameters are set ∆x ⋆ = 0.01 and ∆t ⋆ = 2 • 10 -3 . The numerical model NM 1, based on an Euler explicit approach for the temporal discretisation, is also used with ∆x ⋆ = 0.01 and

λ k a kp b kp Θ ν k µ k s , ( 4

Results

Figures 2(a), 2(c) and 2(e) show the evolution in time of the fields at x ⋆ ∈ 0.1 , 0.5 , 0.9 , which profiles are illustrated in Figures 2(b), 2(d) and 2(f). For the sake of compactness, only the results of NM 3 are presented. A perfect agreement is noted between the solutions of the numerical models based on the Scharfetter-Gummel and two-step Runge-Kutta approaches and the Chebfun one. Since the solutions are overlaid, it highlights that the numerical approximations of the numerical models are very small. As noticed in Table 2, the error with the reference solution is very satisfactory for the three numerical models NM 3b, NM 3c and NM 1. A slight difference can be observed on the error of the field w due to the increase of the time discretisation parameter between NM 1 and NM 3. The computational time required by the models to compute the solution is reported in Table 2. For the same level of accuracy, the numerical model NM 3b enables to compute the solution almost thirteen times faster than NM 1. These results are consistent since the number of operations to perform between two-step Runge-Kutta approaches for s = 2 and for s = 3 is multiplied by two, for fixed discretisation parameters. Moreover, for the same level of accuracy, the two-step Runge-Kutta scheme enables to relax the stability condition ∆t = 10 -4 (with ∆x ⋆ = 0.01) by twenty times with ∆t = 2 • 10 -3 . This validation case enhances the possibility of using the proposed numerical models for real applications with a very satisfactory accuracy and a reduced computational time compared to more conventional approaches.

Comparison of the numerical predictions with experimental data

Previous section aimed at validating the results of the numerical model with several reference solutions. The proposed numerical model showed a high accuracy with a relaxed stability condition compared to standard approaches. In other words, we verified that the differences due to numerical approximations of the mathematical model by the numerical one are minors. The numerical predictions of the numerical model are now compared with experimental data. The issue is to evaluate the physical approximations introduced by the definition of the mathematical model.

Experimental set-up

The experimental data was obtained at the University of Savoie Mont Blanc, laboratory LOCIE (Laboratory of Optimisation of the Conception and Engineering of the Environment) in the frame of the research project HYGRO-BAT [2,46]. The experimental benchmark provides data for one-dimensional heat and mass transfer in wood fiberboard.

The site is located in Le Bourget-du-Lac, France. A picture of the PASSYS cell [START_REF] Jensen | The PASSYS Project Phase 1: Subgroup Model Validation and Development : Final Report 1986-1989[END_REF] is shown in Figure 3(a). The wall is composed of a 2 cm outside coating and two 8 cm wood fiberboards, as illustrated in Figure 3(b). As mentioned in Section 2.8, at the interface between the materials, the continuities of the air, moisture and heat fluxes are assumed. The temperature and relative humidity in the cell are controlled by the air handling unit. During the first 7 days, the inside temperature and relative humidity are set to approximately 24 • C and 0.4 . After this period, the relative humidity is increased to 0.7 , while the temperature is maintained constant. The exterior boundary conditions are imposed by the climate.

Several SHT75 Sensirion sensors are evenly placed inside the material as shown in Figure 3(b). They provide temperature and relative humidity measurements with an uncertainty σ meas T = 0.3 • C and σ meas φ = 0.018 , respectively [37]. The sensors are shifted (in y and z directions) to avoid perturbations of the transfer by themselves. In addition, the sensors at x = 4 , 12 cm are inserted within the layer by perforating a hole, which is then fulfilled with wood fiber. According to this design, the uncertainty on the position of the sensors scales with σ pos = 1 cm for the sensors at x = 4 , 12 cm and σ pos = 0.5 cm for the others. The difference air pressure between the inside and outside parts of the cell is measured using Furness Controls FC0332 sensors with and uncertainty certified to σ meas P = 0.5 % . Data measurement is taken with an interval of 5 min .

For the comparison, only the experimental data in the wood fiberboard are considered. It enables to avoid additional uncertainties from the climate-driven boundary conditions, particularly wind-driven rain (that was not measured) and from the unknown surface convective coefficient. Moreover, during the research project, only the material properties Table 3. Material properties of wood fiber and physical constants.

Material properties of wood fiber. of the wood fiber have been fully determined in [34,[START_REF] Vololonirina | Characterization of hygrothermal properties of wood-based products -Impact of moisture content and temperature[END_REF]. The properties used for this case are synthesized in Table 3 with other physical constants. For the saturation pressure P sat , associating the vapor pressure P 1 and the relative humidity φ , the following expression are considered:

φ = P 1 P sat T , P sat T = P • sat T -T A T B α , (5.1a)
P • sat = 997.3 Pa , T A = 159.5 • K , T B = 120.6 • K , α = 8.275 . (5.1b)
The total sequence of experimental data corresponds to 14 days and the physical domain is defined for x ∈ 0 , L with L = 16 cm . Accordingly, for the three fields, Dirichlet boundary conditions are considered. Measurements provided by sensors at x = 0 cm and x = 16 cm are used to prescribe both temperature and vapor pressure at the boundaries. For the air pressure P , at x = 0 cm , the pressure is imposed using measurements from a differential sensor and considering the reference pressure of 1 bar . At x = 16 cm , the air pressure is assumed as constant in time and equal to the initial condition. This assumption is equivalent to consider a null air velocity at the interface between the coating and the wood fiber (x = 16 cm). This hypothesis is very likely since the air permeability of the coating is 100 times lower than It should be noted there was a long period (around 6 months) between the installation of the sensors and the presented set of experimental data. However, there are no evidences that gradients of the temperature and the vapor pressure are well established within the wall. An open question to perform numerical predictions with the proposed model is the initial condition in the wall. Here, we used the third-order polynomial interpolation of the measurements at x = 0, 4, 8, 12, 16 cm , at the first time instant. For the air pressure, the linear interpolation is supposed, assuming an established gradient of P . The initial condition of the numerical model are shown in Figure 5. An absolute error ε 2 of 5.21 Pa and 0.29 • C and a relative error error ε 2,,r of 0.08 and 0.04 are obtained between the experimental data and the interpolated initial profiles of vapor pressure and temperature, respectively. Some complementary information on the experimental design is given in [START_REF] Gasparin | An adaptive simulation of nonlinear heat and moisture transfer as a boundary value problem[END_REF].

Discussion on some hypotheses regarding the mathematical model

Before comparing the numerical predictions with the experimental observations, some possible hypotheses on the formulation of the mathematical model are discussed. The analysis is based on the relative error ε 2 ,r defined in Section 3.5 and its probability density functions computed using the kernel smoothing function for the extrapolation [START_REF] Hill | Kernel estimation of a distribution function[END_REF].

First of all, the investigations concern the constant gas law R 13 . The exact definition is provided in Eq. (2.6). Since in building physics application the dry air pressure is higher than the vapor pressure P 3 ≫ P 1 , the constant gas law can be approximated by R 13 ≈ R 3 . The probability of the relative error between both expressions is shown in Figure 6. One may conclude that for this case study an error lower than ∼ 4 % is committed using the approximation of R 13 . Another approximation may be taken in the expression of the latent heat of vaporization r 12 . The expression provided by Kelvin's law in Eq. (2.10) can be approximated by the vaporization latent heat at the reference temperature r 12 ≈ r • 12 . As shown in Figure 6, the relative error regarding the effect of temperature on r 12 is lower than 10 -2 . Thus, the approximation might be also satisfactory for this case study.

In Eq. (2.14b), the volumetric specific heat c q is defined by summing up the capacity of the dry material and that of each constituting phase (liquid, vapor and dry air). Commonly in literature, as for instance in [START_REF] Belleudy | Experimental and numerical investigations of the effects of air leakage on temperature and moisture fields in porous insulation[END_REF], the thermal capacity is approximated by the sum of the one for the dry-basis material and the one for the liquid phase c q ≈ c 0 ρ 0 + c 2 w 12 . It should be remarked that in the approximation, the liquid water content is provided by the water sorption curve w 12 . In the exact expression, the liquid water content is given by the expression of w 2 in Eq. (2.2). As noticed in Figure 6, the relative error reaches ∼ 30 % indicating that the approximation may not provide reliable results for this case study.

The significance of the two transient termsr 12 c qv ∂P 1 ∂t + r 12 c qs ∂σ ∂t in the righthand side of the heat Equation (2.14b) is studied. The relative error on the right-hand side is computed for both cases. As shown in Figure 6, the maximum error reaches ∼ 21 % with an average value of ∼ 4 % . Moreover, the probability of the relative error has an important standard deviation. It indicated that the approximation may impact locally the reliability of the numerical predictions. But overall, the average prediction of the temperature may be satisfactory.

An important hypothesis in the mathematical model concerns the expression of the volumetric vapor source as presented in Section 2.5. In [31], Luikov neglects the time variation of the vapor mass. The importance of this approximation is evaluated by comparing both literature expressions: Exact expression:

I 1 = ∂w 1 ∂t + ∇ • j c , 1 ,
Approximation from the literature:

I 1 ≈ ∇ • j c , 1 .
The relative error is shown in Figure 7(a) and is of the order of 50 % . Therefore, the importance of the transient term is non-negligible for this case. Important discrepancies are remarked in Figure 7(b) between the two approximations.

The last hypothesis involves the gradient of enthalpy in the heat Equation (2.14b). As mentioned in Section 3.3, the term 3 i=1 ∇ c i T • j c , i is included into the advective heat flux by freezing the convective flux at each time iteration. The importance of this term is investigated by comparing the two following expressions:

Exact expression: j a = a q T + 3 i=1 c i j c , i T , Approximation from the literature: j a ≈ a q T .

In Figure 6(b), a relative error of 40 % on the advective heat flux is committed with this approximation. The term with the gradient of enthalpy is particularly important when the vapor pressure increases within the wall, i.e., for t ∈ 8 , 14 days.

Results of the comparison and discussion

If the experimental data was already used in previous works [2, [START_REF] Rouchier | Hygric characterization of wood fiber insulation under uncertainty with dynamic measurements and Markov Chain Monte-Carlo algorithm[END_REF]46], only the uncertainties of the sensor measurement were presented. Here, we propose to evaluate the propagation of uncertainties due to sensor measurement and sensor location. The total uncertainty on measurement of temperature T and vapor pressure P 1 are denoted by σ T and σ P 1 , respectively. At the measurement point x 0 = 4 , 8 , 12 cm, they are computed by:

σ 2 T ( x 0 , t ) = σ meas T 2 + ∂T ∂x x = x 0 σ pos 2 , σ 2 P 1 ( x 0 , t ) = σ meas P 1 2 + ∂P 1 ∂x x = x 0 σ pos 2 ,
where, as mentioned in Section 5.1, σ pos is the sensor position uncertainty. The quantities σ meas are the sensor measurement uncertainties. For the temperature, σ meas T is directly given by the sensor manufacturer. For the vapor pressure, σ meas Eq. (5.1):

σ meas P 1 = P sat T σ meas φ + φ P sat T T -T A σ meas T ,
where σ meas φ is directly given by the sensor manufacturer. The quantities ∂T ∂x x = x 0 and ∂P 1 ∂x x = x 0 are evaluated using the results of the numerical model at the points of probes

x 0 = 4 , 8 , 12 cm.
Using the results obtained in Section 4 on the efficiency of the proposed numerical models in terms of accuracy and reduced computational time, the solution of the present case study is solved using NM 3c. The discretisation parameters are set to ∆x ⋆ = 4 • 10 -2 and ∆t ⋆ = 5 • 10 -4 , corresponding to ∆x = 6.4 mm and ∆t = 1.8 s in physical dimensions. According to the numerical values of the parameters, we have the following stability conditions: Eq. (2.14c) : ∆t 0.6 s , Eq. (2.14b) : ∆t 1.8 s , Eq. (2.14c) : ∆t 0.03 s .

Thus, the most restrictive conditions are due to to Eq. (2.14c). It can be remarked that the stability conditions are relaxed by a factor of ∼ 50 due to the efficiency of the proposed innovative method. It enables to save important computational efforts. The model needs 117 min to compute the solution with the given discretisation parameters. It corresponds to a ratio CPU time / time horizon of simulation of 8.3 min./day for the proposed numerical model and a ratio of 139.3 min./day for the standard approach based on Euler explicit scheme with central finite-differences. In other words, the NM 3c is sixteen faster to predict the physical phenomena during occurring during one day. It is also important to remark that with the discretisation parameter ∆x ⋆ = 4 • 10 -2 , the solution is computed only at 24 points in the discrete space domain. The fields at the points of observation are computed using the exact interpolation equations for solutions u , v and w as detailed in Section 3.1. Figure 8 gives an overview of the comparison of the experimental observations at x = 4 , 8 , 12 cm with the numerical predictions. The numerical predictions remain within the uncertainty range of the experimental observations. The relative error is shown in Figures 9(a) and 9(b) for vapor pressure and temperature, respectively. One may deduce there is an overall satisfactory agreement among the observations and predictions. Some moderate discrepancies can be noticed in the temperature field for t ∈ 8 , 12 days, particularly at x = 8 cm . It corresponds to the increase of the vapor pressure. It highlights that there might be uncertainties in the material properties given in Table 3 and/or that some physical phenomena are not considered in the actual model. It may arise from the interface at x = 8 cm between the two layers of wood fiberboard. In the literature [START_REF] De Freitas | Moisture migration in building walls -Analysis of the interface phenomena[END_REF][START_REF] Guimaraes | Interface influence on moisture transport in buildings[END_REF], some works question this hypothesis of flux continuity at the interface, particularly for mass transfer. During the construction of the tested wall, the intention was to obtain a good contact between the two layers. Although the experiments were designed to avoid perturbation due to the contact resistance, some uncertainties sources might still exist.

The difference air pressure P ( x , t ) -P ( L , t ) and the mass average velocity inside the material are depicted in Figure 10. The model enables to compute a velocity varying with time and space. It provides a more complex representation of the physical phenomena than the model proposed in [5,[START_REF] Berger | On the Solution of Coupled Heat and Moisture Transport in Porous Material[END_REF] where the air velocity is assumed as a constant in the whole material. Unfortunately, no experimental observations of air pressure were carried out for this case study so the reliability of the predictions cannot be evaluated. The gradient of the pressure is well established and it remains almost invariant in time. The air velocity is negative, indicating that the air flux is directed from the outside to the inside of the cell. In addition, at x = 8 cm , interface between the two layers of wood fiber, the velocity is almost constant. As mentioned before, to explain the discrepancies on the temperature at the interface, this hypothesis may be reevaluated. The introduction of sensors inside the wall, with electronic connection passing through the layers, may disturb the air leakage at this interface.

The mass flux driven by advection and diffusion are shown in Figures 11(a) and 11(b), respectively. The advection driven flux is around hundred times lower than the diffusion one. Similar observations can be made for the heat flux in Figures 11(c) and 11(d). The air transfer has a reduced impact on the heat and mass ones. It is due to a relatively low air permeability of the material k 13 = 1.1 • 10 -13 m 2 . In addition, the boundary conditions for the air pressure are almost invariant in time as can be seen in Figure 4(c). In Figures 11(c) and 11(d), it can also be remarked that the latent heat flux is smaller than the sensible one.

If the reliability of the proposed model is very satisfactory, most approaches from the literature do not consider the air transport equation (2.14c). For instance, in [START_REF] Gasparin | An adaptive simulation of nonlinear heat and moisture transfer as a boundary value problem[END_REF] the reliability is proven for the same case study but considering a model with only heat and moisture transport by diffusion processes. A natural question raises on the importance of considering air transport to better represent the whole physical phenomenon. To answer this question, the numerical predictions from the model in [START_REF] Gasparin | An adaptive simulation of nonlinear heat and moisture transfer as a boundary value problem[END_REF], denoted as HM (Heat and Moisture) model are compared to the one obtained with the proposed model, denoted as HAM (Heat, Air and Moisture) model. Figure 12 shows the probability density function of the relative error for each model between the numerical predictions and the experimental observations. If both models admit a very satisfactory reliability, the fidelity is slightly better for the proposed model. Particularly, the relative error of the predictions of temperature scales with 0.3% and with 4% for the HAM and HM models, respectively. The difference between the predictions of the models is reduced for the observations at the middle of the wall, i.e. x = 8 cm .

Conclusion

Within the context of predicting the impact of air transport on the combined heat and moisture transfer in porous building materials, this article proposes a new model with an efficient numerical scheme. In contrast to earlier proposed models in literature, the approach takes into account the transient effects of air convection without being constrained by the high numerical cost induced by the very small characteristic times of the physical phenomena. After describing in details the achievements of the mathematical model to describe the physical phenomena in Section 2, the numerical model is defined in Section 3. A system of three partial differential equations is obtained, being two of them of advective-diffusive nature, while the third one is purely diffusive. The numerical model is based on the Scharfetter-Gummel numerical scheme for the spatial discretisation of the two coupled advection-diffusion equations. It has important numerical properties as well-balanced and asymptotically preserving, ensuring the accuracy of the computed solutions. The scheme is explicit in time avoiding costly sub-iterations at each time step to handle the nonlinearities of the problem. In addition, as demonstrated in [START_REF] Gosse | L-Splines and Viscosity Limits for Well-Balanced Schemes Acting on Linear Parabolic Equations[END_REF], the CFL stability condition of the scheme ∆t ≃ ∆x is relaxed compared to standard finite-difference approaches, enabling to save extra computational efforts. For the time discretisation, a two-step Runge-Kutta approach is employed. With this scheme, the stability condition ∆t ∼ ∆x is even more relaxed. For the exclusive diffusion equation, the explicit and unconditionally stable Du Fort-Frankel scheme is used. It circumvents the strong restrictions on the stability of the numerical model coming from the very small time characteristics of the air transfer. Finally, an efficient numerical scheme is proposed to compute an accurate solution. The computational efforts are saved since the restrictions on the choice of the time discretisation steps are weakened by the use of robust numerical schemes. The features of the numerical model are enhanced in Section 4 with one case study, where all coefficients are nonlinear and the Robin-type boundary conditions are applied. The reference solution is computed using a pseudo-spectral approach. Results show a very accurate solution with a reduced computational time.

After validating the efficiency of the numerical model, Section 5 intends to evaluate the reliability of the numerical predictions by comparing them to experimental observations. The latter are achieved by measuring temperature and vapor pressure in a wall composed of two layers of wood fiberboard. The internal boundary conditions are controlled, with a sudden isothermal increase on the vapor pressure imposed by an air handling unit. For the external side, the boundary conditions are driven by climate variations. There is a very satisfactory agreement between the numerical predictions and the experimental observations. The results indicate a high reliability of the numerical model to represent the physical phenomena. The stability condition of the problem is reduced by an order of 50 with the innovative numerical model, providing a reduction of the computational effort by 16 compared to standard approach based on Euler explicit scheme with central finite differences.

Further research should be carried out on the characterization of the conditions at the interface of two materials. As noticed in Section 5, some small discrepancies were observed at the interface between two layers. Some works in the literature suggest how to modify the conditions in the model [START_REF] De Freitas | Moisture migration in building walls -Analysis of the interface phenomena[END_REF][START_REF] Guimaraes | Interface influence on moisture transport in buildings[END_REF] but the resistances at the interfaces between porous materials is still an open question. To conclude, we can say the air transport has been considered important since the HAMSTAD project (2002) [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], but a little research has been presented in the literature regarding the importance of advection in porous material and on the calculation of the velocity and pressure profiles. Therefore, to reinforce the research in the field, it is important to promote on one hand detailed experimental research to include, e.g., pressure differences along the thickness as well as microscopic simulations using methods such as the Lattice Boltzmann model including energy and momentum balances and phase change terms.

Nomenclature

Figure 1 .

 1 Figure 1. Stencils of the Scharfetter-Gummel(a) and Du Fort-Frankel(b) schemes.

  3.2.1 Equations (3.1a) and (3.1b), moisture mass and energy balances For the description of the temporal discretisation, Eqs. (3.1a) and (3.1b) are written in the following form:

β

  i j c , i ∂v ∂xis taken into account by freezing the flux j n c , i at each time step t n :

Figure 2 .

 2 Figure 2. (a, c, e) Time evolution and (b, d, f ) profiles at t ⋆ = 8 , 16 , 24 of the fields u , v and w .

Sorption curve kg/m 3 wAir permeability m 2 k 13 = 1 . 1 • 10 - 13 Porosity - Π = 0. 9 6

 31311101396 12 = 70.63 φ 3 -73.6 φ 2 + 41.05 φ + 0.26 Vapor permeability s k 1 = 3.28 • 10 -11 + 4.85 • 10 -11 φ Liquid permeability s k 2 = 0 Specific heat J/(kg.K) c 0 = 1103 Dry-basis specific mass kg/m 3 ρ 0 = 146 Physical constants Vapor gas constant J/(kg.K) R 1 = 462 Dry air gas constant J/(kg.K) R 3 = 287 Latent heat of vaporization J/kg r 12 = 2.5 • 10 Vapor specific heat J/(kg.K) c 1 = 1870 Liquid specific heat J/(kg.K) c 2 = 4180 Dry air specific heat J/(kg.K) c 3 = 1006 Dynamic viscosity Pa.s µ = 1.8 • 10 -5

Figure 3 .

 3 Figure 3. (a) Picture of the PASSYS cell and (b) illustration of the wall with the implementation of the sensors.

Figure 4 .

 4 Figure 4. Measured boundary conditions used for the numerical model.

Figure 5 .

 5 Figure 5. Measured and interpolated initial conditions used for the numerical model.

T and σ meas P 1

 1

Figure 6 .

 6 Figure 6. Probability density function of the relative error for some hypothesis in the formulation of the mathematical model (a) for the latent heat of vaporization r 12 , for the volumetric heat capacity c q and for the transient term in the right hand side of the heat Equation (2.14b) and (b) for the for the constant gas law R 13 and the gradient of enthalpy.

Figure 7 .

 7 Figure 7. Importance of the approximation of the volumetric vapor source I 1 .

Figure 8 .

 8 Figure 8. Comparison of the experimental observations with the numerical predictions for vapor pressure (a, c, e) and temperature (b, d, f ). The grey shadows represent the measurement uncertainties.

Figure 9 .Figure 10 .

 910 Figure 9. Relative error between the experimental observations and the numerical predictions for vapor pressure (a) and temperature (b).

Figure 11 .

 11 Figure 11. Time evolution of the mass flux by advection (a) and diffusion (b). Time evolution of the heat flux by advection (c) and diffusion (d).

Figure 12 .

 12 Figure 12. Comparison of the probability density function of the relative error between the experimental observations and the numerical predictions for vapor pressure (a,c,e) and temperature (b,d,f ), obtained with the proposed model and the purely diffusive one.

Table 1 .

 1 Synthesis of the numerical scheme used to build the numerical model.

	Equation	Eqs. (3.1a) and (3.1b)	Eq. (3.1c)
	Type	Advection-diffusion	Diffusion
	Space discr.	Scharfetter-Gummel	Central Finite-Differences
	Time discr.	Two-step Runge-Kutta	Du Fort-Frankel
	Features	Explicit, extended stability conditions ∆t ∼ ∆x Explicit, unconditionally stable

Table 2 .

 2 Computer run time required by the numerical models to solve the case study (Section 4.1), with t 0 = 136 s , and variation of the error ε ∞ of the solutions u , v and w. • 10 -3 0.08 • t 0 7.36 • 10 -3 4.76 • 10 -3 6.48 • 10 -3 NM 3c 10 -2 2 • 10 -3 0.10 • t 0 7.35 • 10 -3 4.76 • 10 -3 6.48 • 10 -3 ∆t ⋆ = 10 -4 , which is consistent with the CFL condition. It can be noted that for higher time discretisation (with ∆t ⋆ = 0.01), the model NM 1 does not compute a bounded solution. The reference solution is the pseudo-spectral one obtained with the Matlab TM open source toolbox Chebfun [12].

	Numerical Model (NM) ∆x	∆t	CPU time ε ∞ for u	ε ∞ for v	ε ∞ for w
	NM 1	10 -2 10 -4	t 0	7.35 • 10 -3 4.75 • 10 -3 2.1 • 10 -3
	NM 3b	10 -2 2

  /m 2 ] k 1 , k 2 , k v , k m vapor, liquid, moisture permeability [s] P , P 1 , P 2 , P 3 pressure [Pa] R 1 , R 3 , R 13 specie gaz constant [J/(kg • K)]

		Greek letters
	α m surface vapour transfer coefficient [s/m]
	α q surface heat transfer coefficient	[W/(m 2 • K)]
	φ relative humidity	[-]
	ρ specific mass	[kg/m 3 ]
	Π material porosity	[-]
	σ saturation rate	[-]
	µ dynamic viscosity	[Pa • s]
		Latin letters
	a q	heat advection coefficient	[W/(m 2 • K)]
	a v	vapor advection coefficient	[s/m]
	c	specific heat		[J/(kg • K)]
	c at	storage coefficient		[kg • Pa/(J • K)]
	c m	moisture storage capacity	[kg/(m 3 • Pa)]
	c q	volumetric heat capacity	[J/(m 3 • K)]
	c qs , c as	storage coefficient		[kg • Pa/J]
	c qv , c a , c av	storage coefficient		[kg/J]
	g	liquid flux		[kg/(s • m 2 )]
	h	specific enthalpy		[J/kg]
	I	volumetric capacity of source/sink	[kg/(m 3 .s)]
	j	mass flow		[kg/(s • m 2 )]
	j q	heat flow		[W/m 2 ]
	k a [s 2 k q permeability coefficient thermal conductivity [W/(m • K)]
	k 13	intrinsic permeability		[m 2 ]
	L	length		[m]
	m	mass		[kg]
	r 12	latent heat of evaporation	[J/kg]
	T	temperature		[K]
	t	time coordinate		[s]
	x	space coordinate		[m]
	V	volume		[m 3 ]
	v	mass average velocity		[m/s]
	w	volumetric concentration	[kg/m 3 ]

Acknowledgments

The authors acknowledge the Junior Chair Research program "Building performance assessment, evaluation and enhancement" from the University of Savoie Mont Blanc in collaboration with The French Atomic and Alternative Energy Center (CEA) and Scientific and Technical Center for Buildings (CSTB). The authors thanks the grants from the Ministry of Education and Science of the Republic of Kazakhstan and the visiting professor grants from the University of Savoie Mont Blanc. The authors also acklowledge the French and Brazilian agencies for their financial supports through the project CAPES-COFECUB, as weel as the CNPQ of the Brazilian Ministry of Education and of the Ministry of Science, Technology and Innovation, respectively, for co-funding.