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

TL;DR

This paper introduces a comprehensive mathematical and numerical model for simulating heat, air, and moisture transfer in porous building materials, achieving accurate results with significantly reduced computational effort.

Contribution

The work develops an innovative numerical scheme combining Du Fort-Frankel, Runge-Kutta, and Scharfetter-Gummel methods to efficiently simulate coupled heat, air, and moisture transfer.

Findings

Model predicts heat, air, moisture transfer accurately.

Simulation is approximately 16 times faster than standard methods.

Validated against experimental data with good agreement.

Abstract

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.

Tables

Table 1. Table 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 $\Delta t\sim \Delta x$	Explicit, unconditionally stable

Table 2. Table 2. Computer run time required by the numerical models to solve the case study (Section  4.1 ), with t  0 = 136 ​ [ 𝗌 ] subscript 𝑡  0 136 delimited-[] 𝗌 t_{\,0}\ =\ 136\ \bigl{[}\,\mathsf{s}\,\bigr{]} , and variation of the error ε ∞ subscript 𝜀 \varepsilon_{\,\infty} of the solutions u 𝑢 u\, , v 𝑣 v and w 𝑤 w .

Numerical Model (NM)	$\Delta x$	$\Delta t$	CPU time	${\epsilon }_{\infty }$ for $u$	${\epsilon }_{\infty }$ for $v$	${\epsilon }_{\infty }$ for $w$
NM $1$	${10}^{-2}$	${10}^{-4}$	$t_{\mathrm{\hspace{0.17em}0}}$	$7.35\cdot {10}^{-3}$	$4.75\cdot {10}^{-3}$	$2.1\cdot {10}^{-3}$
NM $3$b	${10}^{-2}$	$2\cdot {10}^{-3}$	$0.08\cdot t_{\mathrm{\hspace{0.17em}0}}$	$7.36\cdot {10}^{-3}$	$4.76\cdot {10}^{-3}$	$6.48\cdot {10}^{-3}$
NM $3$c	${10}^{-2}$	$2\cdot {10}^{-3}$	$0.10\cdot t_{\mathrm{\hspace{0.17em}0}}$	$7.35\cdot {10}^{-3}$	$4.76\cdot {10}^{-3}$	$6.48\cdot {10}^{-3}$

Table 3. Table 3. Material properties of wood fiber and physical constants.

Material properties of wood fiber.
Sorption curve $\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{\hspace{0.17em}3}}\right]$	$w_{\mathrm{\hspace{0.17em}12}}=70.63{\varphi }^{\mathrm{\hspace{0.17em}3}}-73.6{\varphi }^{\mathrm{\hspace{0.17em}2}}+41.05\varphi +0.26$
Vapor permeability $\left[𝗌\right]$	$k_{\mathrm{\hspace{0.17em}1}}=3.28\cdot {10}^{-11}+4.85\cdot {10}^{-11}\varphi$
Liquid permeability $\left[𝗌\right]$	$k_{\mathrm{\hspace{0.17em}2}}=0$
Air permeability $\left[{𝗆}^{\mathrm{\hspace{0.17em}2}}\right]$	$k_{\mathrm{\hspace{0.17em}13}}=1.1\cdot {10}^{-13}$
Porosity $[-]$	$\Pi =0.9$
Specific heat $[𝖩/(\mathrm{𝗄𝗀}.𝖪)]$	$c_{\mathrm{\hspace{0.17em}0}}=1103$
Dry–basis specific mass $\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{\hspace{0.17em}3}}\right]$	${\rho }_{\mathrm{\hspace{0.17em}0}}=146$
Physical constants
Vapor gas constant $[𝖩/(\mathrm{𝗄𝗀}.𝖪)]$	$R_{\mathrm{\hspace{0.17em}1}}=462$
Dry air gas constant $[𝖩/(\mathrm{𝗄𝗀}.𝖪)]$	$R_{\mathrm{\hspace{0.17em}3}}=287$
Latent heat of vaporization $\left[𝖩/\mathrm{𝗄𝗀}\right]$	$r_{\mathrm{\hspace{0.17em}12}}=2.5\cdot {10}^{\mathrm{\hspace{0.17em}6}}$
Vapor specific heat $[𝖩/(\mathrm{𝗄𝗀}.𝖪)]$	$c_{\mathrm{\hspace{0.17em}1}}=1870$
Liquid specific heat $[𝖩/(\mathrm{𝗄𝗀}.𝖪)]$	$c_{\mathrm{\hspace{0.17em}2}}=4180$
Dry air specific heat $[𝖩/(\mathrm{𝗄𝗀}.𝖪)]$	$c_{\mathrm{\hspace{0.17em}3}}=1006$
Dynamic viscosity $\left[\mathrm{𝖯𝖺}.𝗌\right]$	$\mu =1.8\cdot {10}^{-5}$

Equations

$$\displaystyle\mathsf{w}_{\,12}\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\frac{m_{\,12}}{V}\ =\ \mathsf{f}\,\bigl{(}\,\phi\,\bigr{)}\,,$$

$$\displaystyle\sigma\ =\ \frac{1}{\rho_{\,2}\,R_{\,1}\,T\ -\ P_{\,1}}\Biggl{(}\,\frac{R_{\,1}\,T}{\Pi}\,\mathsf{w}_{\,12}\ -\ P_{\,1}\,\Biggr{)}\,.$$

$${\mathsf{w}}_1$$

$$\displaystyle\mathsf{w}_{\,2}\ =\ \mathsf{w}_{\,12}\ -\ \mathsf{w}_{\,1}\ =\ \mathsf{w}_{\,12}\ -\ \frac{P_{\,1}}{R_{\,1}\,T}\ \Pi\ \bigl{(}\,1\ -\ \sigma\,\bigr{)}\,.$$

$$\frac{\partial {\mathsf{w}}_i}{\partial t}=-\nabla \cdot {\mathbf{j}}_{c,i}+I_i,$$

$$I_1+I_2=0.$$

$${\mathbf{j}}_c={\mathbf{j}}_d+{\mathbf{j}}_a.$$

$$\displaystyle\frac{\partial}{\partial t}\Bigl{(}\,\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,2}\,\Bigr{)}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\Bigl{(}\,\boldsymbol{j}_{\,c\,,\,1}\ +\ \boldsymbol{j}_{\,c\,,\,2}\,\Bigr{)}\,.$$

$$\frac{\partial {\mathsf{w}}_{12}}{\partial t}=\frac{\partial {\mathsf{w}}_{12}}{\partial \varphi }\frac{1}{P_{\mathrm{s}}}\frac{\partial P_1}{\partial t},$$

$$\displaystyle\frac{\partial\mathsf{w}_{\,12}}{\partial\phi}\,\frac{1}{P_{\,\mathrm{s}}}\,\frac{\partial P_{\,1}}{\partial t}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\Bigl{(}\,\boldsymbol{j}_{\,c\,,\,1}\ +\ \boldsymbol{j}_{\,c\,,\,2}\,\Bigr{)}\,.$$

$$\displaystyle\frac{\partial}{\partial t}\Bigl{(}\,\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,3}\,\Bigr{)}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\Bigl{(}\,\boldsymbol{j}_{\,c\,,\,1}\ +\ \boldsymbol{j}_{\,c\,,\,3}\,\Bigr{)}\ +\ I_{\,1}\,,$$

$${\mathsf{w}}_1+{\mathsf{w}}_3=\frac{P}{R_{13}T}\Pi (1-\sigma ),$$

$$\displaystyle R_{\,13}\ =\ \frac{R_{\,3}}{1\ -\ \frac{P_{\,1}}{P_{\,3}}\,\Bigl{(}\,1\ -\ \nicefrac{{R_{\,3}}}{{R_{\,1}}}\Bigr{)}}\,.$$

$$\displaystyle\partial\,\bigl{(}\,\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,3}\,\bigr{)}\ =\ \frac{\Pi\,(\,1\ -\ \sigma\,)}{R_{\,13}\,T}\,\partial P\ +\ \frac{\Pi\,P}{R_{\,13}\,T}\,\partial(\,1\ -\ \sigma\,)\ -\ \frac{\Pi\,(\,1\ -\ \sigma\,)\,P}{R_{\,13}\,T^{\,2}}\,\partial T\,.$$

$$\displaystyle\frac{\Pi\,(\,1\ -\ \sigma\,)}{R_{\,13}\,T}\,\frac{\partial P}{\partial t}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\Bigl{(}\,\boldsymbol{j}_{\,c\,,\,1}\ +\ \boldsymbol{j}_{\,c\,,\,3}\,\Bigr{)}\ +\ I_{\,1}\ +\ \frac{\Pi\,P}{R_{\,13}\,T}\,\frac{\partial\sigma}{\partial t}\ +\ \frac{\Pi\,(\,1\ -\ \sigma\,)\,P}{R_{\,13}\,T^{\,2}}\,\frac{\partial T}{\partial t}\,.$$

$$\displaystyle\frac{\partial}{\partial t}\biggl{(}\,h_{\,0}\,\rho_{\,0}\ +\ \sum_{i=1}^{3}\,h_{\,i}\,\mathsf{w}_{\,i}\,\biggr{)}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\Bigl{(}\,\boldsymbol{j}_{\,q}\ +\ \sum_{i=1}^{3}\,h_{\,i}\,\boldsymbol{j}_{\,c\,,\,i}\,\Bigr{)}\,,$$

$$\biggl{(}\,c_{\,0}\,\rho_{\,0}\ +\ \sum_{i=1}^{3}\,c_{\,i}\,\mathsf{w}_{\,i}\,\biggr{)}\frac{\partial T}{\partial t}\ +\ \sum_{i=1}^{3}\,h_{\,i}\,\frac{\partial\mathsf{w}_{\,i}}{\partial t}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\boldsymbol{j}_{\,q}\\ \ -\ \sum_{i=1}^{3}\,h_{\,i}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\boldsymbol{j}_{\,c\,,\,i}\ -\ \sum_{i=1}^{3}\,\boldsymbol{\nabla}\,\bigl{(}\,h_{\,i}\,\bigr{)}\boldsymbol{\cdot}\boldsymbol{j}_{\,c\,,\,i}\,.$$

$$\sum _{i=1}^3{h}_i\frac{\partial {\mathsf{w}}_i}{\partial t}=\text{-}\sum _{i=1}^3{h}_i\nabla \cdot {\mathbf{j}}_{c,i}+h_1{I}_1+h_2{I}_2.$$

$$\displaystyle\biggl{(}\,c_{\,0}\,\rho_{\,0}\ +\ \sum_{i=1}^{3}\,c_{\,i}\,\mathsf{w}_{\,i}\,\biggr{)}\frac{\partial T}{\partial t}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\boldsymbol{j}_{\,q}\ -\ h_{\,1}\,I_{\,1}\ -\ h_{\,2}\,I_{\,2}\ -\ \sum_{i=1}^{3}\,\boldsymbol{\nabla}\,\bigl{(}\,c_{\,i}\,T\,\bigr{)}\boldsymbol{\cdot}\boldsymbol{j}_{\,c\,,\,i}\,.$$

$$\displaystyle r_{\,12}\ =\ r_{\,12}^{\,\circ}\ +\ \bigl{(}\,c_{\,1}\ -\ c_{\,2}\,\bigr{)}\cdot\bigl{(}\,T\ -\ T_{\,\mathrm{ref}}\,\bigr{)}\ -\ R_{\,1}\,T\,\ln\,\phi\,.$$

$$\displaystyle\biggl{(}\,c_{\,0}\,\rho_{\,0}\ +\ \sum_{i=1}^{3}\,c_{\,i}\,\mathsf{w}_{\,i}\,\biggr{)}\frac{\partial T}{\partial t}\ =\ -\,\boldsymbol{\nabla}\boldsymbol{\cdot}\,\boldsymbol{j}_{\,q}\ -\ r_{\,12}\ I_{\,1}\ -\ \sum_{i=1}^{3}\,\boldsymbol{\nabla}\,\bigl{(}\,c_{\,i}\,T\,\bigr{)}\boldsymbol{\cdot}\boldsymbol{j}_{\,c\,,\,i}\,.$$

$$I_1=\nabla \cdot j_{c,1}.$$

$$\displaystyle I_{\,1}\ =\ -\,\biggl{(}\,\frac{\partial\mathsf{w}_{\,2}}{\partial t}\ +\ \boldsymbol{\nabla}\boldsymbol{\cdot}\,j_{\,c\,,\,2}\,\biggr{)}\,.$$

$$I_1=\frac{\partial {\mathsf{w}}_1}{\partial t}+\nabla \cdot j_{c,1},$$

$$\displaystyle\mathsf{w}_{\,1}\ =\ \frac{\Pi\,\,\bigl{(}\,1\ -\ \sigma\,\bigr{)}}{R_{\,1}\,T}\,.$$

$$\partial {\mathsf{w}}_1=\frac{\Pi (1-\sigma )}{R_{1}T}\partial P_1+\frac{\Pi P_1}{R_{1}T}\partial (1-\sigma )-\frac{\Pi (1-\sigma )P_1}{R_1{T}^2}\partial T.$$

$$I_1=\frac{\Pi (1-\sigma )}{R_{1}T}\frac{\partial P_1}{\partial t}-\frac{\Pi P_1}{R_{1}T}\frac{\partial \sigma }{\partial t}-\frac{\Pi (1-\sigma )P_1}{R_1{T}^2}\frac{\partial T}{\partial t}+\nabla \cdot j_{c,1},$$

$$\mathsf{v}=-\frac{k_{13}}{\mu }\nabla P,$$

$$\displaystyle\boldsymbol{j}_{\,c\,,\,1}\ +\ \boldsymbol{j}_{\,c\,,\,3}\ =\ \frac{\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,3}}{\Pi\,\bigl{(}\,1\ -\ \sigma\,\bigr{)}}\,\boldsymbol{\mathsf{v}}\,,$$

$$\displaystyle\boldsymbol{j}_{\,c\,,\,1}\ +\ \boldsymbol{j}_{\,c\,,\,3}\ =\ -\,\frac{\bigl{(}\,\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,3}\,\bigr{)}}{\Pi\,\bigl{(}\,1\ -\ \sigma\,\bigr{)}}\,\frac{k_{\,13}}{\mu}\,\boldsymbol{\nabla}P\,.$$

Full text

Key words and phrases: heat and mass transfer; porous material; Benchmarking with experimental data; Scharfetter–Gummel numerical scheme; Du Fort–Frankel numerical scheme

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Formulation of the physical phenomena

3. 2.1 Porous medium basics

4. 2.2 Moisture mass balance

5. 2.3 Moist air mass balance

6. 2.4 Energy balance

7. 2.5 Volumetric vapor source

8. 2.6 Expression of the fluxes

9. 2.7 Governing Equations

10. 2.8 Boundary, interface and initial conditions

11. 3 Elaborating an efficient numerical model

12. 3.1 Spatial discretisation

13. 3.1.1 Equation (3.1a), moisture mass balance

14. 3.1.2 Equation (3.1b), energy balance

15. 3.1.3 Equation (3.1c), moist air mass balance

16. 3.2 Temporal discretisation

17. 3.2.1 Equations (3.1a) and (3.1b), moisture mass and energy balances

18. 3.2.2 Equation (3.1c), moist air mass balance

19. 3.3 Extension to nonlinear coefficients

20. 3.4 Important features of the numerical scheme

21. 3.5 Efficiency and reliability of a numerical model

22. 4 Validation of the numerical model

23. 4.1 Description of the case

24. 4.2 Results

25. 5 Comparison of the numerical predictions with experimental data

26. 5.1 Experimental set–up

27. 5.2 Discussion on some hypotheses regarding the mathematical model

28. 5.3 Results of the comparison and discussion

29. 6 Conclusion

30. Acknowledgments

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 [7]. 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 [32]. 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 [43, 26, 10].

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 [1, 33], the commercial COMSOL™ 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 [40] or [30] for one-dimensional transfer. The former references uses an explicit scheme provided by the commercial COMSOL™ 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 [3, 4] for the simulation of two–dimensional transfer in building structures. It is also based on an explicit scheme from COMSOL™ 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 [36] with very recent theoretical results in [19, 18]. In the context of building porous media, it is successfully applied in [5] to water transport and then in [6] 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 \mathsf{w}_{\,i}\ \bigl{[}\,\mathsf{kg/m^{\,3}}\,\bigr{]} of specie(s) $i\,$, which varies in space and time, is defined by \mathsf{w}_{\,i}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \lim\limits_{V\,\rightarrow\,0}\dfrac{m_{\,i}\,(\,V\,\bigr{)}}{V}\,. We denot by m_{\,i}\,(\,V\,\bigr{)}\ \bigl{[}\,\mathsf{kg}\,\bigr{]} the mass of the specie(s) in the representative volume V\ \bigl{[}\,\mathsf{m^{\,3}}\,\bigr{]}\,. 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_{\,\mathrm{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:

[TABLE]

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 [39], 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}\ \ll\ 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 \sigma\ \bigl{[}\,\mathsf{-}\,\bigr{]} 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}\ \approx\ m_{\,2}$ is considered, the saturation rate is evaluated directly by the sorption curve. Presently, the saturation rate is given by:

[TABLE]

Assuming the perfect gas law, the volumetric vapor and dry air content are given by:

[TABLE]

Last, the liquid volumetric mass content is given by:

[TABLE]

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:

[TABLE]

where I_{\,i}\ \bigl{[}\,\mathsf{kg/(m^{\,3}.\,s)}\,\bigr{]} 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:

[TABLE]

The flow of mass is denoted by \boldsymbol{j}_{\,c\,,\,i}\ \bigl{[}\,\mathsf{kg/(m^{\,2}.\,s)}\,\bigr{]}\,. Due to the airflow occurring through the porous medium, the mechanism of mass transfer is driven by both diffusion $\boldsymbol{j}_{\,d}$ and advection $\boldsymbol{j}_{\,a}$ flows:

[TABLE]

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\ \leftarrow\ 1$) and Eq. (2.3) ($i\ \leftarrow\ 2$), the moisture transfer equation is obtained:

[TABLE]

The term $\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,2}$ is equal to the moisture content $\mathsf{w}_{\,12}\,$. Using the expression of the sorption isotherm curve Eq. (2.1), the partial differential of the moisture content becomes:

[TABLE]

where $P_{\,\mathrm{s}}$ is the saturation pressure, depending on temperature, and \displaystyle\frac{\partial}{\partial\phi}\,\Bigl{(}\,\mathsf{w}_{\,1}\ +\ \mathsf{w}_{\,2}\,\Bigr{)} is the derivative of the sorption curve of the material. In this way, the moisture transfer equation becomes:

[TABLE]

Moist air mass balance

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

[TABLE]

where $\mathsf{w}_{\,1}\ +\ \mathsf{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$:

[TABLE]

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:

[TABLE]

Assuming the porosity to be invariant, the following equation is obtained by computing the differential of Eq. (2.5):

[TABLE]

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\ -\ \sigma$ and of the temperature $T\,$. Thus, Eq.(2.4) becomes:

[TABLE]

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\ \bigl{[}\,\mathsf{J/kg}\,\bigr{]} equals the divergence of the conduction and enthalpy flux:

[TABLE]

where \rho_{\,0}\ \bigl{[}\,\mathsf{kg/m^{\,3}}\,\bigr{]} is the dry material density and c_{\,i}\ \bigl{[}\,\mathsf{J/(kg.K)}\,\bigr{]} the specific heat of each species. The expression of the volumetric water content $\mathsf{w}_{\,i}$ is detailed in Section 2.1. By assuming a constant volume, Eq (2.8) becomes:

[TABLE]

Since $I_{\,3}\ =\ 0$, by summing Eq. (2.3) ($i\ \leftarrow\ 1$), multiplied by $h_{\,1}$, Eq. (2.3) ($i\ \leftarrow\ 2$), multiplied by $h_{\,2}$ and Eq. (2.3) ($i\ \leftarrow\ 3$), multiplied by $h_{\,3}$, we obtain that:

[TABLE]

Thus, Eq. (2.9) becomes:

[TABLE]

We denote by r_{\,12}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}h_{\,1}\ -\ h_{\,2}\ \bigl{[}\,\mathsf{J/kg}\,\bigr{]} as the latent heat of vaporization. According to Kelvin’s law, it is defined as [39]:

[TABLE]

Finally, the energy conservation equation yields:

[TABLE]

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

[TABLE]

Similar approach has been adopted in [29, 40]. This assumption was discussed in [39], for the case of rammed earth material, by comparing the results using expression (2.11) and the following one:

[TABLE]

The latter expression assumes that the moisture concentration $\mathsf{w}_{\,2}$ is obtained by the sorption isotherm $w_{\,12}\ \approx\ w_{\,2}\,$. This assumption is not considered in this work, a third expression of the source is proposed:

[TABLE]

where the vapor concentration $\mathsf{w}_{\,1}$ is expressed as:

[TABLE]

The differential of $\mathsf{w}_{\,1}$ gives:

[TABLE]

Then, we obtain the final expression of the volumetric source term $I_{\,1}\,$:

[TABLE]

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, $\boldsymbol{\mathsf{v}_{\,3}}\ =\ \boldsymbol{\mathsf{v}_{\,2}}\ =\ \boldsymbol{\mathsf{v}}\,$, and can be expressed by the Darcy law, remembering that the gravity effects are neglected:

[TABLE]

where k_{\,13}\ \bigl{[}\,\mathsf{m^{\,2}}\,\bigr{]} is the intrinsic permeability of the material and $\mu$ 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:

[TABLE]

which can be rewritten using Eq. (2.13) for the velocity:

[TABLE]

The vapor flow is driven by diffusion and advection:

[TABLE]

where $k_{\,1}$ is the vapor permeability of the material and R_{\,1}\ \bigl{[}\,\mathsf{J/(kg.K)}\,\bigr{]} 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:

[TABLE]

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}$ [6]. Thus, the total moisture flux yields:

[TABLE]

The heat flux is expressed as:

[TABLE]

It is driven by conduction and advection of the moist air.

Governing Equations

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

[TABLE]

along with the storage coefficients:

[TABLE]

and the permeability coefficients:

[TABLE]

Finally, the set of partial differential governing equations for the physical problem of heat, air and moisture transfer through porous building material becomes:

[TABLE]

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:

[TABLE]

where g_{\,m\,,\infty}\ \bigl{[}\,\mathsf{kg/(m^{\,2}.s)}\,\bigr{]} 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:

[TABLE]

where \alpha_{\,m}\ \bigl{[}\,\mathsf{s/m}\,\bigr{]} is the surface vapor transfer coefficient and $\boldsymbol{n}$ is the surface normal unitary vector. Therefore, for the moisture transfer equation, the following boundary conditions are stated:

[TABLE]

Similarly, for the heat transfer, the flow is proportional to the surface and ambient conditions. An additional heat flux $\boldsymbol{g}_{\,q\,,\,\infty}$ is introduced to take into account the radiative and convective heat transfer:

[TABLE]

Therefore, the boundary condition for the heat transfer equation can be written as:

[TABLE]

Finally, the air pressure at the interface between the material and the ambient air is imposed:

[TABLE]

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 [9, 20]:

[TABLE]

As initial conditions, the fields may depend on the space coordinate $x\,$:

[TABLE]

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:

[TABLE]

Thus, when using an Euler explicit scheme, since we have a system of equations coupled through the advection coefficient $a_{\,v}^{\,\star}\,$, 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 [6]. 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 $\Delta x$ for the space and $\Delta t$ for the time. The spatial cell \mathcal{C}\ =\ \bigl{[}\,x_{\,j\ -\ \frac{1}{2}}\,,\,x_{\,j\ +\ \frac{1}{2}}\,\bigr{]} is illustrated in Figure 1(a). The discrete values of function $u\,(\,x\,,\,t\,)$ are denoted by $u_{\,j}^{\,n}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ u(\,x_{\,j}\,,\,t^{\,n}\,)$ with j\ \in\ \bigl{\{}\,1\,,\ldots\,,N\,\bigr{\}} and n\ \in\ \bigl{\{}\,1\,,\ldots\,,N_{\,t}\,\bigr{\}}\,.

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\ \in\ \bigl{[}\,0\,,\,1\,\bigr{]} and $t\ >\ 0\,$:

[TABLE]

where

[TABLE]

The term \displaystyle\sum_{i=1}^{3}\,\boldsymbol{\nabla}\,\bigl{(}\,c_{\,i}\,T\,\bigr{)}\boldsymbol{\cdot}\boldsymbol{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:

[TABLE]

The Scharfetter–Gummel scheme was first proposed in [36]. It assumes that the numerical flux is constant on the dual cell \mathcal{C}^{\,\star}\ =\ \Bigl{[}\,x_{\,j}\,,\,x_{\,j+1}\,\Bigr{]}\,, which is illustrated in Figure 1(a). Additional theoretical results have been proposed recently in [17, 18]. The description of the scheme for advection–diffusion transport in porous material are provided in [5, 6]. 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:

[TABLE]

As mentioned in [6], one important advantage of the Scharfetter–Gummel scheme is the possibility of computing the flux $f_{\,j+\frac{1}{2}}$ as wells as the exact solution of $u\,$:

[TABLE]

where

[TABLE]

is a defined constant.

3.1.2  Equation (3.1b), energy balance

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

[TABLE]

The source term $S_{\,j}$ is evaluated using Eq. (3.2):

[TABLE]

By analogy, using the Scharfetter–Gummel approach, the numerical flux $g_{\,j+\frac{1}{2}}$ is assumed constant on the dual cell $\mathcal{C}^{\,\star}\,$. The exact solution of $u\,(\,x\,)$ is already known for $x\,\in\,\mathcal{C}^{\,\star}$ so that the numerical flux $g_{\,j+\frac{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™ file. It should be noted that the exact solution $v\,(\,x\,)\,\in\,\mathcal{C}^{\,\star}$ can also be computed.

3.1.3  Equation (3.1c), moist air mass balance

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

[TABLE]

where the numerical flux $h_{\,j+\frac{1}{2}}$ is given by:

[TABLE]

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 $\mathcal{C}^{\,\star}$, 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:

[TABLE]

where the source term

[TABLE]

is evaluated using Eq. (3.2). For the treatment of the boundary conditions, interested reader is invited to consult [6] 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 [24]. 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.

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:

[TABLE]

The two-step Runge–Kutta(TSRK) methods give the following fully discrete scheme [8]:

[TABLE]

and similarly for $v\,$:

[TABLE]

where $\theta\,$, $\nu_{\,k}\,$, $\mu_{\,k}\,$, $\lambda_{\,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 \Bigl{(}\,\,u^{\,n+1}\,,\,v^{\,n+1}\,\Bigr{)} using the previous values \Bigl{(}\,\,u^{\,n}\,,\,v^{\,n}\,\Bigr{)} and \Bigl{(}\,\,u^{\,n-1}\,,\,v^{\,n-1}\,\Bigr{)} at $t^{\,n}$ and $t^{\,n-1}\,$, respectively.

3.2.2  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:

[TABLE]

where

[TABLE]

Using the Du Fort–Frankel approach [13], the numerical scheme given for Eq. (3.8) is expressed as:

[TABLE]

where the term $2\,w_{\,j}^{\,n}$ has been replaced by $w_{\,j}^{\,n+1}+w_{\,j}^{\,n-1}\,$. The stencil is illustrated in Figure 1(b). Re-arranging Eq. (3.9) to compute $w_{\,j}^{\,n+1}\,$, it gives the following discrete system:

[TABLE]

where

[TABLE]

According to the standard von Neumann stability analysis, the Du Fort–Frankel scheme is unconditionally stable [16, 41]. Further details and examples of applications of this approach are presented in [16, 15].

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 $\mathcal{C}^{\,\star}\,$ is adopted. Thus, the flux $f_{\,j\ +\ \frac{1}{2}}$ is written as:

[TABLE]

where the coefficients are computed as:

[TABLE]

with

[TABLE]

Then, the numerical flux $f_{\,j\pm\frac{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\pm\frac{1}{2}}$ and $h_{\,j\pm\frac{1}{2}}\,$.

The term $\displaystyle\sum_{i=1}^{3}\;\beta_{\,i}\;j_{\,c\,,\,i}\;\frac{\partial v}{\partial x}$ is taken into account by freezing the flux $j_{\,c\,,\,i}^{\,n}$ at each time step $t^{\,n}\,$:

[TABLE]

Then, the flux $g$ in Equation (3.1b) is written as:

[TABLE]

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 [6]. 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:

[TABLE]

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 [19]:

[TABLE]

The stability condition (3.11) is nonlinear in $\Delta x\,$. It can be observed that for large space discretisation $\Delta x\,$, the conditions is equivalent to $\Delta t\,\leqslant\ C_{\,1}\,\Delta 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): $\Delta t\,\leqslant\,C_{\,2}\,\Delta 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 $\mathcal{C}^{\,\star}$. 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\,)\,,\forall\ x\ \in\ \bigl{[}\,x_{\,j}\,,\,x_{\,j+1}\,\bigr{]}\,.

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 $\Delta t\ \sim\ \Delta 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^{\,\mathrm{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:

[TABLE]

where $N_{\,t}$ is the number of temporal steps. The global uniform error is measured using the norm of the space

[TABLE]

is given by the maximum value of $\varepsilon_{\,2}\,(\,x\,)\,$:

[TABLE]

where $\tau$ is the time horizon of simulation. Another criterion is the significative correct digits of the solution, computed according to [38]:

[TABLE]

The last criterion is the computational (CPU) time required by the numerical model to compute the solution. It is measured using Matlab™ 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:

[TABLE]

where super script ${}^{\mathrm{\,num.}}$ stands for the output field obtained by the numerical model and ${}^{\mathrm{\,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 $\mathcal{L}_{\,2}$ error is computed according to:

[TABLE]

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 $\varepsilon_{\,2\,,r}\,$. The quantity $q^{\mathrm{\,app}}$ designates the quantity approximated by the hypothesis and $q^{\mathrm{\,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™ open source toolbox Chebfun [12].

Description of the case

The dimensionless problem is considered for the validation:

[TABLE]

with the following boundary conditions at the air–material interface:

[TABLE]

where $\eta$ corresponds to the projection of the surface normal unity vector $\boldsymbol{n}$ on the axis $x\,$. The simulation is performed for x^{\,\star}\ \in\ \bigl{[}\,0\,,\,1\,\bigr{]} and t^{\,\star}\ \in\ \bigl{[}\,0\,,\,24\,\bigr{]}\,. At the initial state, all fields are set to zero. The boundary conditions are defined as:

[TABLE]

[TABLE]

with the following numerical values for the Biot numbers:

[TABLE]

The numerical values of other dimensionless numbers are:

[TABLE]

with the following functions for the parameters $c^{\,\star}\,$, $k^{\,\star}$ and $a^{\,\star}\,$:

[TABLE]

Other parameters are equal to zero and $r^{\,\star}\,$ 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\in\bigl{\{}\,2\,,\,3\,\Bigl{\}}\, for the two–step Runge–Kutta method, which coefficients are given in the following tables [8]:

[TABLE]

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 $3$b stands for $s\ =\ 2$ and NM $3$c for $s\ =\ 3\,$. The following discretisation parameters are set $\Delta x^{\,\star}\ =\ 0.01$ and $\Delta t^{\,\star}\ =\ 2\cdot 10^{\,-3}\,$. The numerical model NM $1$, based on an Euler explicit approach for the temporal discretisation, is also used with $\Delta x^{\,\star}\ =\ 0.01$ and $\Delta t^{\,\star}\ =\ 10^{\,-4}\,$, which is consistent with the CFL condition. It can be noted that for higher time discretisation (with $\Delta t^{\,\star}\ =\ 0.01$), the model NM $1$ does not compute a bounded solution. The reference solution is the pseudo–spectral one obtained with the Matlab™ open source toolbox Chebfun [12].

Results

Figures 2(a), 2(c) and 2(e) show the evolution in time of the fields at x^{\,\star}\ \in\ \bigl{\{}\,0.1\,,\,0.5\,,\,0.9\,\bigr{\}}\,, 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 $3$b, NM $3$c 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 $3$b 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 $\Delta t\ =\ 10^{\,-4}$ (with $\Delta x^{\,\star}\ =\ 0.01$) by twenty times with $\Delta t\ =\ 2\cdot 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 [25] is shown in Figure 3(a). The wall is composed of a $2\ \mathsf{cm}$ outside coating and two $8\ \mathsf{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\ \mathsf{days}$, the inside temperature and relative humidity are set to approximately $24\ \mathsf{{}^{\,\circ}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 $\sigma_{\,T}^{\,\mathrm{meas}}\ =\ 0.3\ \mathsf{{}^{\,\circ}C}$ and $\sigma_{\,\phi}^{\,\mathrm{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\ =\ \bigl{\{}\,4\,,\,12\,\bigr{\}}\ \mathsf{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 $\sigma^{\,\mathrm{pos}}\ =\ 1\ \mathsf{cm}$ for the sensors at x\ =\ \bigl{\{}\,4\,,\,12\,\bigr{\}}\ \mathsf{cm} and $\sigma^{\,\mathrm{pos}}\ =\ 0.5\ \mathsf{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 $\sigma_{\,P}^{\,\mathrm{meas}}\ =\ 0.5\ \mathsf{\%}\,$. Data measurement is taken with an interval of $5\ \mathsf{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 of the wood fiber have been fully determined in [34, 42]. The properties used for this case are synthesized in Table 3 with other physical constants. For the saturation pressure $P_{\,\mathrm{sat}}$, associating the vapor pressure $P_{\,1}$ and the relative humidity $\phi\,$, the following expression are considered:

[TABLE]

The total sequence of experimental data corresponds to $14\ \mathsf{days}$ and the physical domain is defined for x\,\in\,\bigl{[}\,0\,,\,L\,\bigr{]} with $L\ =\ 16\ \mathsf{cm}\,$.

Accordingly, for the three fields, Dirichlet boundary conditions are considered. Measurements provided by sensors at $x\ =\ 0\ \mathsf{cm}$ and $x\ =\ 16\ \mathsf{cm}$ are used to prescribe both temperature and vapor pressure at the boundaries. For the air pressure $P\,$, at $x\ =\ 0\ \mathsf{cm}\,$, the pressure is imposed using measurements from a differential sensor and considering the reference pressure of $1\ \mathsf{bar}\,$. At $x\ =\ 16\ \mathsf{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\ \mathsf{cm}$). This hypothesis is very likely since the air permeability of the coating is $100$ times lower than wood fiber one. The boundary conditions with their measurement uncertainties are shown in Figure 4. The increase of vapor pressure on the inside part of the PASSYS cell can be noticed. The maximal uncertainty due to sensor measurements reaches $100\ \mathsf{Pa}$ for the vapor pressure and $3\ \mathsf{Pa}$ for the air pressure difference.

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\ =\ \bigl{\{}\,0,\,4,\,8,\,12,\,16\,\bigr{\}}\,\mathsf{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 $\varepsilon_{\,2}$ of $5.21\ \mathsf{Pa}$ and $0.29\ \mathsf{{}^{\,\circ}C}$ and a relative error error $\varepsilon_{\,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 [14].

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 $\varepsilon_{\,2\,,r}$ defined in Section 3.5 and its probability density functions computed using the kernel smoothing function for the extrapolation [22].

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}\,\gg\,P_{\,1}\,$, the constant gas law can be approximated by $R_{\,13}\,\approx\,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 $\sim\,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}\,\approx\,r_{\,12}^{\,\circ}\,$. 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 [3], 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}\,\approx\,c_{\,0}\,\rho_{\,0}\ +\ c_{\,2}\,\mathsf{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 $\sim\,30\,\%$ indicating that the approximation may not provide reliable results for this case study.

The significance of the two transient terms \displaystyle\biggl{(}\,-\ r_{\,12}\,c_{\,qv}\,\frac{\partial P_{\,1}}{\partial t}\ +\ r_{\,12}\,c_{\,qs}\,\frac{\partial\sigma}{\partial t}\,\biggr{)} in the right–hand 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 $\sim\,21\,\%$ with an average value of $\sim\,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:

[TABLE]

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 \sum_{i=1}^{3}\,\boldsymbol{\nabla}\,\bigl{(}\,c_{\,i}\,T\,\bigr{)}\boldsymbol{\cdot}\boldsymbol{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:

[TABLE]

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\,\in\,\bigl{[}\,8\,,\,14\,\bigr{]}\ \mathsf{days}.

Results of the comparison and discussion

If the experimental data was already used in previous works [35, 2, 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 $\sigma_{{}_{\,T}}$ and $\sigma_{{}_{\,P_{\,1}}}\,$, respectively. At the measurement point x_{\,0}\ =\ \bigl{\{}\,4\,,\,8\,,\,12\,\bigr{\}}\ \mathsf{cm}, they are computed by:

[TABLE]

where, as mentioned in Section 5.1, $\sigma^{\,\mathrm{pos}}$ is the sensor position uncertainty. The quantities $\sigma_{{}_{\,T}}^{\,\mathrm{meas}}$ and $\sigma_{{}_{\,P_{\,1}}}^{\,\mathrm{meas}}$ are the sensor measurement uncertainties. For the temperature, $\sigma_{{}_{\,T}}^{\,\mathrm{meas}}$ is directly given by the sensor manufacturer. For the vapor pressure, $\sigma_{{}_{\,P_{\,1}}}^{\,\mathrm{meas}}$ is evaluated using Eq. (5.1):

[TABLE]

where $\sigma_{\,\phi}^{\,\mathrm{meas}}$ is directly given by the sensor manufacturer. The quantities \displaystyle\frac{\partial T}{\partial x}\,\biggl{|}_{\,x\ =\ x_{\,0}} and \displaystyle\frac{\partial P_{\,1}}{\partial x}\,\biggl{|}_{\,x\ =\ x_{\,0}} are evaluated using the results of the numerical model at the points of probes x_{\,0}\ =\ \bigl{\{}\,4\,,\,8\,,\,12\,\bigr{\}}\ \mathsf{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 $3$c. The discretisation parameters are set to $\Delta x^{\,\star}\ =\ 4\cdot 10^{\,\,-2}$ and $\Delta t^{\,\star}\ =\ 5\,\cdot 10^{\,-4}\,$, corresponding to $\Delta x\ =\ 6.4\ \mathsf{mm}$ and $\Delta t\ =\ 1.8\ \mathsf{s}$ in physical dimensions. According to the numerical values of the parameters, we have the following stability conditions:

[TABLE]

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 $\sim\,50$ due to the efficiency of the proposed innovative method. It enables to save important computational efforts. The model needs $117\ \mathsf{min}$ to compute the solution with the given discretisation parameters. It corresponds to a ratio CPU time / time horizon of simulation of $8.3\ \mathsf{min./day}\,$ for the proposed numerical model and a ratio of $139.3\ \mathsf{min./day}$ for the standard approach based on Euler explicit scheme with central finite-differences. In other words, the NM $3$c is sixteen faster to predict the physical phenomena during occurring during one day. It is also important to remark that with the discretisation parameter $\Delta x^{\,\star}\ =\ 4\cdot 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\ =\ \bigl{\{}\,4\,,\,8\,,\,12\,\bigr{\}}\ \mathsf{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\,\in\,\bigl{[}\,8\,,\,12\,\bigr{]}\ \mathsf{days}, particularly at $x\ =\ 8\ \mathsf{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\ \mathsf{cm}$ between the two layers of wood fiberboard. In the literature [9, 20], 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, 6] 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\ \mathsf{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\cdot 10^{\,-13}\ \mathsf{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 [14] 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 [14], 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\ \mathsf{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 [19], the CFL stability condition of the scheme $\Delta t\,\simeq\,\Delta 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 $\Delta t\,\sim\,\Delta 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 [9, 20] 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) [21], 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

[TABLE]

[TABLE]

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.