Suelen Gasparin
Pontifical Catholic University of Paraná, Brazil

Julien Berger
Pontifical Catholic University of Paraná, Brazil

Denys Dutykh

arXiv:1701.07059v3 [] 28 Mar 2017

CNRS, Université Savoie Mont Blanc, France

Nathan Mendes
Pontifical Catholic University of Paraná, Brazil

Stable explicit schemes for
simulation of nonlinear moisture
transfer in porous materials

arXiv.org / hal

Last modified: March 29, 2017

Stable explicit schemes for simulation of
nonlinear moisture transfer in porous
materials
Suelen Gasparin∗, Julien Berger, Denys Dutykh, and Nathan Mendes

Abstract. Implicit schemes have been extensively used in building physics to compute
the solution of moisture diffusion problems in porous materials for improving stability
conditions. Nevertheless, these schemes require important sub-iterations when treating
nonlinear problems. To overcome this disadvantage, this paper explores the use of improved explicit schemes, such as Dufort–Frankel, Crank–Nicolson and hyperbolisation approaches. A first case study has been considered with the hypothesis of linear
transfer. The Dufort–Frankel, Crank–Nicolson and hyperbolisation schemes were
compared to the classical Euler explicit scheme and to a reference solution. Results have
shown that the hyperbolisation scheme has a stability condition higher than the standard Courant-Friedrichs-Lewy (CFL) condition. The error of this schemes depends
on the parameter τ representing the hyperbolicity magnitude added into the equation.
The Dufort–Frankel scheme has the advantages of being unconditionally stable and
is preferable for nonlinear transfer, which is the three others cases studies. Results have
shown the error is proportional to O(∆t) . A modified Crank–Nicolson scheme has
been also studied in order to avoid sub-iterations to treat the nonlinearities at each time
step. The main advantages of the Dufort–Frankel scheme are (i) to be twice faster
than the Crank–Nicolson approach; (ii) to compute explicitly the solution at each time
step; (iii) to be unconditionally stable and (iv) easier to parallelise on high-performance
computer systems. Although the approach is unconditionally stable, the choice of the
time discretisation ∆t remains an important issue to accurately represent the physical
phenomena.
Key words and phrases: moisture diffusion; numerical methods; finite differences; explicit schemes; CFL condition; Dufort–Frankel 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. moisture diffusion; numerical methods; finite differences; explicit schemes;
CFL condition; Dufort–Frankel scheme.
∗
Corresponding author.

Contents
1

Introduction

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2

Moisture transfer in porous materials . . . . . . . . . . . . . . . . . . . . . . . 6

3

Numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.1 The Crank–Nicolson scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 The Euler explicit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 Improved explicit scheme: Dufort–Frankel method . . . . . . . . . . . . . . . . . . . . . . . 10
3.4 Hyperbolisation scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Dispersion relation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Error estimate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.5 Validation of the numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4

Numerical application: linear transfer . . . . . . . . . . . . . . . . . . . . . . 14

5

Extension for nonlinear transfer . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5.1 The modified Crank–Nicolson scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
5.2 The Dufort–Frankel scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5.3 Numerical application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5.4 Further nonlinear case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Driving rain case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
Capillary adsorption case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

6

Conclusion

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

A

Stability analysis of the Dufort–Frankel scheme . . . . . . . . . . . . . . . . 30

B

Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

B.1 Case from Section 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
B.2 Case from Section 5.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
B.3 Case from Section 5.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

References

4 / 34

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5 / 34

Stable explicit schemes

1. Introduction
Excessive levels of moisture may lead to mould growth and may affect the indoor air
quality, the thermal comfort of the occupants and HVAC energy consumption and demand. Moreover, it can deteriorate building façades and decrease envelope durability
[3, 11]. Assessment of relative humidity in constructions is also important for management
and performance of HVAC systems.
Models for moisture transfer in porous building materials have been implemented as
building simulation tools since the nineties in software such as Delphin [2], MATCH
[25], MOIST [4], WUFI [9] and UMIDUS [17, 19] among others. Moisture models have
also been implemented in whole-building simulation tools and tested in the frame of the
International Energy Agency Annex 41, which reported on most of detailed models and
their successful applications for accurate assessment of hygrothermal transfer in buildings
[32].
Moisture transfer is represented by a diffusion equation, formulated as:
∂u
= ∇ · ( ν ∇u ) ,
∂t
associated to boundary and initial conditions, where ν is the diffusion of the material and
where u (x, t) is the moisture potential being diffused in the spatial domain Ω x during the
time interval Ω t .We denote ∆x and ∆t the spatial and time discretisation within those
the domains. Due to the nonlinearities of the material properties and due to the nonperiodicity of the boundary conditions, the models use numerical techniques to compute
the moisture content from the partial differential governing equation.
Due to its property of unconditional stability, the Euler implicit scheme has been used
in many works reported in the literature [2, 9, 12, 13, 18, 26, 28]. However, it has the order
of accuracy O (∆t + ∆x 2 ) , while the Crank–Nicolson scheme can be used to increase
the accuracy to O (∆t 2 + ∆x 2 ) . For the same time and space discretisation, numerical
results obtained with this scheme are more accurate than those obtained from Euler implicit scheme. The Crank–Nicolson scheme has been implemented for instance in [31].
Nevertheless, at every time step, one has to use a tridiagonal solver to invert the linear system of equations to determine the solution value at the following time layer. For instance
in [18], a multi-tridiagonal matrix algorithm has been developed to compute the solution of
coupled equations of nonlinear heat and moisture transfer, using an Euler implicit scheme.
Furthermore, when dealing with nonlinearities of the material properties for instance, one
has to perform sub-iterations to linearise the system, increasing thus the total number of
iterations. In [12], thousands of iterations are required to converge to the solution of a
mass diffusion problem.
On the other hand, an explicit scheme enables a direct computation of the the solution
at the following time layer. Some examples of works based on explicit schemes can be
found in the literature as [16, 29]. Nevertheless, this scheme is conditionally stable under

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

6 / 34

the Courant–Friedrichs–Lewy (CFL) condition:
1
∆t 6
∆x 2 .
2ν
The CFL condition is restrictive for fine discretisations, explaining why few works have
used this approach in building physics.
This paper is devoted to explore the use of improved explicit schemes to overcome the
instability limitation of the standard explicit scheme. The proposed schemes are evaluated
to solve nonlinear transfer of moisture in porous material. The advantages and drawbacks are discussed for two numerical applications. Next Section aims at describing the
physical phenomena of moisture transfer in porous material. In Section 3, basics of the
Dufort–Frankel and the hyperbolisation explicit schemes, are detailed. Then, linear
and nonlinear moisture transfer cases are considered to verify the features of the proposed
approaches.

2. Moisture transfer in porous materials
The physical problem involves one-dimension moisture diffusion through a porous material defined by the spatial domain Ω x = [ 0, L ]. The moisture transfer occurs according
to the liquid and vapour diffusion. The physical problem can be formulated as [12]:


∂P c
∂
∂P v
∂ρ l+v
kl
,
(2.1)
=
+ kv
∂t
∂x
∂x
∂x

where ρ l+v is the volumetric moisture content of the material and k v and k l , the vapour
and liquid permeabilities.
Eq. (2.1) can be written using the vapour pressure P v as the driving potential. For this,
we consider the physical relation, known as the Kelvin equation, between P v and P c :


Pv
P c = R v · T · ln
P s (T )
Rv T
∂P c
=
.
∂P v
Pv
Thus we have:
∂P c ∂P v
∂P c ∂T
∂P c
=
·
+
·
.
∂x
∂P v ∂x
∂T ∂x
As we consider the mass transfer under isothermal conditions, the second term vanishes
and we obtain:
∂P c
R v T ∂P v
=
·
.
∂x
Pv
∂x
In addition, we have:
∂ρ l+v ∂φ ∂P v
∂ρ l+v ∂T
∂ρ l+v ∂φ ∂P v
∂ρ l+v
=
·
·
+
·
≃
·
·
.
∂t
∂φ
∂P v ∂t
∂T
∂t
∂φ
∂P v ∂t

7 / 34

Stable explicit schemes

Considering the relation ρ l+v = f(φ) = f(P v , T ) , obtained from material properties and
from the relation between the vapour pressure P v and the relative humidity φ , we get:
∂ρ l+v
1 ∂P v
= f ′ (P v )
.
∂t
P s ∂t
Eq. (2.1) can be therefore rewritten as:

 ∂P 
1 ∂P v
∂  Rv T
v
f (P v )
kl
.
=
+ kv
P s ∂t
∂x
Pv
∂x
′

(2.2)

The material properties f ′ (P v ), k l and k v depend on the vapour pressure P v . At the
material bounding surfaces, Robin-type boundary conditions are considered:



Rv T
kl
+ kv
Pv


Rv T
+ kv
− kl
Pv


∂P v
= h vL · P v − P vL − g lL ,
∂x

∂P v
= h vR · P v − P vR − g lR ,
∂x

x = 0,

(2.3)

x = L,

(2.4)

where P vL and P vR are the vapour pressure of the ambient air, g lL and g lR are the liquid
flow (driving rain) at the two bounding surfaces. We consider a uniform vapour pressure
distribution as initial condition:
P v = P vi ,

t = 0.

(2.5)

While performing a mathematical and numerical analysis of a given practical problem, it
is of capital importance to obtain a unitless formulation of governing equations, due to a
number of good reasons. First of all, it enables to determine important scaling parameters
(Biot numbers for instance). Henceforth, solving one dimensionless problem is equivalent
to solve a whole class of dimensional problems sharing the same scaling parameters. Then,
dimensionless equations allow to estimate the relative magnitude of various terms, and thus,
eventually to simplify the problem using asymptotic methods [21]. Finally, the floating
point arithmetics is designed such as the rounding errors are minimal if you manipulate
the numbers of the same magnitude [15]. Moreover, the floating point numbers have the
highest density in the interval (0, 1) and their density decays exponentially when we move
further away from zero. So, it is always better to manipulate numerically the quantities at
the order of O(1) to avoid severe round-off errors and to likely improve the conditioning of
Rv T
+ k v as a global moisture
the problem in hands. Therefore, we denote d m = k l ·
Pv
1
as the moisture storage coefficient and define the
transport coefficient, c m = f ′ (P v )
Ps

8 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

following dimensionless quantities:
u =

Pv
,
P vi

uR =

P vR
,
P vi

uL =

P vL
,
P vi

x⋆ =

t⋆ =

t
,
t0

⋆
cm
=

c m · L2
,
d 0m · t 0

⋆
dm
=

dm
,
d 0m

Bi Lv =

Bi R
v =

hR
v ·L
,
d 0m

g l,⋆ L =

g lL · L
,
0 ·Pi
dm
v

g l,⋆ R =

x
,
L
h Lv · L
,
d 0m

g lR · L
.
0 ·Pi
dm
v

In this way, the dimensionless governing equations are then written as:
c ⋆m
⋆
dm

−d ⋆m

∂
∂u
=
⋆
∂t
∂x ⋆



d ⋆m

∂u
∂x ⋆



,


∂u
L
L
=
Bi
·
u
−
u
− g l,⋆ L ,
v
⋆
∂x

∂u
R
R
=
Bi
·
u
−
u
− g l,⋆ R ,
v
⋆
∂x
u = 1,

t⋆ > 0 ,

x⋆ ∈

t⋆ > 0 ,

x⋆ = 0 ,

(2.6b)

t⋆ > 0 ,

x⋆ = 1 ,


x ⋆ ∈ 0, 1 .

(2.6c)

t⋆ = 0 ,




0, 1 ,

(2.6a)

(2.6d)

3. Numerical schemes
In order to describe numerical schemes, let us consider a uniform discretisation of the
interval Ω x
Ωh :
Ωh =

