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

by   Suelen Gasparin, et al.

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 non-linear 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 non-linear transfer, which is the second case study. Results have shown the error is proportional to O(Δ t). A modified Crank-Nicolson scheme has been proposed in order to avoid sub-iterations to treat the non-linearities 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.



page 1


An improved explicit scheme for whole-building hygrothermal simulation

Although implicit methods require extra calculation, they have been larg...

Two-level schemes for the advection equation

The advection equation is the basis for mathematical models of continuum...

Fast and stable schemes for non-linear osmosis filtering

We consider a non-linear variant of the transport-diffusion osmosis mode...

Explicit Stencil Computation Schemes Generated by Poisson's Formula for the 2D Wave Equation

A new approach to building explicit time-marching stencil computation sc...

Monotonicity considerations for stabilized DG cut cell schemes for the unsteady advection equation

For solving unsteady hyperbolic conservation laws on cut cell meshes, th...

Finite element solution of nonlocal Cahn-Hilliard equations with feedback control time step size adaptivity

In this study, we evaluate the performance of feedback control-based tim...

Optimized Runge-Kutta (LDDRK) timestepping schemes for non-constant-amplitude oscillations

Finite differences and Runge-Kutta time stepping schemes used in Computa...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 [Harris2001, Berger2015a]. 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 [BauklimatikDresden2011], MATCH [Rode2003], MOIST [Burch1993], WUFI [IBP2005] and UMIDUS [Mendes1997, Mendes1999] among others. Moisture models have also been implemented in whole-building simulation tools and tested in the frame of the International Energy Agency Annex , which reported on most of detailed models and their successful applications for accurate assessment of hygrothermal transfer in buildings [Woloszyn2008].

Moisture transfer is represented by a diffusion equation, formulated as:

associated to boundary and initial conditions, where is the diffusion of the material and where is the moisture potential being diffused in the spatial domain during the time interval .We denote and the spatial and time discretisation within those the domains. Due to the nonlinearities of the material properties and due to the non-periodicity 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 [Mendes2005, IBP2005, BauklimatikDresden2011, Steeman2009, Rouchier2013, Janssen2014, Janssen2007]. However, it has the order of accuracy , while the Crank–Nicolson scheme can be used to increase the accuracy to . For the same time and space discretisation, numerical results obtained with this scheme are more accurate than those obtained from Euler implicit scheme. The CrankNicolson scheme has been implemented for instance in [VanGenuchten1982]. 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 [Mendes2005], 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 [Janssen2014], 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 [Tariku2010, Kalagasidis2007]. Nevertheless, this scheme is conditionally stable under the Courant–Friedrichs–Lewy (CFL) condition:

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 DufortFrankel 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 . The moisture transfer occurs according to the liquid and vapour diffusion. The physical problem can be formulated as [Janssen2014]:


where is the volumetric moisture content of the material and and , the vapour and liquid permeabilities.

Eq. (2.1) can be written using the vapour pressure as the driving potential. For this, we consider the physical relation, known as the Kelvin equation, between and :

Thus we have:

As we consider the mass transfer under isothermal conditions, the second term vanishes and we obtain:

In addition, we have:

Considering the relation , obtained from material properties and from the relation between the vapour pressure and the relative humidity , we get:

Eq. (2.1) can be therefore rewritten as:


The material properties , and depend on the vapour pressure . At the material bounding surfaces, Robin-type boundary conditions are considered:


where and are the vapour pressure of the ambient air, and are the liquid flow (driving rain) at the two bounding surfaces. We consider a uniform vapour pressure distribution as initial condition:


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


. Finally, the floating point arithmetics is designed such as the rounding errors are minimal if you manipulate the numbers of the same magnitude

[Kahan1979]. Moreover, the floating point numbers have the highest density in the interval 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 to avoid severe round-off errors and to likely improve the conditioning of the problem in hands. Therefore, we denote as a global moisture transport coefficient, as the moisture storage coefficient and define the following dimensionless quantities:

In this way, the dimensionless governing equations are then written as:


3. Numerical schemes

In order to describe numerical schemes, let us consider a uniform discretisation of the interval :

The time layers are uniformly spaced as well , , The values of function in discrete nodes will be denoted by .

For the sake of simplicity and without loosing the generality, the numerical schemes are explained considering and as constant, noting and the linear diffusion equation:

(a) Euler Explicit
(b) Crank–Nicolson
(c) Dufort–Frankel
(d) Hyperbolisation
Figure 1. Stencils of the numerical schemes.

3.1. The CrankNicolson 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:


This scheme is 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 .

3.2. The Euler explicit scheme

The standard explicit scheme can be written as:


The stencil of this scheme is depicted in Figure 1(a). This discretisation is completed using the two boundary conditions:

where functions 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

By solving Eq. (3.3) with respect to , we obtain a discrete dynamical system

whose starting value is directly obtained from the initial condition:

It is well-known that scheme (3.3) approximates the continuous operator to order . The explicit scheme is conditionally stable under the following CFL-type condition:


Unfortunately, this condition is too restrictive for sufficiently fine discretisations.

3.3. Improved explicit scheme: DufortFrankel method

