A New Model For Simulating Heat, Air and Moisture Transport in Porous Building Materials

Julien Berger
LOCIE–CNRS, Université Savoie Mont Blanc, France
Denys Dutykh
LAMA–CNRS, Université Savoie Mont Blanc, France
Nathan Mendes
arXiv:1903.05456v1 [physics.app-ph] 7 Mar 2019
Pontifical Catholic University of Paraná, Brazil
Bolatbek Rysbaiuly
IIT University, Almaty, Kazakhstan
A new model for simulating heat,

air and moisture transport in
porous building materials
arXiv.org / hal
Last modified: March 14, 2019
A new model for simulating heat, air and
moisture transport in porous building
materials
Julien Berger∗ , Denys Dutykh, Nathan Mendes, and Bolatbek Rysbaiuly
Abstract. This work presents a detailed mathematical model combined with an innova-
tive 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 im-
pact on the heat and moisture transfer. The achievement of the mathematical model is de-
tailed 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.
Key words and phrases: heat and mass transfer; porous material; Benchmarking with
experimental data; Scharfetter–Gummel numerical scheme; Du Fort–Frankel nu-
merical 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)
Key words and phrases. heat and mass transfer; porous material; Benchmarking with experimental
data; Scharfetter–Gummel numerical scheme; Du Fort–Frankel numerical scheme.
∗
Corresponding author.
A new model for simulating heat, air and moisture transport 3 / 45
Contents
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Formulation of the physical phenomena . . . . . . . . . . . . . . . . . . . . . . 6

2.1 Porous medium basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Moisture mass balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Moist air mass balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Energy balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.5 Volumetric vapor source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.6 Expression of the fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.7 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.8 Boundary, interface and initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3 Elaborating an efficient numerical model . . . . . . . . . . . . . . . . . . . . 13

3.1 Spatial discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Equation (3.1a), moisture mass balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Equation (3.1b), energy balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Equation (3.1c), moist air mass balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Temporal discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Equations (3.1a) and (3.1b), moisture mass and energy balances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Equation (3.1c), moist air mass balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Extension to nonlinear coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.4 Important features of the numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.5 Efficiency and reliability of a numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4 Validation of the numerical model . . . . . . . . . . . . . . . . . . . . . . . . 22

4.1 Description of the case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5 Comparison of the numerical predictions with experimental data . . . . . . 27

5.1 Experimental set–up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.2 Discussion on some hypotheses regarding the mathematical model . . . . . . . . . . . . . . 29
5.3 Results of the comparison and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
1. 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 [10, 26, 43].
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 COMSOLTM 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 trans-
fer. The former references uses an explicit scheme provided by the commercial COMSOLTM
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 COMSOLTM 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
J. Berger, D. Dutykh et al. 6 / 45
standard hourly climate driven boundary conditions. It should be noted that this condi-
tion 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 tran-
sient 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 [18, 19]. In the context of building porous media, it is suc-
cessfully 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 Schar-
fetter–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 math-
ematical 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 compu-
tational 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 ob-
servations. 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.
2. Formulation of the physical phenomena
First, the mathematical model including the governing equations to describe the physical
phenomena is presented.
2.1. 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 .
3
The volumetric concentration w i kg/m of specie(s) i , which varies in space and time,

def mi ( V
is defined by w i := lim . We denot by m i ( V kg the mass of the specie(s) in
V →0 V
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:
def m 12
w 12 := = f φ , (2.1)
V
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 ≪ 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 R1 T
σ = w 12 − P 1 .
ρ2 R1 T − P1 Π
Assuming the perfect gas law, the volumetric vapor and dry air content are given by:
P1 P3
w1 = Π 1 − σ , w3 = Π 1 − σ .
R1 T R3 T
Last, the liquid volumetric mass content is given by:
P1
w 2 = w 12 − w 1 = w 12 − Π 1 − σ . (2.2)
R1 T
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 n o
= − ∇ · j c,i + Ii , i ∈ 1, 2, 3 , (2.3)
∂t

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:
I1 + I2 = 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:
jc = jd + ja.
It is important to mention that the gravity effects are neglected since the transport phe-
nomena will be investigated on a horizontal plane perpendicular to the gravity forces.
2.2. Moisture mass balance
By summing Eq. (2.3) (i ← 1) and Eq. (2.3) (i ← 2), the moisture transfer equation
is obtained:
∂
w1 + w2 = − ∇ · j c,1 + j c,2 .
∂t
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 ∂w 12 1 ∂P 1
= ,
∂t ∂φ P s ∂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 1
= − ∇ · j c,1 + j c,2 .
∂φ P s ∂t
2.3. 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:
∂
w1 + w3 = − ∇ · j c,1 + j c,3 + I1 , (2.4)
∂t
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 :
P
w1 + w3 = Π(1 − σ), (2.5)
R 13 T
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
R 13 = 3 . (2.6)
P1 R3
1 − P3 1 − R1
Assuming the porosity to be invariant, the following equation is obtained by computing

the differential of Eq. (2.5):
Π(1 − σ) ΠP Π(1 − σ)P
∂ w1 + w3 = ∂P + ∂( 1 − σ ) − ∂T . (2.7)
R 13 T R 13 T R 13 T 2
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 − σ ) ∂P Π P ∂σ Π ( 1 − σ ) P ∂T
= − ∇ · j c,1 + j c,3 + I1 + + .
R 13 T ∂t R 13 T ∂t R 13 T 2 ∂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.
2.4. 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:
3 3
∂ X X
h0 ρ0 + hi wi = −∇· jq + hi j c,i , (2.8)
∂t i=1 i=1

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:
3 3
X ∂T X ∂w i
c0 ρ0 + ci wi + hi = −∇ · jq
i=1
∂t i=1
∂t
3
X 3
X
− hi ∇ · j c,i − ∇ h i · j c , i . (2.9)
i=1 i=1
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 3
X ∂w i X
hi = hi ∇ · j c,i + h1 I1 + h2 I2 .
i=1
∂t i=1
Thus, Eq. (2.9) becomes:

3 3
X ∂T X
c0 ρ0 + ci wi = − ∇ · j q − h1 I1 − h2 I2 − ∇ ci T · j c,i .
i=1
∂t i=1
def
We denote by r 12 := h 1 − h 2 J/kg as the latent heat of vaporization. According to
Kelvin’s law, it is defined as [39]:
◦

r 12 = r 12 + c 1 − c 2 · T − T ref − R 1 T ln φ . (2.10)
Finally, the energy conservation equation yields:
3 3
X ∂T X
c0 ρ0 + ci wi = − ∇ · j q − r 12 I 1 − ∇ ci T · j c,i .
i=1
∂t i=1
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.
2.5. 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
I1 = ∇ · jc,1 . (2.11)
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:

∂w 2
I1 = − + ∇ · jc,2 .
∂t
The latter expression assumes that the moisture concentration w 2 is obtained by the sorp-
tion isotherm w 12 ≈ w 2 . This assumption is not considered in this work, a third
expression of the source is proposed:
∂w 1
I1 = + ∇ · jc,1 ,
∂t
where the vapor concentration w 1 is expressed as:

Π 1 − σ
w1 = .
R1 T
The differential of w 1 gives:
Π(1 − σ) Π P1 Π ( 1 − σ ) P1
∂ w1 = ∂P 1 + ∂( 1 − σ ) − ∂T . (2.12)
R1 T R1 T R1 T 2
Then, we obtain the final expression of the volumetric source term I 1 :
Π ( 1 − σ ) ∂P 1 Π P 1 ∂σ Π ( 1 − σ ) P 1 ∂T
I1 = − − + ∇ · jc,1 ,
R1 T ∂t R 1 T ∂t R1 T 2 ∂t
which depends on the vapor pressure, the amount of liquid in the pores, the temperature
and the variations of the vapor fluxes.
2.6. 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:
k 13
v = − ∇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:
w1 + w3
j c,1 + j c,3 = v,
Π 1 − σ
which can be rewritten using Eq. (2.13) for the velocity:

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

w1 P1
j c , 1 = − k 1 ∇P 1 + v = − k 1 ∇P 1 + v,
Π 1 − σ R1 T

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:
P1
j c , 1 + j c , 2 = − k 2 ∇P 2 − k 1 ∇P 1 + v.
R1 T
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:
P1
j c , 1 + j c , 2 = − k m ∇P 1 + v.
R1 T
The heat flux is expressed as:
T
j q = − k q ∇T + w1 c1 + w3 c3 v.
Π 1 − σ
It is driven by conduction and advection of the moist air.

2.7. Governing Equations
For the sake of clarity, we introduce the following advection coefficients:

def w1 c1 + w3 c3 def v
a q := v, a v := ,
Π 1 − σ R1 T
along with the storage coefficients:

def Π ( 1 − σ ) def Π P P1 def Π ( 1 − σ ) P P1
c a := , c as := − , c at := − ,
R 13 T T R 13 R1 T2 R 13 R1
def Π(1 − σ) def Π ( 1 − σ ) def Π P 1
c av := , c qv := , c qs := ,
R1 T R1 T R1 T
3
def ∂w 12 1 def X Π P1 ( 1 − σ )
cm := , c q := c 0 ρ 0 + ci wi − r 12 ,
∂φ P s i=1
R1 T 2
and the permeability coefficients:

def w 1 + w 3 k 13 def
ka := , k v := k 1 .
Π 1 − σ µ
Finally, the set of partial differential governing equations for the physical problem of heat,
air and moisture transfer through porous building material becomes:
∂P 1
cm = ∇ · k m ∇P 1 − a v P 1 , (2.14a)
∂t
∂T
cq = ∇ · k q ∇T − a q T + r 12 ∇ · k v ∇P 1 − a v P 1 (2.14b)
∂t
3
X ∂P 1 ∂σ
− ∇ c i T · j c , i − r 12 c qv + r 12 c qs ,
i=1
∂t ∂t
∂P
ca = ∇· k a ∇P − ∇· k v ∇P 1 − a v P 1 (2.14c)
∂t
∂P 1 ∂T ∂σ
+ c av + c at + c as .
∂t ∂t ∂t
2.8. 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 P1 − P1,∞ · 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 P1 − P1,∞ · n + g m,∞ .
Similarly, for the heat transfer, the flow is proportional to the surface and ambient condi-
tions. 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 [9, 20]:
j c,A ( xi ) = j c,B ( xi ) .
As initial conditions, the fields may depend on the space coordinate x :
T = Ti ( x ) , P = Pi ( x ) , P1 = P1,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.
3. 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:
n o
Fo a ≫ max Fo v , Fo m , Fo q .
Thus, when using an Euler explicit scheme, since we have a system of equations cou-
pled 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 Schar-
fetter–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 ∆x for the space and ∆t for the time. The spatial cell
C = x j − 1 , x j + 1 is illustrated in Figure 1(a). The discrete values of function u ( x , t )
2 2
n def

are denoted by u j := 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 ∂f
= , (3.1a)
∂t ∂x
∂v ∂g
= + S, (3.1b)
∂t ∂x
∂w ∂h
= + Σ, (3.1c)
∂t ∂x
where
∂u
f = k1 − a1 u ,
∂x
∂v ∂u ∂u ∂σ
g = k 22 − a 22 v + k 21 − a 21 u , S = − c 21 + c 2s ,
∂x ∂x ∂t ∂t
∂w ∂u ∂u ∂v ∂σ
h = k 33 − k 31 + a 31 u , Σ = c 31 + c 32 + c 3s .
∂x ∂x ∂t ∂t ∂t
X3

The term ∇ c i T · j c , i in Eq. (2.14b) is first omitted for the description of the
i=1
numerical method. Its inclusion is detailed in Section 3.3.
3.1. 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 1
= f j+ 1 − f j− 1 .
dt ∆x 2 2
The Scharfetter–Gummel scheme was first proposed h in i[36]. 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 [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:
" #
du j d1
= B Θ u j+1 − B Θ + B − Θ u j + B − Θ u j−1 . (3.2)
dt ∆x
As mentioned in [6], one important advantage of the Scharfetter–Gummel scheme is
the possibility of computing the flux f j+ 1 as wells as the exact solution of u :
2
a1 x
1
u(x) = − f 1 + Cj e k1 , ∀ x ∈ C⋆ , (3.3)
a 1 j+ 2
where
a1 xj

def
u j − u j+1 −
C j := e k1
1 − eΘ
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:

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

du j
S j = − c 21 .
dt
By analogy, using the Scharfetter–Gummel approach, the numerical flux g j+ 1 is as-
2
sumed 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 in the dual cell can be computed. It is important
2
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 MapleTM file. It should be noted that the exact solution v ( x ) ∈ C ⋆ 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):

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

2

∂w ∂u
h j+ 1 = k 33 − k 31 + a 31 u j+ 1 . (3.4)
2 ∂x j+ 1 ∂x j+ 1 2
2 2
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 k 33
= w j+1 − 2 w j + w j−1 (3.5)
dt ∆x 2
a 1 x j+ 1 a 1 x j− 1 !
2 2
k 31 a 1 k k
− Cj e 1 − C j−1 e 1 (3.6)
∆x k 1
a 1 x j+ 1 a 1 x j− 1 !
2 2
a 31 k1
− f j+ 1 − f j− 1 + a 31 Cj e − C j−1 e k 1 + Σj ,
a1 2 2
(3.7)
where the source term
du j
Σ j = c 31
dt
is evaluated using Eq. (3.2). For the treatment of the boundary conditions, interested
reader is invited to consult [6] for the details.
3.2. 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
(a) (b)
Figure 1. Stencils of the Scharfetter–Gummel(a) and

Du Fort–Frankel(b) schemes.
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:
∂u ∂f
= F(u), F(u) = ,
∂t ∂x
∂v ∂g
= G( v ) , G( v ) = + S.
∂t ∂x
The two-step Runge–Kutta(TSRK) methods give the following fully discrete scheme [8]:
s h
n+1 n−1 n
X i
U kn−1 U kn

u = θu + 1 − θ u + ∆t νk F + µk F ,
k=1
s h
X i
U kn n−1 n
U pn−1 U pn

= λk u + 1 − λk u + ∆t akp F + bkp F
p=1
and similarly for v :

s h
n+1 n−1 n
X i
ν k G V kn−1 + µ k G V kn

v = θv + 1 − θ v + ∆t ,
k=1
s h
X i
V kn n−1 n
a k p G V pn−1 + b k p G V pn

= λk v + 1 − λk v + ∆t ,
p=1
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
vocabulary in [8], the approach
to the has

n+1 n+1
two–steps since it computes the fields u ,v using the previous values u , v n
n

and u n−1 , v n−1 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:

dw j k 33
= w j+1 − 2 w j + w j−1 + H j , (3.8)
dt ∆x 2
where
a 1 x j+ 1 a 1 x j− 1 !
2 2
k 31 a 1 k1 k
Hj = Cj e − C j−1 e 1
∆x k 1
a 1 x j+ 1 a 1 x j− 1 !
2 2
a 31 k1
− f j+ 1 − f j− 1 + a 31 Cj e − C j−1 e k 1 + Σj .
a1 2 2
Using the Du Fort–Frankel approach [13], the numerical scheme given for Eq. (3.8)
is expressed as:
w jn+1 − w jn−1

k 33 n n−1 n+1
+ w j−1 + H jn ,
n

= 2
w j+1 − w j + wj (3.9)
2 ∆t ∆x
where the term 2 w jn has been replaced by w jn+1 + w jn−1 . The stencil is illustrated in
Figure 1(b). Re-arranging Eq. (3.9) to compute w jn+1 , it gives the following discrete
system:

n+1 1 − δ n−1 δ n n 2 ∆t
wj = wj + u j+1 + w j−1 + H jn ,
1 + δ 1 + δ 1 + δ
where
def 2 ∆t
δ := k 33 .
∆x 2
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 ap-
proach are presented in [15, 16].
3.3. 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 is written as:
2
∂u
f j+ 1 = k 1 , j+ 1 − a 1 , j+ 1 u ,
2 2 ∂x 2
where the coefficients are computed as:

k 1 , j+ 1 = k 1 u j+ 1 , v j+ 1 , w j+ 1 , a 1 , j+ 1 = a 1 u j+ 1 , v j+ 1 , w j+ 1 ,
2 2 2 2 2 2 2 2
with

1 1 1
u j+ 1 = u j+1 + u j , v j+ 1 = v j+1 + v j , w j+ 1 = w j+1 + w j .
2 2 2 2 2 2
Then, the numerical flux f j± 1 is computed using the approach described in Section 3.1.1.
2
By analogy, the approach is extended for the computation of the fluxes g j± 1 and h j± 1 .
2 2
3
X ∂v
The term βi jc,i is taken into account by freezing the flux j cn, i at each time step
i=1
∂x
tn :

∂v ∂v
βi jc,i = β i j cn, i .
∂x t = t n ∂x
Then, the flux g in Equation (3.1b) is written as:
3
∂v X
n ∂u
g = k 22 − a 22 + β i j c , i v + k 21 − a 21 u ,
∂x i=1
∂x
and the whole description of the Scharfetter–Gummel approach for the spatial dis-
cretisation can be applied.
3.4. 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 Schar-
fetter–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. There-
fore, 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 equa-
tions 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:
∆x 2

1
∆t 6 min min . (3.10)
2 1 6 k , l 6 2 1 6 j 6 N k kl , j
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]:
( −1 )
a kl , j a kl , j ∆x
∆t max max k kl , j max tanh 6 ∆x . (3.11)
16k,l62 16j 6N 16j 6N k kl , j 2 k kl , j
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 6 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 6 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 tem-
poral 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 un-
conditionally 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
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 ∆t ∼ ∆x Explicit, unconditionally stable
the temporal discretisation of these two equations. It should be noted that the Schar-
fetter–Gummel scheme cannot be used for the diffusion-type Equation (3.1c).
3.5. 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:
v
u Nt
def t 1 X
u 2
ε 2 ( x ) := u ( x , t n ) − u ref ( x , t n ) ,
Nt n=1
where N t is the number of temporal steps. The global uniform error is measured using the
norm of the space

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

def
ε ∞ := sup
ε2 ( x ) ,
x ∈ 0,L
where τ is the time horizon of simulation. Another criterion is the significative correct
digits of the solution, computed according to [38]:
u ( x , τ ) − u ref ( x , τ )

def
scd ( u ) := − log 10
.
u ref ( x , τ )
∞
The last criterion is the computational (CPU) time required by the numerical model to
compute the solution. It is measured using MatlabTM 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:
def T num ( x 0 , t ) − T obs ( x 0 , t )

ε rT ( x 0 , t ) := ,
max T obs ( x 0 , t ) − min T obs ( x 0 , t )
t t
def P 1num ( x 0 , t ) − P 1obs ( x 0 , t )

ε rP 1 ( x 0 , t ) := ,
P 1obs ( 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:
v
u Nx 2
uX
app ref
u
u q ( x j , t ) − q ( x j , t )
def u j =1
ε 2 ,r q ( t ) := uu Nx
.
X 2
q ref ( x j , t )
u
t
j=1
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.
4. 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 ref-
erence solution is computed using a numerical pseudo–spectral approach obtained with
the MatlabTM open source toolbox Chebfun [12].
4.1. Description of the case
The dimensionless problem is considered for the validation:

⋆ ∂u ∂ ⋆ ∂u ⋆
c m ⋆ = Fo m km − Pe m a m u ,
∂t ∂x ⋆ ∂x ⋆

⋆ ∂v ∂ ⋆ ∂v ⋆ ⋆ ∂ ⋆ ∂u ⋆
[4pt]c q ⋆ = Fo q kq − Pe q a q v + Fo q γ r kv − Pe v a v u
∂t ∂x ⋆ ∂x ⋆ ∂x ⋆ ∂x ⋆
3
∂u ⋆ ∂σ
X ∂v
− Ko qv c ⋆qv
r ⋆
⋆ ⋆
+ Ko qs c qs r ⋆
− Fo q βi jc,i ,
∂t ∂t i=1
∂x

⋆ ∂w ∂ ⋆ ∂w ∂ ⋆ ∂u ⋆
c a ⋆ = Fo a ka − Fo a δ kv − Pe v a v u
∂t ∂x ⋆ ∂x ⋆ ∂x ⋆ ∂x ⋆
∂u ⋆ ∂σ ⋆ ∂v
+ Ko av c ⋆av + Ko as c as + Ko at c at ,
∂t ⋆ ∂t ⋆ ∂t ⋆
with the following boundary conditions at the air–material interface:
⋆ ∂u

⋆ ⋆
km ⋆
− Pe m a m u = Bi m η u ∞ − u + η g m ,∞ ,
∂x

⋆ ∂v ⋆ ∂u

⋆ ⋆ ⋆
kq − Pe a
q q v + γ r k v − Pe a
v v u = Bi q η v ∞ − v
∂x ⋆ ∂x ⋆

+ γ 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:
! 2
2 2π 2π
u ∞ , L = 0.6 sin π t + sin t , u ∞ , R = 0.9 sin t ,
24 6
2 ! 2
2π 2π 2π
v ∞ , L = 1.2 sin t + sin t , v ∞ , R = 0.5 sin t ,
5 24 3
2 2 !
2π 2π 2π
w( 0 , t ) = 1.3 sin t , w( 1 , t ) = 0.7 sin t + sin t ,
6 4 24
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 ,
⋆
cm = 1 + 0.05 · u 2 + 0.03 · v 2 + 0.01 · w 2 , a ⋆m = 1 + 0.03 · u ,
⋆
km = 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
o 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]:
λ k a kp b kp
Θ νk µk , (4.1a)
s
0.308343 0.154172 0.154172

−0.113988 −0.556994 0.00838564 1.43462
, (4.1b)
0 −0.526458 0.0787465 1.47417 −0.0264581
s = 2
0.0213802 0.264446 −0.507512 0.264446 0 0 0

0.0991119 0.240292 −0.631471 0.597963 0.392328 0 0
0.353437 0.0969341 −0.578148 0.98944 0.582927 0.262283 0 . (4.1c)
0 0.197972 −0.543131 0.521515 0.0867149 0.621047 0.115883
s = 3
The proposed numerical model will be compared to the reference solution. The abbre-
viation 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 Eu-
ler explicit approach for the temporal discretisation, is also used with ∆x ⋆ = 0.01 and
Table 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.
Numerical Model (NM) ∆x ∆t CPU time ε ∞ for u ε ∞ for v ε ∞ for w

NM 1 10 −2 10 −4 t0 7.35 · 10 −3 4.75 · 10 −3 2.1 · 10 −3
NM 3b 10 −2 2 · 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 MatlabTM
open source toolbox Chebfun [12].
4.2. 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.
0.8
0.8 0.48
0.46
0.6 0.44
0.42
0.6 0.18 0.19 0.2
0.4
0.4 0.2
0
0.2
-0.2
0
-0.4
0 5 10 15 20 0 0.2 0.4 0.6 0.8 1
(a) (b)
1.5
1
1
0.5
0.5
0 0
0.16
0.14
-0.5 0.12
-0.5
0.28 0.29 0.3 0.31 0.32
0 5 10 15 20 0 0.2 0.4 0.6 0.8 1
(c) (d)
0.5
0.5
0
0
-0.5 -0.26
-0.5
-0.28
0.75 0.76 0.77
0 5 10 15 20 0 0.2 0.4 0.6 0.8 1
(e) (f)
Figure2. (a, c, e) Time evolution and (b, d, f ) profiles at

t ⋆ = 8 , 16 , 24 of the fields u , v and w .
5. Comparison of the numerical predictions with experi-

mental 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.
5.1. Experimental set–up
The experimental data was obtained at the University of Savoie Mont Blanc, labora-
tory LOCIE (Laboratory of Optimisation of the Conception and Engineering of the En-
vironment) 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 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 un-
certainty σ Tmeas = 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 σ Pmeas = 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.

Sorption curve kg/m 3
w 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 k2 = 0

Air permeability m 2 k 13 = 1.1 · 10 −13

Porosity − Π = 0.9

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 6

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
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 sat ,
associating the vapor pressure P 1 and the relative humidity φ , the following expression
are considered:
!α
P1 ◦ T − T A
φ = , P sat T = P sat , (5.1a)
P sat T TB
◦
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. Mea-
surements 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
(a) (b)
Figure 3. (a) Picture of the PASSYS cell and (b) illustration of the wall
with the implementation of the sensors.
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 Pa for the
vapor pressure and 3 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 = 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 [14].
5.2. 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 [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
(a) (b)
(c)
Figure 4. Measured boundary conditions used for the numerical model.
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 [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 ≈ c 0 ρ 0 + c 2 w 12 . It
1300
24
22
1250
20
1200 18
16
1150 14
12
1100 10
8
1050
0 0.05 0.1 0.15 0 0.05 0.1 0.15
(a) (b)
0
-10
-20
-30
-40
-50
-60
-70
0 0.05 0.1 0.15
(c)
Figure 5. Measured and interpolated initial conditions used for the

numerical model.
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.
∂P 1 ∂σ
The significance of the two transient terms − r 12 c qv + r 12 c qs in the right–
∂t ∂t
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 ∼ 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 volu-
metric 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:
∂w 1
Exact expression: I1 = + ∇ · jc,1 ,
∂t
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 Pgradient of enthalpy in the heat Equation (2.14b). As
3
mentioned in Section 3.3, the term 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:
3
X
Exact expression: ja = aq T + ci jc,i T ,
i=1
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.
5.3. Results of the comparison and discussion
If the experimental data was already used in previous works [2, 35, 46], only the un-
certainties of the sensor measurement were presented. Here, we propose to evaluate the
propagation of uncertainties due to sensor measurement and sensor location. The total un-
certainty on measurement of temperature T and vaporpressure 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
2 ∂T
σ T2 ( x 0 , t ) = σ Tmeas + σ pos ,
∂x x = x 0
!2
2 ∂P 1
σ P2 ( x 0 , t ) = σ Pmeas + σ pos ,
1 1 ∂x x = x 0

where, as mentioned in Section 5.1, σ pos is the sensor position uncertainty. The quantities
σ Tmeas and σ Pmeas are the sensor measurement uncertainties. For the temperature, σ Tmeas is
1
directly given by the sensor manufacturer. For the vapor pressure, σ Pmeas is evaluated using
1
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0.0001 0.001 0.01 0.1 0.3 0.001 0.01 0.1 0.3
(a) (b)
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.
0.506 1
0.504
0.5
0.502
0
0.5
-0.5
0.498
0.496 -1
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(a) (b)
Figure 7. Importance of the approximation of the volumetric vapor source I 1 .
Eq. (5.1):

meas P sat T
σ φmeas σ meas ,

σP = P sat T + φ
1 T − TA T

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

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 6 0.6 s ,
Eq. (2.14b) : ∆t 6 1.8 s ,
Eq. (2.14c) : ∆t 6 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 [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 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 [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 tem-
perature 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 .
6. 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
25
2200
24
2000
23
1800
22
1600
21
1400 20
1200 19
1000 18
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(a) x = 4 cm (b) x = 4 cm
2200
22
2000
1800 20
1600
18
1400
1200 16
1000
14
800
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(c) x = 8 cm (d) x = 8 cm
22
2000
20
1800
18
1600 16
1400 14
1200 12
10
1000
8
800
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(e) x = 12 cm (f) x = 12 cm
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.
(a) (b)
Figure 9. Relative error between the experimental observations and the

numerical predictions for vapor pressure (a) and temperature (b).
10 -3
10 -1
cm
0
cm
-10 -1.5
-20
-2
-30
-40
-2.5
v
-50
-60 -3
-70
-80 -3.5
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(a) (b)
Figure 10. Variation of the differential pressure (a) and of the mass
average velocity (b), both obtained with the numerical predictions (no
experiments available).
efficient numerical scheme. In contrast to earlier proposed models in literature, the ap-
proach 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 physi-
cal 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 Sec-
tion 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
10 -9 10 -7
-3 8
-4 6
-5 4
-6 2
-7 0
-8 -2
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(a) (b)
6
0
-0.01 5
-0.02 4
-0.03
3
-0.04
2
-0.05
1
-0.06
-0.07 0
-0.08 -1
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
(c) (d)
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).
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 solu-
tions. 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 sta-
bility 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 uncon-
ditionally 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
1.2 1
1 0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
-0.05 0 0.05 0.1 0.15 0.2 0 1 2 3 4 5
(a) x = 4 cm (b) x = 4 cm
1.2 1
1 0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
-0.05 0 0.05 0.1 0.15 0.2 0 1 2 3 4 5
(c) x = 8 cm (d) x = 8 cm
1.2 1
1 0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
-0.05 0 0.05 0.1 0.15 0.2 0 1 2 3 4 5
(e) x = 12 cm (f) x = 12 cm
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.
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
Latin letters
aq heat advection coefficient [W/(m 2 · K)]
av vapor advection coefficient [s/m]
c specific heat [J/(kg · K)]
c at storage coefficient [kg · Pa/(J · K)]
cm moisture storage capacity [kg/(m3 · Pa)]
cq 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 )]
jq heat flow [W/m 2 ]
ka permeability coefficient [s 2 /m 2 ]
k1 , k2 , kv , km vapor, liquid, moisture permeability [s]
kq thermal conductivity [W/(m · K)]
k 13 intrinsic permeability [m 2 ]
L length [m]
m mass [kg]
P , P1 , P2 , P3 pressure [Pa]
R 1 , R 3 , R 13 specie gaz constant [J/(kg · K)]
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 ]
Greek letters
αm surface vapour transfer coefficient [s/m]
αq surface heat transfer coefficient [W/(m2 · K)]
φ relative humidity [−]
ρ specific mass [kg/m3 ]
Π material porosity [−]
σ saturation rate [−]
µ dynamic viscosity [Pa · s]
Acknowledgments
The authors acknowledge the Junior Chair Research program “Building performance
assessment, evaluation and enhancement” from the University of Savoie Mont Blanc in col-
laboration with The French Atomic and Alternative Energy Center (CEA) and Scientific
and Technical Center for Buildings (CSTB). The authors thanks the grants from the Min-
istry 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.
References
[1] K. Abahri, R. Bennacer, and R. Belarbi. Sensitivity analyses of convective and diffusive driv-
ing potentials on combined heat air and mass transfer in hygroscopic materials. Numerical
Heat Transfer, Part A: Applications, 69(10):1079–1091, may 2016. 5
[2] ANR Project HYGRO-BAT. Vers une méthode de conception HYGRO-thermique des BA-
Timents performants, 2014. 27, 32
[3] C. Belleudy, A. Kayello, M. Woloszyn, and H. Ge. Experimental and numerical investigations
of the effects of air leakage on temperature and moisture fields in porous insulation. Building
and Environment, 94:457–466, dec 2015. 5, 30
[4] C. Belleudy, M. Woloszyn, M. Chhay, and M. Cosnier. A 2D model for coupled heat, air,
and moisture transfer through porous media in contact with air channels. Int. J. Heat Mass
Transfer, 95:453–465, apr 2016. 5
[5] J. Berger, S. Gasparin, D. Dutykh, and N. Mendes. Accurate numerical simulation of mois-
ture front in porous material. Building and Environment, 118:211–224, jun 2017. 6, 14, 15,
35
[6] J. Berger, S. Gasparin, D. Dutykh, and N. Mendes. On the Solution of Coupled Heat and
Moisture Transport in Porous Material. Transport in Porous Media, 121(3):665–702, feb
2018. 6, 11, 14, 15, 16, 19, 20, 35
[7] J. Berger, S. Guernouti, M. Woloszyn, and C. Buhe. Factors governing the development of
moisture disorders for integration into building performance simulation. J. Building Eng.,
3:1–15, sep 2015. 5
[8] J. Chollom and Z. Jackiewicz. Construction of two-step Runge-Kutta methods with large
regions of absolute stability. J. Comp. Appl. Math., 157(1):125–137, aug 2003. 14, 17, 18, 24
[9] V. P. De Freitas, V. Abrantes, and P. Crausse. Moisture migration in building walls - Analysis
of the interface phenomena. Building and Environment, 31(2):99–108, mar 1996. 13, 34, 40
[10] T. Z. Desta, J. Langmans, and S. Roels. Experimental data set for validation of heat, air
and moisture transport models of building envelopes. Building and Environment, 46(5):1038–
1046, may 2011. 5
[11] G. H. dos Santos and N. Mendes. Heat, air and moisture transfer through hollow porous
blocks. Int. J. Heat Mass Transfer, 52(9-10):2390–2398, apr 2009. 5
[12] T. A. Driscoll, N. Hale, and L. N. Trefethen. Chebfun Guide. Pafnuty Publications, Oxford,
2014. 22, 25
[13] E. C. DuFort and S. P. Frankel. Stability Conditions in the Numerical Treatment of Parabolic
Differential Equations. Mathematical Tables and Other Aids to Computation, 7(43):135, jul
1953. 18
[14] S. Gasparin, J. Berger, D. Dutykh, and N. Mendes. An adaptive simulation of nonlinear
heat and moisture transfer as a boundary value problem. International Journal of Thermal
Sciences, 133:120–139, nov 2018. 29, 35
[15] S. Gasparin, J. Berger, D. Dutykh, and N. Mendes. An improved explicit scheme for whole-
building hygrothermal simulation. Building Simulation, 11(3):465–481, jun 2018. 19
[16] S. Gasparin, J. Berger, D. Dutykh, and N. Mendes. Stable explicit schemes for simulation
of nonlinear moisture transfer in porous materials. J. Building Perf. Simul., 11(2):129–144,
2018. 19
[17] L. Gosse. Computing Qualitatively Correct Approximations of Balance Laws: Exponential-
Fit, Well-Balanced and Asymptotic-Preserving, volume 2 of SIMAI Springer Series. Springer
Milan, Milano, 1 edition, 2013. 15
[18] L. Gosse. Viscous equations treated with L-splines and Steklov-Poincare operator in two di-
mensions. In L. Gosse and R. Natalini, editors, Innovative Algorithms and Analysis. Springer
International Publishing, Milano, 2017. 6, 15
[19] L. Gosse. L-Splines and Viscosity Limits for Well-Balanced Schemes Acting on Linear Para-
bolic Equations. Acta Applicandae Mathematicae, 153(1):101–124, feb 2018. 6, 20, 38
[20] A. S. Guimaraes, J. M. P. Q. Delgado, A. C. Azevedo, and V. P. de Freitas. Interface influence
on moisture transport in buildings. Construction and Building Materials, 162:480–488, feb
2018. 13, 34, 40
[21] C.-E. Hagentoft, A. S. Kalagasidis, B. Adl-Zarrabi, S. Roels, J. Carmeliet, H. Hens,
J. Grunewald, M. Funk, R. Becker, D. Shamir, O. Adan, H. Brocken, K. Kumaran, and
R. Djebbar. Assessment Method of Numerical Prediction Models for Combined Heat, Air
and Moisture Transfer in Building Components: Benchmarks for One-dimensional Cases. J.
Building Phys., 27(4):327–352, apr 2004. 40
[22] P. D. Hill. Kernel estimation of a distribution function. Communications in Statistics -
Theory and Methods, 14(3):605–620, jan 1985. 29
[23] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S.

Woodward. SUNDIALS. ACM Transactions on Mathematical Software, 31(3):363–396, sep
2005. 5
[24] Z. Jackiewicz, R. Renaut, and A. Feldstein. Two-Step Runge-Kutta Methods. SIAM J. Num.
Anal., 28(4):1165–1182, aug 1991. 17
[25] S. O. Jensen. The PASSYS Project Phase 1: Subgroup Model Validation and Development
: Final Report 1986-1989. Technical report, Commission of the European Communities:
Directorate-General for Science, Research, and Development, Brussels, Belgium, 1989. 27
[26] T. Kalamees and J. Kurnitski. Moisture Convection Performance of External Walls and
Roofs. Journal of Building Physics, 33(3):225–247, jan 2010. 5
[27] A. C. Khanduri, T. Stathopoulos, and C. Bédard. Wind-induced interference effects on
buildings - a review of the state-of-the-art. Engineering Structures, 20(7):617–630, jul 1998.
6
[28] H. O. Kreiss and G. Scherer. Method of lines for hyperbolic equations. SIAM Journal on
Numerical Analysis, 29:640–646, 1992. 18
[29] H. M. Kuenzel and J. Kiessl. Simultaneous heat and moisture transport in building compo-
nents: one- and two-dimensional calculation using simple parameters. IRB Verlag, Stuttgart,
1995. 10
[30] J. Langmans, A. Nicolai, R. Klein, and S. Roels. A quasi-steady state implementation of
air convection in a transient heat and moisture building component model. Building and
Environment, 58:208–218, dec 2012. 5
[31] A. V. Luikov. Heat and mass transfer in capillary-porous bodies. Pergamon Press, New York,
1966. 5, 7, 10, 32
[32] N. Mendes, M. Chhay, J. Berger, and D. Dutykh. Numerical methods for diffusion phenomena
in building physics. PUCPRess, Curitiba, Parana, 1 edition, 2017. 5
[33] F. Mnasri, K. Abahri, G. El, R. Bennacer, and S. Gabsi. Numerical analysis of heat, air, and
moisture transfers in a wooden building material. Thermal Science, 21(2):785–795, 2017. 5
[34] H. Rafidiarison, R. Rémond, and E. Mougel. Dataset for validating 1-D heat and mass
transfer models within building walls with hygroscopic materials. Building and Environment,
89:356–368, jul 2015. 28
[35] S. Rouchier, T. Busser, M. Pailha, A. Piot, and M. Woloszyn. Hygric characterization of
wood fiber insulation under uncertainty with dynamic measurements and Markov Chain
Monte-Carlo algorithm. Building and Environment, 114:129–139, mar 2017. 32
[36] D. L. Scharfetter and H. K. Gummel. Large-signal analysis of a silicon Read diode oscillator.
IEEE Transactions on Electron Devices, 16(1):64–77, jan 1969. 6, 15
[37] Sensirion. Digital Humidity Sensor SHT7x (RH/T), 2019. 27
[38] G. Söderlind and L. Wang. Evaluating numerical ODE/DAE methods, algorithms and soft-
ware. J. Comp. Appl. Math., 185(2):244–260, jan 2006. 21
[39] L. Soudani, A. Fabbri, J.-C. Morel, M. Woloszyn, P.-A. Chabriac, H. Wong, and A.-C. Grillet.
Assessment of the validity of some common assumptions in hygrothermal modeling of earth
based materials. Energy and Buildings, 116:498–511, mar 2016. 7, 10
[40] F. Tariku, K. Kumaran, and P. Fazio. Transient model for coupled heat, air and moisture
transfer through multilayered porous media. Int. J. Heat Mass Transfer, 53(15-16):3035–
3044, jul 2010. 5, 10
[41] P. J. Taylor. The stability of the Du Fort-Frankel method for the diffusion equation with
boundary conditions involving space derivatives. The Computer Journal, 13(1):92–97, jan
1970. 19
[42] O. Vololonirina, M. Coutand, and B. Perrin. Characterization of hygrothermal properties
of wood-based products - Impact of moisture content and temperature. Construction and
Building Materials, 63:223–233, jul 2014. 28
[43] L. Wang and H. Ge. Effect of air leakage on the hygrothermal performance of highly insulated
wood frame walls: Comparison of air leakage modelling methods. Building and Environment,
123:363–377, oct 2017. 5
[44] S. Whitaker. Flow in porous media I: A theoretical derivation of Darcy’s law. Transport in
Porous Media, 1(1):3–25, 1986. 7
[45] S. Whitaker. Flow in porous media II: The governing equations for immiscible, two-phase
flow. Transport in Porous Media, 1(2):105–125, 1986. 7
[46] M. Woloszyn, N. Le Pierrès, Y. Kedowidé, J. Virgone, A. Trabelsi, Z. Slimani, E. Mougel,
P. Reymond, H. Rafidiarison, P. Perré, F. Pierre, R. Belarbi, N. Issaadi, K. Abahri, T. Bejat,
A. Piot, E. Wurtz, T. Duforestel, M. Colmet-Daage, B. Perrin, M. Coutand, O. Vololoni-
rina, W. W. Jomaa, J.-S. Lauffer, P. Thiriet, R. Diss, N. Rémond, and O. Legrand. Vers
une méthode de conception HYGRO-thermique des BATiments performants : démarche du
projet HYGRO-BAT. In IBPSA France, page 8, 2014. 27, 32
J. Berger: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LOCIE, 73000
Chambéry, France and LOCIE, UMR 5271 CNRS, Université Savoie Mont Blanc, Campus
Scientifique, F-73376 Le Bourget-du-Lac Cedex, France
E-mail address: [email protected]
URL: https://www.researchgate.net/profile/Julien_Berger3/
D. Dutykh: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAMA, 73000
Chambéry, France and LAMA, UMR 5127 CNRS, Université Savoie Mont Blanc, Campus
Scientifique, F-73376 Le Bourget-du-Lac Cedex, France
URL: http://www.denys-dutykh.com/
N. Mendes: Thermal Systems Laboratory, Mechanical Engineering Graduate Pro-

gram, Pontifical Catholic University of Paraná, Rua Imaculada Conceição, 1155, CEP:
80215-901, Curitiba – Paraná, Brazil
URL: https://www.researchgate.net/profile/Nathan_Mendes/
B. Rysbaiuly: International Information Technology University, Almaty, Manas Str.

34/1, Kazakhstan

A New Model For Simulating Heat, Air and Moisture Transport in Porous Building Materials

Uploaded by

Copyright:

Available Formats

A New Model For Simulating Heat, Air and Moisture Transport in Porous Building Materials

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

A New Model For Simulating Heat, Air and Moisture Transport in Porous Building Materials

Uploaded by

Copyright:

Available Formats

Julien Berger

LOCIE–CNRS, Université Savoie Mont Blanc, France

Pontifical Catholic University of Paraná, Brazil

A new model for simulating heat,

MSC: [2010] 35R30 (primary), 35K05, 80A20, 65M32 (secondary)

2 Formulation of the physical phenomena . . . . . . . . . . . . . . . . . . . . . . 6

3 Elaborating an efficient numerical model . . . . . . . . . . . . . . . . . . . . 13

3.2 Temporal discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.3 Extension to nonlinear coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4 Validation of the numerical model . . . . . . . . . . . . . . . . . . . . . . . . 22

5 Comparison of the numerical predictions with experimental data . . . . . . 27

Moisture is a key factor on durability and performance of buildings. An excessive level

2. Formulation of the physical phenomena

2.1. Porous medium basics

2.2. Moisture mass balance

2.3. Moist air mass balance

Assuming the porosity to be invariant, the following equation is obtained by computing

2.4. Energy balance

Since I 3 = 0, by summing Eq. (2.3) (i ← 1), multiplied by h 1 , Eq. (2.3) (i ← 2),

Thus, Eq. (2.9) becomes:

2.5. Volumetric vapor source

In order to write the system of governing differential equations, it is important to express

2.6. Expression of the fluxes

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

The vapor flow is driven by diffusion and advection:

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

2.7. Governing Equations

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

and the permeability coefficients:

2.8. Boundary, interface and initial conditions

3. Elaborating an efficient numerical model

3.1. Spatial discretisation

3.1.1 Equation (3.1a), moisture mass balance

3.1.2 Equation (3.1b), energy balance

The source term S j is evaluated using Eq. (3.2):

3.1.3 Equation (3.1c), moist air mass balance

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

where the numerical flux h j+ 1 is given by:

3.2. Temporal discretisation

Figure 1. Stencils of the Scharfetter–Gummel(a) and

and similarly for v :

where θ , ν k , µ k , λ k , a kp and b kp are numerical coefficients provided in [8] depending on

3.2.2 Equation (3.1c), moist air mass balance

3.3. Extension to nonlinear coefficients

where the coefficients are computed as:

3.4. Important features of the numerical scheme

Equation Eqs. (3.1a) and (3.1b) Eq. (3.1c)

3.5. Efficiency and reliability of a numerical model

is given by the maximum value of ε 2 ( x ) :

vapor pressure are defined as:

def T num ( x 0 , t ) − T obs ( x 0 , t )

def P 1num ( x 0 , t ) − P 1obs ( x 0 , t )

4. Validation of the numerical model

4.1. Description of the case

The dimensionless problem is considered for the validation:

The numerical values of other dimensionless numbers are:

0.308343 0.154172 0.154172

0.0213802 0.264446 −0.507512 0.264446 0 0 0

Numerical Model (NM) ∆x ∆t CPU time ε ∞ for u ε ∞ for v ε ∞ for w

Figure2. (a, c, e) Time evolution and (b, d, f ) profiles at