N
−1
[

[ x j , x j+1 ] ,

j =0

xj+1 − x j ≡ ∆x ,

∀j ∈


0, 1, . . . , N − 1 .

The time layers are uniformly spaced as well t n = n ∆t , ∆t = const > 0 , n =
def
0, 1, 2, . . . , N t The values of function u(x, t) in discrete nodes will be denoted by u jn := u (x j , t n ) .
For the sake of simplicity and without loosing the generality, the numerical schemes are
d⋆
explained considering d ⋆m and c ⋆m as constant, noting ν = ⋆m and the linear diffusion
cm
equation:
∂u
= ∇ · ( ν ∇u ) .
∂t

(3.1)

9 / 34

Stable explicit schemes

(a) Euler Explicit

(b) Crank–Nicolson

(c) Dufort–Frankel

(d) Hyperbolisation

Figure 1. Stencils of the numerical schemes.

3.1. The Crank–Nicolson scheme
A very useful method was proposed by Crank & Nicolson (CN) and it can be successfully applied to the diffusion equation (3.1) as well:
n+1
n+1
n
u j−1
− 2 u jn + u nj+1
u j−1
− 2 u jn+1 + u j+1
u jn+1 − u jn
= ν
+
ν
,
∆t
2 ∆x 2
2 ∆x 2
j = 1, . . . , N − 1 ,
n ≥ 0 . (3.2)

This scheme is O(∆t2 + ∆x 2 ) accurate and unconditionally stable. That is why numerical
results obtained with the CN scheme will be more accurate than implicit scheme predictions.
The stencil of this scheme is depicted in Figure 1(b). The CN scheme has all advantages
and disadvantages (except for the order of accuracy in time) of the implicit scheme. At
every time step one has to use a tridiagonal solver to invert the linear system of equations
to determine solution value at the following time layer t = t n+1 .

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

10 / 34

3.2. The Euler explicit scheme
The standard explicit scheme can be written as:
n
u jn+1 − u nj
u j−1
− 2 u jn + u nj+1
= ν
,
j = 1, . . . , N − 1 ,
n > 0 . (3.3)
∆t
∆x 2
The stencil of this scheme is depicted in Figure 1(a). This discretisation is completed using
the two boundary conditions:

u 0n+1 = ψ L (t n+1 , u n+1
, ...),
1
n+1
n+1
uN
= ψ R (t n+1 , u N
−1 , . . . ) ,

where functions ψ l, r ( • ) may depend on adjacent values of the solution whose number
depends on the approximation order of the scheme (here we use the second order in space).
For instance, for the left boundary conditions, we have

−3 u 0n+1 + 4 u 1n+1 − u 2n+1
d ⋆m
= Bi Lv · u 0n+1 − u L − g l,⋆ L .
2 ∆x
By solving Eq. (3.3) with respect to u jn+1, we obtain a discrete dynamical system

∆t
n
n
n
u
−
2
u
+
u
,
u jn+1 = u nj + ν
j−1
j
j+1
∆x 2
whose starting value is directly obtained from the initial condition:
u j0 = 1 .
It is well-known that scheme (3.3) approximates the continuous operator to order O(∆t +
∆x 2 ) . The explicit scheme is conditionally stable under the following CFL-type condition:
1
∆t 6
∆x 2 .
(3.4)
2ν
Unfortunately, this condition is too restrictive for sufficiently fine discretisations.

3.3. Improved explicit scheme: Dufort–Frankel method
Using the so-called Dufort–Frankel method, the numerical scheme is expressed as:

n
n
u jn+1 − u jn−1
u j−1
− u jn−1 + u jn+1 + u j+1
= ν
,
j = 1, . . . , N − 1 ,
n > 0,
2 ∆t
∆x 2
(3.5)
where the term 2 u jn is replaced by u jn−1 + u jn+1. The scheme (3.5) has the stencil depicted
in Figure 1(c). At first glance, the scheme (3.5) looks like an implicit scheme, however, it
is not truly the case. Eq. (3.5) can be easily solved for u jn+1 to give the following discrete
dynamical system:

1 − λ n−1
λ
n
u jn+1 =
u nj+1 + u j−1
,
uj
+
1 + λ
1 + λ

11 / 34

Stable explicit schemes

where:
∆t
.
∆x 2
The standard von Neumann stability analysis (detailed in the Appendix) shows that the
Dufort–Frankel scheme is unconditionally stable.
The consistency error analysis of the scheme (3.5) shows the following result:
def

λ := 2 ν

L nj = ν

where

3
∆t 2 ∂ 2 u
∂u
∂2u
2 ∂ u
1
∆t
+
−
ν
+
6
∆x 2 ∂t2
∂t
∂x2
∂t3
 ∆t 4 
4
∂3 ∂2u
2 ∂ u
2
1
1
− 12 ν ∆t ∆x 3 2 + O
, (3.6)
− 12 ν ∆x
∂x4
∂x ∂t
∆x 2


n
u nj−1 − u jn−1 + u n+1
+ u j+1
u jn+1 − u jn−1
j
− ν
.
2 ∆t
∆x 2
So, from the asymptotic expansion for L nj we obtain that the Dufort–Frankel scheme
is second order accurate in time and:
• First order accurate in space if ∆t ∝ ∆x 3/2
• Second order accurate in space if ∆t ∝ ∆x 2
It is important to note that the boundary condition have to be discretised to the second
order of accuracy O(∆x 2 ) to maintain the method features uniformly.
def
L nj :=

3.4. Hyperbolisation scheme
From Eq. (3.6), it can be noted that the Dufort–Frankel scheme is unconditionally
consistent with the so-called hyperbolic heat conduction equation. Thus, the scheme is a
hidden way to add a small amount of ’hyperbolicity’ into the model (3.1). In this Section
we shall invert the order of operations: first, we perturb Eq. (3.1) in an ad-hoc way and
only after we discretise it with a suitable method. We consider the 1−dimension Eq. (3.1)
perturbed by adding a low magnitude term containing the second derivative in time:
∂u
∂2u
∂2u
+
=
ν
.
(3.7)
∂t2
∂t
∂x2
This is the hyperbolic diffusion equation already familiar to us since it appeared in
the consistency analysis of the Dufort–Frankel scheme. Here we perform a singular
perturbation by assuming that
τ

∂2u
∂u
k
≪
k
k.
∂t2
∂t
The last condition physically means that the new term has only limited influence on the
solution of Eq. (3.7). Here τ is a small ad-hoc parameter whose value is in general related
to physical and discretisation parameters τ = τ (ν, ∆x, ∆t) .
One can notice Eq. (3.7) requires two initial conditions to obtain a well-posed initial
value problem. However, the parabolic Eq. (3.1) is only first order in time and it only
kτ

12 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

requires the knowledge of the initial temperature field distribution. When we solve the
hyperbolic Eq. (3.7), the missing initial condition is simply chosen to be
∂u
= 0,
∂t

t = 0.

3.4.1 Dispersion relation analysis
The classical dispersion relation analysis looks at plane wave solutions:
u(x, t) = u 0 e i (κ x

− ω t)

.

(3.8)

By substituting this solution ansatz into Eq. (3.1) we obtain the following relation between
wave frequency ω and wavenumber k :
ω(κ) = − i ν κ 2 .

(3.9)

The last relation is called the dispersion relation even if the diffusion Eq. (3.1) is not
dispersive but dissipative. The real part of ω contains information about wave propagation
properties (dispersive if Re ω(κ)
6= const and non-dispersive otherwise) while the imaginary
κ
part describes how different modes κ dissipate (if Im ω(κ) < 0) or grow (if Im ω(κ) > 0 ).
The dispersion relation (3.9) gives the damping rate of different modes.
The same plane wave ansatz (3.8) can be substituted into the hyperbolic heat Eq. (3.7)
as well to give the following implicit relation for the wave frequency ω :
− τ ω2 − i ω + ν κ2 = 0 .
By solving this quadratic equation with complex coefficients for ω , we obtain two branches:
√
− i ± 4 ν κ2 τ − 1
ω ± (κ) =
.
2τ
This dispersion relation will be analysed asymptotically with τ ≪ 1 being the small
parameter. The branch ω − (κ) is not of much interest to us since it is constantly damped,
i.e.
i
+ O(1) .
ω − (κ) = −
τ
It is much more instructive to look at the positive branch ω + (κ) :


ω + (κ) = − i ν κ 2 1 + ν κ 2 τ + 2 ν 2 κ 4 τ 2 + O(τ 3 ) .

The last asymptotic expansion shows that for small values of parameter τ , we obtain a
valid asymptotic approximation of the dispersion relation (3.9) for the diffusion equation
(3.1).

13 / 34

Stable explicit schemes

3.4.2 Error estimate
It is legitimate to ask the question how far the solutions u h (x, t) of the hyperbolic
equation (3.7) are from the solutions u p (x, t) of the parabolic diffusion equation (3.1) (for
the same initial condition). This question for the initial value problem was studied in [20]
and we shall provide here only the obtained error estimate. Let us introduce the difference
between two solutions:
def
δu (x, t) := u h (x, t) − u p (x, t) .
Then, the following estimate holds
√


4

√
2 
2 π2
| δu (x, t) | ≤ τ M 1 + √
8 2τ +
T ,
π
2
where T > 0 is the time horizon and
∂ 2u p
def
(ξ, ζ) ,
M := sup
∂ t2
Ω ξ, ζ
and the domain Ω ξ, ζ is defined as
n
def
Ω ξ, ζ := (ξ, ζ) : 0 ≤ ζ ≤ t ,

t − ζ
t − ζo
x − √
.
≤ ξ ≤ x + √
τ
τ

3.4.3 Discretisation
Eq. (3.7) is discretised on the stencil depicted in Figure 1(d):
def

L nj := τ

n
n
u jn+1 − u jn−1
u jn+1 − 2 u jn + u jn−1
u j+1
− 2 u jn + u j−1
+
−
ν
= 0,
∆t 2
2 ∆t
∆x 2
j = 1, . . . , N − 1 ,
n > 0 , (3.10)

Using the standard Taylor expansions, it can be proven that the scheme is consistent
with hyperbolic heat Eq. (3.7) to the second order in space and in time O(∆t 2 + ∆x 2 ) .
The stability of the scheme (3.10) was studied in [7] and the following stability condition
was obtained:
r
∆t
τ
6
.
(3.11)
∆x
ν
The choice of parameter τ is therefore an important issue and will be discussed in next
Sections.

3.5. Validation of the numerical solution
One possible comparison of the numerical schemes can be done by computing the L ∞
error between the solution u num and a reference solution u ref :
def

ε :=

u ref − u num

∞

(3.12a)

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

14 / 34

The computation of the reference solution is detailed in further Sections.

4. Numerical application: linear transfer
A first case of linear moisture transfer is considered. From a physical point of view, the
numerical values correspond to a material length L = 0.1 m. The moisture properties are
d m = 1.97 · 10−10 s and c m = 7.09 · 10−3 kg/m3 /Pa. The initial vapour pressure in the
material is considered uniform P vi = 1.16 · 103 Pa , corresponding to a relative humidity of
50%. The reference time is t 0 = 1 h, thus the total time of simulation corresponds to 120
hours, or five days. The boundary conditions, represented by the relative humidity φ are
given in Figure 2. The sinusoidal variations oscillates between dry and moist state during
the 120 hours. The convective vapour coefficients are set to 2 · 10−7 s/m and 3 · 10−8 s/m
for the left and right boundary conditions, respectively. As the readers may be interested
in testing the numerical schemes proposed, in the present paper dimensionless values are
provided in Appendix B.
The solution of the problem has been first computed for a discretisation ∆x = 10−2
and ∆t = 10−4 , respecting the CFL condition ∆t 6 4.3 · 10 −4. For the hyperbolisation
scheme, τ = ∆t. The physical phenomena are thus well represented, as illustrated in
Figure 3(a) with the time evolution of the vapour pressure at x = 0. The variations follow
the ones of the left boundary conditions. It can be noted a good agreement between the
four numerical schemes. Furthermore, the vapour pressure profile is shown in Figure 3(b)
for t = 19 h and t = 52 h , corresponding to the highest and lowest vapour pressure
values. All the numerical methods give accurate results as illustrated with the L2 error
calculated as a function of x in Figure 4. The error for the hyperbolisation method is lower
than the others. Indeed the hyperbolisation numerical scheme is of the order O(∆t2 ) and
the numerical solution is therefore more accurate.
A numerical analysis of the behaviour of the four numerical schemes has been carried
out for different values of the temporal discretisation ∆t. The spatial discretisation is
maintained to ∆x = 10−2 and τ = ∆t for the hyperbolisation scheme. The reference
solution has been computed using the Matlab open source package Chebfun [8]. Using the
function pde23t, it enables to compute a numerical solution of a partial derivative equation
using the Chebyshev functions. Results of the L 2 error ε are shown in Figure 5(a).
Before the CFL limit, the errors of the Euler, Dufort–Frankel and hyperbolisation
schemes are proportional to O(∆t). As expected, the Euler scheme enables to compute
the solution as far the CFL condition is respected. Above this limit, the solution diverges.
After the CFL limit, as unconditionally stable, the Dufort–Frankel scheme computes
the solution. An interesting point is that the error is proportional to O(∆t2 ). The error
of the Dufort–Frankel scheme computed solution becomes too high for ∆t > 3 · 10 −2 .
3
For this case, ∆t ∝ ∆x 2 , therefore, the Dufort–Frankel scheme is first order accurate
in space O(∆x), explaining why the error is lower for the hyperbolisation scheme.
As the Euler scheme, the hyperbolisation scheme has a stability condition to respect as
reported in Eq. (3.11). For the case τ = ∆t, it corresponds to ∆t 6 9 · 10 −4 . This limit

15 / 34

Stable explicit schemes
lef t
right

1.2

Φ (-)

1
0.8
0.6
0.4
0.2
0
0

20

40

60
t (h)

80

100

120

Figure 2. Boundary conditions.
3500
C heb.
Euler
Duf ort
H yp.

3000

2000

t = 19h

2500

P v (Pa)

P v (Pa)

1500
2000
1500

1000

C heb.
Euler
Duf ort
H yp.

1000

500

t = 53h

500
0
0

20

40

60
t (h)

(a)

80

100

120

0
0

0.02

0.04
0.06
x (m)

0.08

0.1

(b)

Figure 3. Vapour pressure time evolution at x = 0 m (a) and profiles for
t ∈ {19, 53} h (b).

is higher than the CFL condition for the Euler scheme. The error of the hyperbolisation
scheme varies with the choice of the parameter τ . Figure 5(b) gives the variation of the
error ε, for the hyperbolisation scheme, as a function of ∆t, for different values of parameter
τ . The error ε reaches a limit lower than the parameter τ . It can be verified that for the
3
choice τ = ν ∆x , the stability condition corresponds to ∆x 2 . For this case, the choice
τ 6 ∆t permits to compute the solution with the best accuracy.

16 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes
−2

10

ε (-)

Euler
Duf ort
H yp.

−3

10

−4

10

0

0.02

0.04
0.06
x (m)

0.08

0.1

Figure 4. L 2 error for fixed ∆t = 10 −4 .
−1

10 1
10

-1

Euler
Dufort
Hyp.
Crank
O(∆t2 )
O(∆t)

−2

10

10 -2

ε (-)

ε (-)

10

0

10

CF L

−3

10

10 -3
τ = 5 · 10 −2
τ = 5 · 10 −3
τ = 10 −3
τ = 10 −4
τ = ∆t
O(∆t)

−4

10

10

-4

10 -5
10 -4

C F L = ∆x3/2

10 -3

10 -2
∆t (-)

(a)

10 -1

10 0

−4

10

−3

10
∆t (-)

−2

10

(b)

Figure 5. L 2 error as a function of ∆t for the Euler, Dufort–Frankel and
hyperbolisation schemes (∆x = 10−2 , τ = ∆t) (a) and for the hyperbolisation
scheme (∆x = 10−2 ).

5. Extension for nonlinear transfer
The previous case study investigated the use of three numerical schemes for computing
the solution of a linear problem of moisture diffusion. This second case study considers
now nonlinear diffusion, due to material properties depending on the moisture content

17 / 34

Stable explicit schemes

⋆
⋆
dm
( u ) and c m
( u ). This case will be investigated via the Dufort–Frankel and the
improved Crank–Nicolson schemes. The hyperbolisation approach is not considered,
as a stability condition was observed previously in the linear case. First, the Dufort–
Frankel and Crank–Nicolson schemes are detailed for the nonlinear case. For this,
Eq. (2.6) is re-called with a simplified notation:
"
#
∂u
∂u
∂
c(u)
d(u)
.
(5.1)
=
∂t
∂x
∂x

5.1. The modified Crank–Nicolson scheme
The straightforward application of the Crank–Nicolson scheme to Eq. (5.1) yields
the following scheme:
"
 n+ 1
 n+ 1 #

n+1
n
2
2
u
−
u
∂u
∂u
1
j
j
c nj
d
− d
,
(5.2)
=
∆t
∆x
∂x j+ 1
∂x j− 1
2

2

with


d

∂u
∂x

 n+ 1

2

1
j+
2

=

1
2

"

d

∂u
∂x

 n+1
1
j+
2

+



d

∂u
∂x

n

j+

1
2

#

#
"




1
n+1
.
=
d n+1 u j+1
− u jn+1 + d n 1 u nj+1 − u jn
j+
2 ∆x j+ 21
2
However, this approach leads to deal with nonlinearities due to the evaluation of quantities
(as d n+11 ) at the upcoming time layer t = t n+1 . To deal with this issue, linearisation
j+

2

techniques as Picard or Newton–Raphson ones [5, 23], requiring a high number of subiterations. To overcome these difficulties, it is possible to evaluate the diffusion coefficient
at the actual time layer instead of the upcoming [1]. Thus, the diffusion flux at the interface
becomes:
"
#
 n+ 1





2
1
∂u
n+1
n
=
d n 1 u n+1
d
.
+ d n 1 u j+1
− u jn
j+1 − u j
j+
∂x j+ 1
2 ∆x j+ 2
2
2

Finally, the modified Crank–Nicolson schemes yields to:
"
#

∆t
∆t
∆t
n+1
u jn+1 −
−
dn 1 + dn 1
d n 1 u j+1
d n u n+1
1 +
2
2
2 j− 1 j−1
j−
j+
j+
2 ∆x
2
∆x
2
∆x
2
2
2
2
"

#
∆t
∆t
∆t
n
= 1 −
dn 1 + dn 1
d n 1 u j+1
dn un .
u jn +
+
2
2
2 j− 1 j−1
j+
j+
j−
2 ∆x
2
∆x
2
∆x
2
2
2
2

18 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

The combination of implicit-explicit (IMEX) approaches clearly appear in this formulation.
The major advantage over the classical Crank–Nicolson scheme is to avoid sub-iterations
in the solution procedure, without loosing the accuracy and the stability.

5.2. The Dufort–Frankel scheme
In the nonlinear case, the Dufort–Frankel numerical schemes is written as:
"

n
n #
n+1
n−1
u
−
u
∂u
∂u
1
j
j
c nj
d
− d
.
=
2 ∆t
∆x
∂x j+ 1
∂x j− 1
2

(5.3)

2

The right-hand side term can be expressed as:

n
n !


 

∂u
∂u
1
1
n
n
n
n
n
n
d
d 1 u j+1 + d 1 u j−1 − d 1 + d 1 u jn
− d
=
j+
j−
j+
j−
∆x
∂x j+ 1
∂x j− 1
∆x2
2
2
2
2
2
2
(5.4)
u n+1
+ u jn−1
j
.
2
Thus, considering Eq. (5.3), the Dufort–Frankel schemes can be expressed as an explicit scheme:
Using the Dufort–Frankel stencil (see Figure 1(c)), the term u jn is replaced by

u n+1
=
j

λ2
λ0 − λ3
λ1
n
n
· u j+1
+
· u j−1
+
· u jn−1 ,
λ0 + λ3
λ0 + λ3
λ0 + λ3

with
def

def

λ 2 :=

2 ∆t n
d 1,
∆x2 j+ 2


∆t
def
n
n
:=
d 1 +d 1 .
j−
j+
∆x2
2
2
def

λ 0 := c nj ,

λ 1 :=

2 ∆t n
d 1,
∆x2 j− 2

λ3

When dealing with the nonlinearities of the material properties, an interesting feature of
explicit schemes is that it does not require any sub-iterations (using Newton–Raphson
approach for instance). At the time iteration n, the material properties c j , dj+ 1 , dj− 1 are
2

2

explicitly calculated at t n . It should be noted that the material properties evaluated at
j + 21 is formulated as:


uj + uæ+1
dj+ 1 = d
.
2
2

5.3. Numerical application
From a physical point of view, the storage and diffusion coefficients are given in Figures
6(a) and 6(b). Their variations with the relative humidity are similar to the load bearing
material from [12]. The initial vapour pressure is uniform P vi = 1.16 · 10 3 Pa. No

19 / 34

Stable explicit schemes

moisture flow is taken into account at the boundaries. The ambient vapour pressure at the
boundaries are illustrated in Figure 7. At the left boundary, it has a fast drop until the
saturation state and at the right boundary, it has a sinusoidal variation. The material is
thus excited until the capillary state. The convective vapour transfer coefficients are set
to 2 · 10−7 s/m and 3 · 10−8 s/m for the left and right boundary conditions, respectively.
The final simulation time is also fixed to 120 hours. As in the previous case study, the
dimensionless values can be found in Appendix B.
The solution of the problem has been computed with the following discretisation parameters: ∆t = 10 −4 and ∆x = 10 −2 . For this, the Dufort–Frankel, the modified
Crank–Nicolson and the standard Crank–Nicolson numerical schemes have been
used. A sufficiently converged solution, computed with a Euler explicit scheme, is taken
as reference. For the latter, the tolerance is set to ǫ 6 0.01 · ∆t2 to ensure the convergence
of the sub-iterations. The time variation of the vapour pressure according for the bounding
points is given in Figure 8(a). The vapour pressure in the material is increasing according
to the variation at the left boundary condition. There is a delay between the vapour pressure at the left (x = 0 m) and right (x = 0.1 m) bounding points. This increase can also
be observed on the four profiles of vapour pressure illustrated in Figure 8(b). Furthermore,
a break in the slope of the increase of the vapour pressure can be noticed at t = 12 h , due
to the nonlinear behaviour of the material. The vapour pressure at x = 0.1 m slowly oscillated according to the right boundary condition. All the solutions, computed with each
different numerical schemes, have good agreement to represent the physical phenomena
and their L 2 error with the reference solution is lower than 10 −3 , as reported in Figure 9.
The solution has been computed for different values of ∆t, maintaining ∆x = 10 −2 . For
each value of ∆t, the L2 error has been computed between the numerical and a sufficiently
converged reference solution. Results are given in Figure 10. The equivalent CFL condition
has been computed as:
∆t 6

∆x2
.
 ⋆
2 max dc ⋆m(( uu ))
m

As for the linear case, the Euler explicit scheme enables to compute the solution
while the CFL condition is respected. The Dufort–Frankel and modified Crank–
Nicolson numerical schemes are unconditionally stable. An interesting observation is
that the error of the Dufort–Frankel and modified Crank–Nicolson schemes are
proportional to O(∆t). The modification of the Crank–Nicolson scheme, in order to
avoid the sub-iterations due to nonlinearities, loose the O(∆t2 ) accuracy.
Even with ∆t increasing, the schemes are able to compute a solution. However, the
choice of the time discretisation is an important issue to represent accurately the physical
phenomenon. Figures 11(a) and 11(b) show the vapour pressure evolution computed with
∆t = 10 −1 . The solution lacks of accuracy comparing to the reference solution. For
instance, at t = 40, the solution does not represent accurately the decrease of the vapour
pressure. At x = 0.1 m , ∆t is too large to follow the dynamic of the boundary condition.
The error due to the time discretisation of a sinusoidal boundary condition can be expressed

20 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

as:
ǫ = A ω cos ( ωn∆t ) − A

sin ( ω ( n + 1 ) ∆t ) − sin ( ωn∆t )
,
∆t

where A and ω are the amplitude and the frequency of the signal. The error ǫ is given in
Figure 12 as a function of ∆t for this case study. If an accuracy of ǫ 6 10 −1 is required, a
time discretisation lower than ∆t 6 2 · 10 −2 is needed. Therefore, for an unconditionally
stable scheme, the choice of the time discretisation depends on the variation of the boundary
conditions, as well as the diffusion time in the material, in order to compute a solution
representing accurately the physical phenomenon.
With N, the number points due to spatial discretisation, at each time iteration t n the
numbers of operations for each schemes scales with:
Euler implicit:
Crank–Nicolson:
Euler explicit:

O(N NL · N) ,

O(N NL · 2 · N) ,
O(N) ,

Dufort–Frankel

O(N) ,

modified Crank–Nicolson:

O(2 · N) .

The standard Euler implicit and Crank–Nicolson schemes require N NL sub-iterations
to treat the nonlinearities (using Newton–Raphson approach for instance). The Crank–
Nicolson approach needs twice more operations as it combines explicit and implicit approaches. The number of operations of the explicit schemes, as Euler and Dufort–
Frankel, scales with the direct computation of the solution u n . The modified Crank–
Nicolson does not require sub-iterations as the coefficients are expressed at the current
time layer. Thus, the number of operations required is O(2 · N). Generally O(N NL ) ≫ 1,
as reported in [12] where the number of sub-iterations scales between 10 to 30 . Therefore,
for the same discretisation, the Crank–Nicolson scheme requires much more operations
per time step than Dufort–Frankel or modified Crank–Nicolson. For this numerical application, the CPU time of each approach, using ∆t = 10 −4 , has been evaluated
using Matlab platform on a computer with Intel i7 CPU and 8GB of RAM and reported
in Table 1. It has been preferred to focus on the ratio of computer run time rather
than on absolute values, that are system-dependent. However as expected, the Dufort–
Frankel scheme is twice faster than the modified Crank–Nicolson one. The average
number of sub-iterations is O(N NL ) = 3 and as expected, the CPU time of the standard
Crank–Nicolson scheme is six time longer. In addition, another computational advantage of explicit schemes over the implicit ones is their ease to be parallelised. They allow
to achieve almost perfect scaling on high-performance computer systems [7].
Resuming these results, the main advantages of the Dufort–Frankel schemes are (i)
to avoid sub-iterations to treat the nonlinearities and thus be faster than the Crank–
Nicolson approach, (ii) to compute explicitly the solution at each time step, (iii) to be
unconditionally stable, as well as (iv) the ease to be parallelised.

21 / 34

Stable explicit schemes
0

10

−7

10

−8

dM (s)

cM (kg/m3 /s)

10
−1

10

−9

10

−10

10

−11

−2

10

10

0

0.2

0.4
0.6
φ (-)

0.8

(a)

0

0.2

0.4
φ (-)

0.6

0.8

(b)

Figure 6. Variation of the moisture storage c m (a) and diffusion d m (b) as a
function of the relative humidity φ .

Uniform discretization approaches are commonly used in building simulation codes such
as Delphin [2] or WUFI [9] to name a few. They allow to meet a prescribed tolerance
parameter. In the case of a conditionally stable scheme, the time step limitation is imposed
by a CFL-type condition to be met. In other cases, the constraint is imposed by the sharp
variations of the field, which has been illustrated in Figs. 5(a) with the choice of the time
step according to the variation of the boundary condition. This limitation comes indirectly
from the solution accuracy requirement. Comparing it to a scheme with a uniform grid (in
space), it might allow to reduce the CPU-time of the code by redistributing the nodes. On
the other hand, non-uniform grids might degrade the solution accuracy as illustrated in [14].
Non-uniform grids can be also combined with the Dufort–Frankel numerical scheme.
However, their implementation requires further development in order to not reduce the
order of accuracy, which is certainly an interesting theme for further investigation.

5.4. Further nonlinear case studies
Previous subsection illustrated the relevancy of using the explicit Dufort–Frankel scheme
to compute the solution of moisture transfer through porous material. The purpose is now
to explore the use of this numerical scheme for further case studies, typical cases of moisture transfer in building materials [12]. The length of the material is fixed to L = 0.1
m . The initial vapour pressure is P vi = 1.16 · 10 3 Pa , equivalent to a relative humidity
φ = 0.5. To test the robustness of the scheme, with strong nonlinearities, the properties

22 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes
lef t
right

1.2

Φ (-)

1
0.8
0.6
0.4
0.2
0
0

20

40

60
t (h)

80

100

120

Figure 7. Boundary conditions.
1600
1550

Euler
Duf ort
mod. C rank
C rank

x = 0m

1500

1600
t = 40 h

1550
1500
1450

1400

P v (Pa)

P v (Pa)

1450

1350
1300

1510

t = 50 h

1400
1350
1300

t = 20 h

1500

1250

1250
1490

1200
1150
0

1200

1480

x = 0.1m
20

1470
40

40

45

60
t (h)

80

(a)

1150

50

100

120

0

Euler
Duf ort
mod. C rank
C rank
0.02

t = 10 h

0.04
0.06
x (m)

0.08

0.1

(b)

Figure 8. Vapour pressure time evolution at x ∈ {0, 0.1} m (a) and profiles for
t ∈ {10, 20, 40, 50} h (b).

of a load bearing material are taken from the HAMSTAD benchmark 4 [10] and are considered for both cases. For each one, the solution is compared to a sufficiently converged
solution obtained with an Euler explicit scheme. Dimensionless values for each case can
also be found in Appendix B.

23 / 34

Stable explicit schemes
Duf ort
mod. C rank
C rank
−3

ε (-)

10

−4

10

0

0.02

0.04
0.06
x (m)

0.08

0.1

Figure 9. L 2 error for a fixed ∆t = 6 · 10 −6 .
1

10

0

10

−1

CFL

ε (-)

10

−2

10

−3

10

Duf ort
Euler
C rank
O(∆t)

−4

10

−4

10

−3

10

−2

10

−1

10

0

10

∆t (-)

Figure 10. L 2 error (a) as a function of ∆t .
5.4.1 Driving rain case
The first additional case represents the increase of moisture in the material caused by
driving rain at one of the bounding surfaces. For this, at x = 0 m a moisture flux
g lL = 3.4 kg/m2 /s is imposed and there is no transfer with the ambient air. The relative
humidity of the ambient air varies according to a sinusoidal variation, with an amplitude of

24 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes
1300

1600

Duf ort
mod. C rank
Ref.
1250

1400

P v (Pa)

P v (Pa)

1500

1300
1200

1150

Duf ort
mod. C rank
Ref.

1100
0

1200

20

40

60
t (h)

80

1100
0

100

(a)

5

10
t (h)

15

20

(b)

Figure 11. Comparison of reference solution and the one computed with
Dufort–Frankel and modified Crank–Nicolson schemes for ∆t = 10 −1 at
x = 0 m (a) and x = 0.1 m (b).
1

10

Error
∆t = 0.025
0

10

−1

ǫ (-)

10

−2

10

−3

10

−3

10

−2

10
∆t (-)

−1

10

Figure 12. Error due to time discretisation for the right boundary condition u R
(A = 0.8 and ω = 2π).

0.2, a frequency of 1h and a mean of 0.5, at x = 0.1 m. The convective vapour coefficient
equals 3 · 10−8 s/m and the final simulation time is 30 h.

25 / 34

Stable explicit schemes

5.4.2 Capillary adsorption case
This case simulates a capillary adsorption of the material. Thus, the left bounding
surface of the material is maintained at a saturated state (φ = 1). In such case, the
Robin type boundary condition Eq. (2.3) is modified in order to get a Dirichlet one:
P v = P vL .
Here, P vL = 2.33 · 10 3 Pa. The vapour pressure of the ambient air is maintained constant
at P vR = 1.16 · 10 3 Pa, at x = 0.1 m, with a convective vapour coefficient set to 3 · 10−8
s/m. The solution is computed for a final simulation time of 1 h.
5.4.3 Results and discussion
The discretisations used to compute the solution of both cases are ∆x = 0.01 and
∆t = 10−5 . The evolution of the vapour pressure at x = 0 m and x = 0.1 m is given in
Figure 13(a) for the driving rain benchmark. At x = 0 m, the vapour pressure increases
due to the constant rain flux g lL imposed at the surface. At x = 0.1 m, for t 6 5 h, the
vapour pressure varies according to the sinusoidal fluctuations of the boundary conditions.
Then, as the moisture from the rain flux has diffused through the material, the vapour
pressure starts increasing after t = 5 h. By t = 30 h, the whole material is saturated.
The vapour pressure profiles at different times are illustrated in Figure 13(b) for the
capillary adsorption benchmark. At x = 0 m, the vapour pressure is fixed to the saturation
state. The pressure diffuses in the material along the time until reaching the saturation
state in the whole material at t = 0.3 h. The dimensionless moisture transfer coefficient
⋆
dm
, represented in Figure 14, highlights the nonlinearities of the material properties. The
coefficient has O(103 ) orders of variation during the simulation.
The L 2 error has been computed between the numerical solution obtained with the
Dufort–Frankel scheme and the reference one for both cases and illustrated in Figure 15. The error is of the order O(10−6 ) and proves the high accuracy of the solution
computed with the explicit Dufort–Frankel scheme. Furthermore, the solution has
been computed using the standard Crank–Nicolson scheme, with a tolerance set to
ǫ 6 10 −2 · ∆t2 , and the CPU times of the schemes are given in Table 1. For both cases,
the algorithm using the Crank–Nicolson scheme requires around 5 sub-iterations per
time step to meet the tolerance. The Dufort–Frankel scheme computes the solution
almost ten time faster. These additional results enhance the relevance of using the explicit
Dufort–Frankel scheme to compute the solution of moisture transfer problems.

6. Conclusion
Most of the Numerical methods applied to mathematical models used in building physics
are commonly based on implicit schemes to compute the solution of diffusion problems. The
main advantage is due to the stability conditions for the choice of the time discretisation ∆t.

26 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

2600
2200

Euler
Duf ort

t = 0.3 h

2400
x = 0m

2200
2000

1800

P v (Pa)

P v (Pa)

2000

1600
1400

1800
1600
t = 0.1 h

1400
1200

1200
x = 0.1m

Euler
Duf ort

1000
0

10

20

1000
0

30

t = 0.01 h
0.02

t = 0.05 h

0.04
0.06
x (m)

t (h)

(a)

0.08

0.1

(b)

Figure 13. Vapour pressure time evolution for the driving rain benchmark (a)
and profiles for the capillary adsorption benchmark (b).

t ⋆ = 0.3
3

10

2

d ⋆m (-)

10

1

10

0

10

t ⋆ = 0.01
0

0.2

t ⋆ = 0.05
0.4

0.6

t ⋆ = 0.1
0.8

1

x ⋆ (-)

Figure 14. Profiles of the dimensionless moisture transfer coefficient d ⋆m for the
capillary adsorption benchmark.

However, implicit schemes require important sub-iterations when treating nonlinear problems. This work was therefore devoted to explore some improved schemes and more specifically, the Dufort–Frankel, the Crank–Nicolson and the hyperbolisation schemes.
The first one is first- or second-order accurate in space, depending on the choice of ∆t and
has the advantage of being unconditionally stable. The second one is also unconditionally

27 / 34

Stable explicit schemes
−4

10

Rain benchm.
C apillary adsorption benchm.

−5

ε (-)

10

−6

10

−7

10

0

0.02

0.04
0.06
x (m)

0.08

0.1

Figure 15. L 2 error for the driving rain and the capillary adsorption benchmarks.
Table 1. Computational time required for the numerical schemes.
Hygroscopic adsorption
Numerical Scheme

CPU time (s) Average number of iterations

Dufort–Frankel

87

0

modified Crank–Nicolson

190

0

Crank–Nicolson

550

3

Driving rain
Numerical Scheme

CPU time (s) Average number of iterations

Dufort–Frankel

284

0

Crank–Nicolson

2220

5

Capillary adsorption
Numerical Scheme

CPU time (s) Average number of iterations

Dufort–Frankel

180

0

Crank–Nicolson

1410

5

stable and second-order accurate in space and time. The latter is second-order accurate in
time and in space O(∆t2 ) and consistent with the hyperbolic diffusion equation.
The first case study considered a linear diffusive transfer through a porous material. The
Dufort–Frankel, Crank–Nicolson and hyperbolisation schemes were compared to
the classical Euler explicit scheme and to a reference solution obtained using Chebyshev

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

28 / 34

functions. Results have shown that the hyperbolisation scheme has a stability condition
higher than the standard CFL. The error of this scheme depends on parameter τ representing the amount of hyperbolicity added in the equation. An optimal choice seems to
be τ = ∆t . As expected, the Dufort–Frankel and Crank–Nicolson schemes are
unconditionally stable and enable to compute the solution for any choice of the time discretisation ∆t . In addition, for ∆t 6 10−3 , the error of the Dufort–Frankel scheme is
first-order accurate in time, and for ∆t > 10−3 second-order accurate in time. The first
∆x 2
conclusions revealed that for ∆t 6
, it is preferable to use the hyperbolisation scheme,
2ν
for its accuracy. For larger ∆t values, or for nonlinear cases, the Dufort–Frankel and
Crank–Nicolson schemes are preferable due to their stability.
The second case study focused on nonlinear transfer model, with material properties
strongly dependent on the vapour pressure field. The extension of the Dufort–Frankel and
Crank–Nicolson schemes were given specially to treat the nonlinearities of the problem.
A modified Crank–Nicolson was proposed in order to avoid sub-iterations at each time
step of the algorithm. Both Dufort–Frankel and modified Crank–Nicolson schemes
were used to compute the solution of the problem. Results have shown that the error is
proportional to O(∆t). The modified Crank–Nicolson is twice longer than the Dufort–
Frankel to compute the solution, due to the operations required to compute the implicit
and explicit parts of the scheme. The main advantages of the Dufort–Frankel schemes
is (i) to avoid sub-iterations to treat the nonlinearities, (ii) to compute explicitly the solution at each time step, (iii) the unconditionally stable property, as well as (iv) the ease
to be parallelised. Additional case-studies with stronger nonlinearities and sharper profiles
were analysed, enhancing the advantages of such approach. Attention should be paid for
every scheme unconditionally stable because the choice of the time discretisation ∆t is an
important issue to represent accurately the physical phenomena. As mentioned in [22]:
An inexperienced user often interprets this [the unconditionally stable property] to imply that a physically realistic solution will result no matter how
large is the time step, and such user is, therefore surprised to encounter oscillatory solutions. The ’stability’ in a mathematical sense simply ensures that
these oscillations will eventually die out, but it does not guarantee physically
plausible solutions.
Some examples of unrealistic solutions and some advices on the choice of ∆t, considering
the time variations of the boundary conditions, were provided in this study. Keeping this
in mind, the Dufort–Frankel scheme is a valuable option to compute the solution of
nonlinear problems of moisture diffusion in porous materials. Results are encouraging for
the use of this approach for treating problems of coupled heat and moisture transfer in
two- or three-dimensions and explore other methods such as proposed by Saulyev in [27]
to integrate parabolic equations.

29 / 34

Stable explicit schemes

Acknowledgements
The authors acknowledge the Brazilian Agencies CAPES of the Ministry of Education
and CNPQ of the Ministry of Science, Technology and Innovation, for the financial support.

Nomenclature
cm
g
hv
k
L
Pc
Ps
Pv
Rv
T
φ
ρ

Latin letters
moisture storage capacity
liquid flux
vapour convective transfer coefficient
permeability
length
capillary pressure
saturation pressure
vapour pressure
water gas constant
temperature
Greek letters
relative humidity
specific mass

[kg/m3 /Pa]
[kg/m2 /s]
[s/m]
[s]
[m]
[Pa]
[Pa]
[Pa]
[J/(kg · K)]
[K]
[−]
[kg/m3 ]

30 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

A. Stability analysis of the Dufort–Frankel scheme
The stability analysis is performed for the Initial Value Problem of the linear diffusion
equation:
∂u
∂ 2u
− ν
= 0,
x ∈ R,
∂t
∂x 2
The Dufort–Frankel numerical scheme yields to:

n
u j−1
− u jn−1 + u jn+1 + u nj+1
u jn+1 − u jn−1
= ν
,
2 ∆t
∆x 2

ν > 0.

j ∈ Z,

(A.1)

n > 0.
(A.2)

Using the von Neumann’s stability analysis [6], we seek for discrete plane wave solutions
of the form:
u nj = α u 0 e i j k ∆x ,

(A.3)

where k ∆x is the wave length (k > 1), u 0 is the initial wave amplitude (t = 0) and α
is the amplification factor between two successive time layers.
In general α = α( k ∆x ) and we shall compute this factor for the Dufort–Frankel scheme,
knowing that:
1 n
n
u jn+1 = α u jn ,
u jn−1 =
u ,
u j±1
= e ± i k ∆x u nj .
α j
Thus, Eq. (A.2) becomes:




1
1
α −
= µ 2 cos(σ) − α −
,
(A.4)
α
α
def

def 2 ν ∆t
.
∆x2

where σ := k ∆x and µ :=

α+
α−

Eq. (A.4) admits two roots:
p
1 − µ 2 sin 2 (σ)
µ cos(σ) +
,
=
1 + µ
p
µ cos(σ) −
1 − µ 2 sin 2 (σ)
=
,
1 + µ

In Figures 16(a) and 16(b), we show the dependence of the absolute value of the amplification factors α ± on cos(σ) for two fixed values of the parameter µ — one below
and one above the critical CFL value µ = 1 . It is important to notice that α ± 6 1
for all values of the parameter σ . It can be shown (see [24]) that this property holds for
all values of µ . This observation completes the proof of the unconditional stability of the
DF scheme. The Boundary Value Problem (BVP) for equation Eq. (A.1) was studied in
[30]. Namely, it was shown that the same conclusions on the unconditional stability of the
Dufort–Frankel scheme hold for BVP demonstrated that the boundary conditions are
appropriately discretised.

31 / 34

Stable explicit schemes

1

1
- -α 2 -

0.8

- -α 1 -

- -α 2 -

0.8

- -α-

0.6

- -α-

0.6

- -α 1 -

0.4

0.4

0.2

0.2

0

0

-1

-0.5

0
cos(σ)

0.5

1

-1

-0.5

(a)

0
cos(σ)

0.5

1

(b)

Figure 16. Amplification factors for the Dufort–Frankel scheme for
µ = 0.95 (a) and µ = 1.125 (b).

B. Dimensionless parameters
B.1. Case from Section 4
Problem (2.6) is considered with g l,⋆ L = g l,⋆ R = 0 . The dimensionless properties of
the material are equal to d ⋆m = 1 and c ⋆m = 8.6 . The final simulation time is fixed
to τ ⋆ = 120, while Biot numbers are Bi Lv = 101.5 and Bi R
v = 15.2 . The boundary
conditions are expressed as:




1
1
2π t ⋆
2π t ⋆
L ⋆
u (t ) = 1 +
+
,
sin
sin
2
24
2
4


2π t ⋆
4
R ⋆
.
sin
u (t ) = 1 +
5
12

B.2. Case from Section 5.3
The dimensionless properties of the materials are:
h
2 i
,
d ⋆m = 1 + 0.91 u + 600 · exp −10 u − 1.9
i
h
2
.
c ⋆m = 900 − 656 u + 104 · exp −5 u − 1.9

The Biot numbers are also assumed to be Bi Lv = 101.5 and Bi R
v = 15.2 . The ambient
water vapour pressure at the boundaries are different from the previous case study. At the

32 / 34

S. Gasparin, J. Berger, D. Dutykh & N. Mendes



left boundary, u L has a fast drop until the saturation state u L = 2, ∀t ∈ 10, 40 and at
the right boundary, u R = 1 + 0.8 sin ( 2π t ⋆ ) .

B.3. Case from Section 5.4
For both case-studies, an analytic expression of the material properties has been fitted,
which dimensionless formulation is:
h
2 i
⋆
dm
= 0.85 u −0.71 + 900 exp −8 u − 2
,
h
2 i
⋆
2 −0.53
c m = 1.69 · 10 u
+ 3 exp −9 u − 1.3
.
The dimensionless parameters for the driving rain case (Section 5.4.1) are:

g l,⋆ L = 14.7 , Bi Lv = 0 , Bi R
u R = 1 + 0.4 sin 2 π t ⋆ , τ ⋆ = 30 ,
v = 15.2 ,
while, for the capillary adsorption case (Section 5.4.2), they are:
uL = 2 ,

Bi R
v = 15.2 ,

uR = 1 ,

τ ⋆ = 1.

References
[1] U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton. Implicit-Explicit methods for timedependent partial differential equations. SIAM J. Numer. Anal., 32(3):797–823, 1995. 17
[2] B. Bauklimatik Dresden.
Simulation program for the calculation of coupled
heat, moisture, air, pollutant, and salt transport.
http://www.bauklimatikdresden.de/delphin/index.php?aLa=en, 2011. 5, 21
[3] 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
[4] D. M. Burch. An Analysis of Moisture Accumulation in Walls Subjected to Hot and Humid
Climates. ASHRAE Transactions, 93(16):429–439, 1993. 5
[5] E. Cajori. Historical note on the Newton-Raphson method ofapproximation. Amer. Math.,
18:29–32, 1911. 17
[6] J. G. Charney, R. Fjörtoft, and J. Neumann. Numerical Integration of the Barotropic Vorticity Equation. Tellus, 2(4):237–254, nov 1950. 30
[7] B. N. Chetverushkin and A. V. Gulin. Explicit schemes and numerical simulation using
ultrahigh-performance computer systems. Doklady Mathematics, 86(2):681–683, sep 2012.
13, 20
[8] T. A. Driscoll, N. Hale, and L. N. Trefethen. Chebfun Guide. Pafnuty Publications, Oxford,
2014. 14
[9] I. Fraunhofer. Wufi. http://www.hoki.ibp.fhg.de/wufi/wufi_frame_e.html, 2005. 5, 21
[10] 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

Stable explicit schemes

[11]
[12]
[13]

[14]
[15]
[16]

[17]
[18]

[19]

[20]
[21]
[22]
[23]

[24]
[25]
[26]

[27]
[28]

33 / 34

and Moisture Transfer in Building Components: Benchmarks for One-dimensional Cases. J.
Building Phys., 27(4):327–352, apr 2004. 22
S. Y. Harris. Building Pathology: Deterioration, Diagnostics, and Intervention. Wiley, New
York, 2001. 5
H. Janssen. Simulation efficiency and accuracy of different moisture transfer potentials.
Journal of Building Performance Simulation, 7(5):379–389, sep 2014. 5, 6, 18, 20, 21
H. Janssen, B. Blocken, and J. Carmeliet. Conservative modelling of the moisture and heat
transfer in building components under atmospheric excitation. Int. J. Heat Mass Transfer,
50(5-6):1128–1140, mar 2007. 5
L. Jianchun, G. A. Pope, and K. Sepehrnoori. A high-resolution finite-difference scheme for
nonuniform grids. Appl. Math. Model., 19(3):162–172, mar 1995. 21
W. Kahan and J. Palmer. On a proposed floating-point standard. ACM SIGNUM Newsletter,
14(si-2):13–21, oct 1979. 7
A. S. Kalagasidis, P. Weitzmann, T. R. Nielsen, R. Peuhkuri, C.-E. Hagentoft, and C. Rode.
The International Building Physics Toolbox in Simulink. Energy and Buildings, 39(6):665–
674, jun 2007. 5
N. Mendes. Models for prediction of heat and moisture transfer through porous building
elements. PhD thesis, Federal University of Santa Catarina - UFSC, 1997. 5
N. Mendes and P. C. Philippi. A method for predicting heat and moisture transfer through
multilayered walls based on temperature and moisture content gradients. Int. J. Heat Mass
Transfer, 48(1):37–51, 2005. 5
N. Mendes, I. Ridley, R. Lamberts, P. C. Philippi, and K. Budag. Umidus: A PC program
for the Prediction of Heat and Mass Transfer in Porous Building Elements. In IBPSA 99,
pages 277–283, Japan, 1999. International Conference on Building Performance Simulation.
5
E. E. Myshetskaya and V. F. Tishkin. Estimates of the hyperbolization effect on the heat
equation. Comp. Math. Math. Phys., 55(8):1270–1275, aug 2015. 13
A. Nayfeh. Perturbation Methods. Wiley-VCH, New York, 1 edition, 2000. 7
S. V. Patankar. Numerical Heat Transfer and Fluid Flow. CRC Press, United States of
America, 1980. 28
J. Raphson. Analysis aequationum universalis seu adaequationes algebraicas resolvendas
methodus generalis, et expedita, ex nova infinitarum serierum doctrina deducta ac demonstrata. Microfilm copy: University Microfilms, Ann Arbor(MI), 1690. 17
R. D. Richtmyer and K. W. Morton. Difference Methods for Initial-Value Problems. Interscience Publishers, New York, 1967. 30
C. Rode and K. Grau. Whole Building Hygrothermal Simulation Model. ASHRAE Transactions, 109(1):572–582, 2003. 5
S. Rouchier, M. Woloszyn, G. Foray, and J.-J. Roux. Influence of concrete fracture on the
rain infiltration and thermal performance of building facades. Int. J. Heat Mass Transfer,
61:340–352, jun 2013. 5
V. K. Saulyev. Integration of parabolic equations by the grid method. Fizmatgiz, Moscow,
1960. 28
H.-J. Steeman, M. Van Belleghem, A. Janssens, and M. De Paepe. Coupled simulation
of heat and moisture transport in air and porous materials for the assessment of moisture
related damage. Building and Environment, 44(10):2176–2184, oct 2009. 5

S. Gasparin, J. Berger, D. Dutykh & N. Mendes

34 / 34

[29] 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
[30] 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. 30
[31] M. T. van Genuchten. A comparison of numerical solutions of the one-dimensional unsaturated - saturated flow and mass transport equations. Adv. in Water Ressources, 5(1):47–55,
mar 1982. 5
[32] M. Woloszyn and C. Rode. Tools for performance simulation of heat, air and moisture
conditions of whole buildings. Building Simulation, 1(1):5–24, mar 2008. 5
S. Gasparin: Thermal Systems Laboratory, Mechanical Engineering Graduate Program, Pontifical Catholic University of Paraná, Rua Imaculada Conceição, 1155, CEP:
80215-901, Curitiba – Paraná, Brazil
E-mail address: suelengasparin@hotmail.com
URL: https://www.researchgate.net/profile/Suelen_Gasparin/
J. Berger: Thermal Systems Laboratory, Mechanical Engineering Graduate Program,
Pontifical Catholic University of Paraná, Rua Imaculada Conceição, 1155, CEP: 80215901, Curitiba – Paraná, Brazil
E-mail address: Julien.Berger@pucpr.edu.br
URL: https://www.researchgate.net/profile/Julien_Berger3/
D. Dutykh: LAMA, UMR 5127 CNRS, Université Savoie Mont Blanc, Campus Scientifique, 73376 Le Bourget-du-Lac Cedex, France
E-mail address: Denys.Dutykh@univ-savoie.fr
URL: http://www.denys-dutykh.com/
N. Mendes: Thermal Systems Laboratory, Mechanical Engineering Graduate Program, Pontifical Catholic University of Paraná, Rua Imaculada Conceição, 1155, CEP:
80215-901, Curitiba – Paraná, Brazil
E-mail address: Nathan.Mendes@pucpr.edu.br
URL: https://www.researchgate.net/profile/Nathan_Mendes/