Using the so-called DufortFrankel method, the numerical scheme is expressed as:


where the term is replaced by . 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 to give the following discrete dynamical system:


The standard von Neumann stability analysis (detailed in the Appendix) shows that the DufortFrankel scheme is unconditionally stable.

The consistency error analysis of the scheme (3.5) shows the following result:



So, from the asymptotic expansion for we obtain that the DufortFrankel scheme is second order accurate in time and:

  • First order accurate in space if

  • Second order accurate in space if

It is important to note that the boundary condition have to be discretised to the second order of accuracy to maintain the method features uniformly.

3.4. Hyperbolisation scheme

From Eq. (3.6), it can be noted that the DufortFrankel 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 dimension Eq. (3.1) perturbed by adding a low magnitude term containing the second derivative in time:


This is the hyperbolic diffusion equation already familiar to us since it appeared in the consistency analysis of the DufortFrankel scheme. Here we perform a singular perturbation by assuming that

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 .

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

3.4.1. Dispersion relation analysis

The classical dispersion relation analysis looks at plane wave solutions:


By substituting this solution ansatz into Eq. (3.1) we obtain the following relation between wave frequency and wavenumber :


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 and non-dispersive otherwise) while the imaginary part describes how different modes dissipate (if ) or grow (if ). 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 :

By solving this quadratic equation with complex coefficients for , we obtain two branches:

This dispersion relation will be analysed asymptotically with being the small parameter. The branch is not of much interest to us since it is constantly damped, i.e.

It is much more instructive to look at the positive branch :

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).

3.4.2. Error estimate

It is legitimate to ask the question how far the solutions of the hyperbolic equation (3.7) are from the solutions of the parabolic diffusion equation (3.1) (for the same initial condition). This question for the initial value problem was studied in [Myshetskaya2015] and we shall provide here only the obtained error estimate. Let us introduce the difference between two solutions:

Then, the following estimate holds

where is the time horizon and

and the domain is defined as

3.4.3. Discretisation

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


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 .

The stability of the scheme (3.10) was studied in [Chetverushkin2012] and the following stability condition was obtained:


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 error between the solution and a reference solution :


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 . The moisture properties are and . The initial vapour pressure in the material is considered uniform , corresponding to a relative humidity of %. The reference time is , 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 and 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 and , respecting the CFL condition . For the hyperbolisation scheme, . The physical phenomena are thus well represented, as illustrated in Figure 3(a) with the time evolution of the vapour pressure at . 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 and , corresponding to the highest and lowest vapour pressure values. All the numerical methods give accurate results as illustrated with the error calculated as a function of in Figure 4. The error for the hyperbolisation method is lower than the others. Indeed the hyperbolisation numerical scheme is of the order 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 . The spatial discretisation is maintained to and for the hyperbolisation scheme. The reference solution has been computed using the Matlab open source package Chebfun [Driscoll2014]. Using the function , it enables to compute a numerical solution of a partial derivative equation using the Chebyshev functions. Results of the error are shown in Figure 5(a). Before the CFL limit, the errors of the Euler, DufortFrankel and hyperbolisation schemes are proportional to . 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 DufortFrankel scheme computes the solution. An interesting point is that the error is proportional to . The error of the DufortFrankel scheme computed solution becomes too high for . For this case, , therefore, the DufortFrankel scheme is first order accurate in space , 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 , it corresponds to . This limit 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 , for different values of parameter . The error reaches a limit lower than the parameter . It can be verified that for the choice , the stability condition corresponds to . For this case, the choice permits to compute the solution with the best accuracy.

Figure 2. Boundary conditions.
Figure 3. Vapour pressure time evolution at (a) and profiles for (b).
Figure 4. error for fixed .
Figure 5. error as a function of for the Euler, DufortFrankel and hyperbolisation schemes (, ) (a) and for the hyperbolisation scheme ().

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 and . This case will be investigated via the DufortFrankel and the improved CrankNicolson schemes. The hyperbolisation approach is not considered, as a stability condition was observed previously in the linear case. First, the DufortFrankel and CrankNicolson schemes are detailed for the nonlinear case. For this, Eq. (2.6) is re-called with a simplified notation:


5.1. The modified CrankNicolson scheme

The straightforward application of the CrankNicolson scheme to Eq. (5.1) yields the following scheme:



However, this approach leads to deal with nonlinearities due to the evaluation of quantities (as ) at the upcoming time layer . To deal with this issue, linearisation techniques as Picard or NewtonRaphson ones [Raphson1690, Cajori1911], requiring a high number of sub-iterations. To overcome these difficulties, it is possible to evaluate the diffusion coefficient at the actual time layer instead of the upcoming [Ascher1995]. Thus, the diffusion flux at the interface becomes:

Finally, the modified CrankNicolson schemes yields to:

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

5.2. The DufortFrankel scheme

In the nonlinear case, the DufortFrankel numerical schemes is written as:


The right-hand side term can be expressed as:


Using the DufortFrankel stencil (see Figure 1(c)), the term is replaced by . Thus, considering Eq. (5.3), the DufortFrankel schemes can be expressed as an explicit scheme: