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

TL;DR

This paper investigates improved explicit numerical schemes like Dufort-Frankel, Crank-Nicolson, and hyperbolisation for simulating nonlinear moisture transfer in porous materials, aiming to enhance stability and computational efficiency over traditional implicit methods.

Contribution

It introduces and compares explicit schemes for nonlinear moisture transfer, highlighting the advantages of Dufort-Frankel in stability, speed, and parallelization, with a modified Crank-Nicolson scheme to handle non-linearities.

Findings

Dufort-Frankel scheme is unconditionally stable and faster than Crank-Nicolson.

Hyperbolisation scheme has a higher stability condition than CFL.

Modified Crank-Nicolson avoids sub-iterations for non-linear problems.

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 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 \tau 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(\Delta 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 $\Delta t$ remains an important issue to accurately represent the physical phenomena.

Tables

Table 1. Table 1. Computational time required for the numerical schemes.

Hygroscopic adsorption
Numerical Scheme	CPU time ($𝗌$)	Average number of iterations
Dufort–Frankel	87	0
modified Crank–Nicolson	190	0
Crank–Nicolson	550	3
Driving rain
Numerical Scheme	CPU time ($𝗌$)	Average number of iterations
Dufort–Frankel	284	0
Crank–Nicolson	2220	5
Capillary adsorption
Numerical Scheme	CPU time ($𝗌$)	Average number of iterations
Dufort–Frankel	180	0
Crank–Nicolson	1410	5

Equations

$$\frac{\partial u}{\partial t}=\nabla \cdot \left(\nu \nabla u\right),$$

$$\Delta t⩽\frac{1}{2\nu }\Delta x^2.$$

$$\frac{\partial {\rho }_{l+v}}{\partial t}=\frac{\partial }{\partial x}\left(k_l\frac{\partial P_c}{\partial x}+k_v\frac{\partial P_v}{\partial x}\right),$$

$$P_c$$

$$\frac{\partial P_c}{\partial P_v}$$

$$\frac{\partial P_c}{\partial x}=\frac{\partial P_c}{\partial P_v}\cdot \frac{\partial P_v}{\partial x}+\frac{\partial P_c}{\partial T}\cdot \frac{\partial T}{\partial x}.$$

$$\frac{\partial P_c}{\partial x}=\frac{R_{v}T}{P_v}\cdot \frac{\partial P_v}{\partial x}.$$

$$\frac{\partial {\rho }_{l+v}}{\partial t}=\frac{\partial {\rho }_{l+v}}{\partial \varphi }\cdot \frac{\partial \varphi }{\partial P_v}\cdot \frac{\partial P_v}{\partial t}+\frac{\partial {\rho }_{l+v}}{\partial T}\cdot \frac{\partial T}{\partial t}\simeq \frac{\partial {\rho }_{l+v}}{\partial \varphi }\cdot \frac{\partial \varphi }{\partial P_v}\cdot \frac{\partial P_v}{\partial t}.$$

$$\frac{\partial {\rho }_{l+v}}{\partial t}={\mathrm{f}}^{\prime }(P_v)\frac{1}{P_s}\frac{\partial P_v}{\partial t}.$$

$$\displaystyle\mathrm{f}^{\,\prime}(P_{\,v})\;\frac{1}{P_{\,s}}\;\frac{\partial P_{\,v}}{\partial t}\ =\ \frac{\partial}{\partial x}\biggl{[}\,\Bigl{(}\,k_{\,l}\,\frac{R_{\,v}\,T}{P_{\,v}}\ +\ k_{\,v}\,\Bigr{)}\,\frac{\partial P_{\,v}}{\partial x}\,\biggr{]}\,.$$

$$\left(k_l\frac{R_{v}T}{P_v}+k_v\right)\frac{\partial P_v}{\partial x}$$

$$-\left(k_l\frac{R_{v}T}{P_v}+k_v\right)\frac{\partial P_v}{\partial x}$$

$$P_v$$

$$u=\frac{P_v}{P_v^i},$$

$$t^{\star }=\frac{t}{t^0},$$

$${\mathrm{Bi}}_v^{\mathrm{R}}=\frac{h_v^{\mathrm{R}}\cdot L}{d_m^0},$$

$$c_m^{\star }\frac{\partial u}{\partial t^{\star }}$$

$$d_m^{\star }\frac{\partial u}{\partial x^{\star }}$$

$$-d_m^{\star }\frac{\partial u}{\partial x^{\star }}$$

$$u$$

$$\Omega_{\,h}\ =\ \bigcup_{j\,=\,0}^{N-1}[\,x_{\,j},\;x_{\,j+1}\,]\,,\qquad x_{j+1}\ -\ x_{\,j}\ \equiv\ \Delta x\,,\quad\forall j\ \in\ \bigl{\{}0,\,1,\,\ldots,\,N-1\bigr{\}}\,.$$

$$\frac{\partial u}{\partial t}=\nabla \cdot \left(\nu \nabla u\right).$$

$$\begin{gathered} \frac{u_j^{n+1}-u_j^n}{\Delta t}=\nu \frac{u_{j-1}^n-2u_j^n+u_{j+1}^n}{2\Delta x^2}+\nu \frac{u_{j-1}^{n+1}-2u_j^{n+1}+u_{j+1}^{n+1}}{2\Delta x^2}, \\ j=1,\dots ,N-1,n\ge 0. \end{gathered}$$

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=\nu \frac{u_{j-1}^n-2u_j^n+u_{j+1}^n}{\Delta x^2},j=1,\dots ,N-1,n⩾0.$$

$$u_0^{n+1}={\psi }_{\mathrm{L}}(t^{n+1},u_1^{n+1},\dots ),$$

$$u_N^{n+1}={\psi }_{\mathrm{R}}(t^{n+1},u_{N-1}^{n+1},\dots ),$$

$$d_m^{\star }\frac{-3u_0^{n+1}+4u_1^{n+1}-u_2^{n+1}}{2\Delta x}={\mathrm{Bi}}_v^{\mathrm{L}}\cdot \left(u_0^{n+1}-u^{\mathrm{L}}\right)-g_{l,\mathrm{L}}^{\star }.$$

$$\displaystyle u_{\,j}^{\,n+1}\ =\ u_{\,j}^{\,n}\ +\ \nu\;\frac{\Delta t}{\Delta x^{\,2}}\;\bigl{(}\,u_{\,j-1}^{\,n}\ -\ 2\,u_{\,j}^{\,n}\ +\ u_{\,j+1}^{\,n}\,\bigr{)}\,,$$

$$u_j^0=1.$$

$$\Delta t⩽\frac{1}{2\nu }\Delta x^2.$$

$$\displaystyle\frac{u_{\,j}^{\,n+1}\ -\ u_{\,j}^{\,n-1}}{2\,\Delta t}\ =\ \nu\;\frac{u_{\,j-1}^{\,n}\ -\ \bigl{(}u_{\,j}^{\,n-1}\ +\ u_{\,j}^{\,n+1}\bigr{)}\ +\ u_{\,j+1}^{\,n}}{\Delta x^{\,2}}\,,\qquad j\ =\ 1,\,\ldots,\,N-1\,,\qquad n\ \geqslant\ 0\,,$$

$$\displaystyle u_{\,j}^{\,n+1}\ =\ \frac{1\ -\ \lambda}{1\ +\ \lambda}\;u_{\,j}^{\,n-1}\ +\ \frac{\lambda}{1\ +\ \lambda}\;\bigl{(}u_{\,j+1}^{\,n}\ +\ u_{\,j-1}^{\,n}\bigr{)}\,,$$

$$\lambda \stackrel{\mathrm{def}}{:=}2\nu \frac{\Delta t}{\Delta x^2}.$$

$$\mathcal{L}_{\,j}^{\,n}\ =\ \nu\;\frac{\Delta t^{\,2}}{\Delta x^{\,2}}\;\frac{\partial^{2}u}{\partial t^{2}}\ +\ \frac{\partial u}{\partial t}\ -\ \nu\,\frac{\partial^{2}u}{\partial x^{2}}\ +\ {\textstyle{1\over 6}}\,\Delta t^{\,2}\,\frac{\partial^{3}u}{\partial t^{3}}\\ \ -\ {\textstyle{1\over 12}}\,\nu\,\Delta x^{\,2}\,\frac{\partial^{4}u}{\partial x^{4}}\ -\ {\textstyle{1\over 12}}\,\nu\,\Delta t^{\,2}\,\Delta x\,\frac{\partial^{3}}{\partial x^{3}}\frac{\partial^{2}u}{\partial t^{2}}\ +\ \mathcal{O}\Bigl{(}\frac{\Delta t^{\,4}}{\Delta x^{\,2}}\Bigr{)}\,,$$

Full text

Key words and phrases: moisture diffusion; numerical methods; finite differences; explicit schemes; CFL condition; Dufort–Frankel scheme

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Moisture transfer in porous materials

3. 3 Numerical schemes

4. 3.1 The Crank–Nicolson scheme

5. 3.2 The Euler explicit scheme

6. 3.3 Improved explicit scheme: Dufort–Frankel method

7. 3.4 Hyperbolisation scheme

8. 3.4.1 Dispersion relation analysis

9. 3.4.2 Error estimate

10. 3.4.3 Discretisation

11. 3.5 Validation of the numerical solution

12. 4 Numerical application: linear transfer

13. 5 Extension for nonlinear transfer

14. 5.1 The modified Crank–Nicolson scheme

15. 5.2 The Dufort–Frankel scheme

16. 5.3 Numerical application

17. 5.4 Further nonlinear case studies

18. 5.4.1 Driving rain case

19. 5.4.2 Capillary adsorption case

20. 5.4.3 Results and discussion

21. 6 Conclusion

22. A Stability analysis of the Dufort–Frankel scheme

23. B Dimensionless parameters

24. B.1 Case from Section 4

25. B.2 Case from Section 5.3

26. B.3 Case from Section 5.4

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 [11, 3]. 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:

[TABLE]

associated to boundary and initial conditions, where $\nu$ is the diffusion of the material and where $u\,(\boldsymbol{x},t)$ is the moisture potential being diffused in the spatial domain $\Omega_{\,x}$ during the time interval $\Omega_{\,t}$.We denote $\Delta x$ and $\Delta t$ 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 [18, 9, 2, 28, 26, 12, 13]. However, it has the order of accuracy $\mathcal{O}\,(\Delta t\ +\ \Delta x^{\,2})\,$, while the Crank–Nicolson scheme can be used to increase the accuracy to $\mathcal{O}\,(\Delta t^{\,2}\ +\ \Delta 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 [29, 16]. Nevertheless, this scheme is conditionally stable under the Courant–Friedrichs–Lewy (CFL) condition:

[TABLE]

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.

Moisture transfer in porous materials

The physical problem involves one-dimension moisture diffusion through a porous material defined by the spatial domain $\Omega_{\,x}\ =\ [\,0,\,L\,]$. The moisture transfer occurs according to the liquid and vapour diffusion. The physical problem can be formulated as [12]:

[TABLE]

where $\rho_{\,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}$:

[TABLE]

Thus we have:

[TABLE]

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

[TABLE]

In addition, we have:

[TABLE]

Considering the relation $\rho_{\,l+v}\ =\ \mathrm{f}(\phi)\ =\ \mathrm{f}(P_{\,v},T)\,$, obtained from material properties and from the relation between the vapour pressure $P_{\,v}$ and the relative humidity $\phi\,$, we get:

[TABLE]

Eq. (2.1) can be therefore rewritten as:

[TABLE]

The material properties $\mathrm{f}^{\,\prime}(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:

[TABLE]

where $P_{\,v}^{\,\mathrm{L}}$ and $P_{\,v}^{\,\mathrm{R}}$ are the vapour pressure of the ambient air, $g_{\,l}^{\,\mathrm{L}}$ and $g_{\,l}^{\,R}$ are the liquid flow (driving rain) at the two bounding surfaces. We consider a uniform vapour pressure distribution as initial condition:

[TABLE]

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 $\mathcal{O}(1)$ to avoid severe round-off errors and to likely improve the conditioning of the problem in hands. Therefore, we denote $d_{\,m}\ =\ k_{\,l}\cdot\dfrac{R_{\,v}\,T}{P_{\,v}}\ +\ k_{\,v}$ as a global moisture transport coefficient, $c_{\,m}\ =\ \mathrm{f}^{\,\prime}(P_{\,v})\;\dfrac{1}{P_{\,s}}$ as the moisture storage coefficient and define the following dimensionless quantities:

[TABLE]

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

[TABLE]

Numerical schemes

In order to describe numerical schemes, let us consider a uniform discretisation of the interval $\Omega_{\,x}\ \rightsquigarrow\ \Omega_{\,h}\,$:

[TABLE]

The time layers are uniformly spaced as well $t^{\,n}\ =\ n\,\Delta t\,$, $\Delta t\ =\ \mathrm{const}\ >\ 0\,$, $n\ =\ 0,\,1,\,2,\,\ldots,\,N_{\,t}$ The values of function $u(x,\,t)$ in discrete nodes will be denoted by $u_{\,j}^{\,n}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ u\,(x_{\,j},\,t^{\,n}\,)\,$.

For the sake of simplicity and without loosing the generality, the numerical schemes are explained considering $d_{\,m}^{\,\star}$ and $c_{\,m}^{\,\star}$ as constant, noting $\nu=\dfrac{d_{\,m}^{\,\star}}{c_{\,m}^{\,\star}}$ and the linear diffusion equation:

[TABLE]

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:

[TABLE]

This scheme is $\mathcal{O}(\Delta t^{2}\ +\ \Delta 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}\,$.

The Euler explicit scheme

The standard explicit scheme can be written as:

[TABLE]

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

[TABLE]

where functions $\psi_{\,\mathrm{l},\,\mathrm{r}}\,(\,\bullet\,)$ 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

[TABLE]

By solving Eq. (3.3) with respect to $u_{\,j}^{\,n+1}$, we obtain a discrete dynamical system

[TABLE]

whose starting value is directly obtained from the initial condition:

[TABLE]

It is well-known that scheme (3.3) approximates the continuous operator to order $\mathcal{O}(\Delta t\ +\ \Delta x^{\,2})\,$. The explicit scheme is conditionally stable under the following CFL-type condition:

[TABLE]

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

Improved explicit scheme: Dufort–Frankel method

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

[TABLE]

where the term 2\,u_{\,j}^{\,n}\ is replaced by $u_{\,j}^{\,n-1}\ +\ u_{\,j}^{\,n+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_{\,j}^{\,n+1}$ to give the following discrete dynamical system:

[TABLE]

where:

[TABLE]

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:

[TABLE]

where

[TABLE]

So, from the asymptotic expansion for $\mathcal{L}_{\,j}^{\,n}$ we obtain that the Dufort–Frankel scheme is second order accurate in time and:

• •

First order accurate in space if $\Delta t\ \propto\ \Delta x^{\,3/2}$

• •

Second order accurate in space if $\Delta t\ \propto\ \Delta x^{\,2}$

It is important to note that the boundary condition have to be discretised to the second order of accuracy $\mathcal{O}(\Delta x^{\,2})$ to maintain the method features uniformly.

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:

[TABLE]

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

[TABLE]

The last condition physically means that the new term has only limited influence on the solution of Eq. (3.7). Here $\tau$ is a small ad-hoc parameter whose value is in general related to physical and discretisation parameters $\tau\ =\ \tau\,(\nu,\,\Delta x,\,\Delta 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 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

[TABLE]

3.4.1  Dispersion relation analysis

The classical dispersion relation analysis looks at plane wave solutions:

[TABLE]

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

[TABLE]

The last relation is called the dispersion relation even if the diffusion Eq. (3.1) is not dispersive but dissipative. The real part of $\omega$ contains information about wave propagation properties (dispersive if $\ \frac{\operatorname{Re}\omega(\kappa)}{\kappa}\ \neq\ \mathrm{const}$ and non-dispersive otherwise) while the imaginary part describes how different modes $\kappa$ dissipate (if $\ \operatorname{Im}\omega(\kappa)\ <\ 0$) or grow (if $\ \operatorname{Im}\omega(\kappa)\ >\ 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 $\omega\,$:

[TABLE]

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

[TABLE]

This dispersion relation will be analysed asymptotically with $\tau\ \ll\ 1$ being the small parameter. The branch $\omega_{\,-}\,(\kappa)$ is not of much interest to us since it is constantly damped, i.e.

[TABLE]

It is much more instructive to look at the positive branch $\omega_{\,+}\,(\kappa)\,$:

[TABLE]

The last asymptotic expansion shows that for small values of parameter $\tau\,$, 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 $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:

[TABLE]

Then, the following estimate holds

[TABLE]

where $T\ >\ 0$ is the time horizon and

[TABLE]

and the domain $\Omega_{\,\xi,\,\zeta}$ is defined as

[TABLE]

3.4.3  Discretisation

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

[TABLE]

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 $\mathcal{O}(\Delta t^{\,2}\ +\ \Delta x^{\,2})\,$.

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

[TABLE]

The choice of parameter $\tau$ is therefore an important issue and will be discussed in next Sections.

Validation of the numerical solution

One possible comparison of the numerical schemes can be done by computing the $\mathcal{L}_{\,\infty}$ error between the solution $u_{\,\mathrm{num}}$ and a reference solution $u_{\,\mathrm{ref}}\,$:

[TABLE]

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

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$ $\mathsf{m}$. The moisture properties are $d_{\,m}\ =\ 1.97\cdot 10^{-10}$ $\mathsf{s}$ and $c_{\,m}\ =\ 7.09\cdot 10^{-3}$ $\mathsf{kg/m^{3}/Pa}$. The initial vapour pressure in the material is considered uniform $P_{\,v}^{\,i}\ =\ 1.16\cdot 10^{3}$ $\mathsf{Pa\,}$, corresponding to a relative humidity of $50$%. The reference time is $t^{\,0}\ =\ 1$ $\mathsf{h}$, thus the total time of simulation corresponds to 120 hours, or five days. The boundary conditions, represented by the relative humidity $\phi$ 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\cdot 10^{-7}$ $\mathsf{s/m}$ and $3\cdot 10^{-8}$ $\mathsf{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 $\Delta x\ =\ 10^{-2}$ and $\Delta t\ =\ 10^{-4}$, respecting the CFL condition $\Delta t\,\leqslant\,4.3\cdot 10^{\,-4}$. For the hyperbolisation scheme, $\tau\ =\ \Delta 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$ $\mathsf{h}$ and $t\ =\ 52$ $\mathsf{h\,}$, corresponding to the highest and lowest vapour pressure values. All the numerical methods give accurate results as illustrated with the $\mathcal{L}_{2}$ 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 $\mathcal{O}(\Delta t^{2})$ 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 $\Delta t$. The spatial discretisation is maintained to $\Delta x\ =\ 10^{-2}$ and $\tau\ =\ \Delta t$ for the hyperbolisation scheme. The reference solution has been computed using the Matlab open source package Chebfun [8]. Using the function $\mathrm{pde23t}$, it enables to compute a numerical solution of a partial derivative equation using the Chebyshev functions. Results of the $\mathcal{L}_{\,2}$ error $\varepsilon$ are shown in Figure 5(a). Before the CFL limit, the errors of the Euler, Dufort–Frankel and hyperbolisation schemes are proportional to $\mathcal{O}(\Delta 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 $\mathcal{O}(\Delta t^{2})$. The error of the Dufort–Frankel scheme computed solution becomes too high for $\Delta t\,\geqslant\,3\cdot 10^{\,-2}$. For this case, $\Delta t\ \propto\ \Delta x^{\,\frac{3}{2}}$, therefore, the Dufort–Frankel scheme is first order accurate in space $\mathcal{O}(\Delta 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 $\tau\ =\ \Delta t$, it corresponds to $\Delta t\ \leqslant\ 9\cdot 10^{\,-4}\,$. 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 $\tau$. Figure 5(b) gives the variation of the error $\varepsilon$, for the hyperbolisation scheme, as a function of $\Delta t$, for different values of parameter $\tau$. The error $\varepsilon$ reaches a limit lower than the parameter $\tau\,$. It can be verified that for the choice $\tau\ =\ \nu\,\Delta x\,$, the stability condition corresponds to $\Delta x^{\,\frac{3}{2}}\,$. For this case, the choice $\tau\ \leqslant\ \Delta t$ permits to compute the solution with the best accuracy.

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 $d_{\,m}^{\,\star}\left(\,u\,\right)$ and $c_{\,m}^{\,\star}\left(\,u\,\right)$. 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:

[TABLE]

The modified Crank–Nicolson scheme

The straightforward application of the Crank–Nicolson scheme to Eq. (5.1) yields the following scheme:

[TABLE]

with

[TABLE]

However, this approach leads to deal with nonlinearities due to the evaluation of quantities (as $d_{\,j+{\textstyle{1\over 2}}}^{\,n+1}$) at the upcoming time layer $t\ =\ t^{\,n+1}$. To deal with this issue, linearisation techniques as Picard or Newton–Raphson ones [23, 5], 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 [1]. Thus, the diffusion flux at the interface becomes:

[TABLE]

Finally, the modified Crank–Nicolson schemes yields to:

[TABLE]

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.

The Dufort–Frankel scheme

In the nonlinear case, the Dufort–Frankel numerical schemes is written as:

[TABLE]

The right-hand side term can be expressed as:

[TABLE]

Using the Dufort–Frankel stencil (see Figure 1(c)), the term $u_{\,j}^{\,n}$ is replaced by $\dfrac{u_{\,j}^{\,n+1}\ +\ u_{\,j}^{\,n-1}}{2}\,$. Thus, considering Eq. (5.3), the Dufort–Frankel schemes can be expressed as an explicit scheme:

[TABLE]

with

[TABLE]

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}$, $d_{j+{\textstyle{1\over 2}}}$, $d_{j-{\textstyle{1\over 2}}}$ are explicitly calculated at $t^{\,n}$. It should be noted that the material properties evaluated at $j\ +\ {\textstyle{1\over 2}}$ is formulated as:

[TABLE]

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_{\,v}^{\,i}\ =\ 1.16\cdot 10^{\,3}$ $\mathsf{Pa}$. No 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\cdot 10^{-7}$ $\mathsf{s/m}$ and $3\cdot 10^{-8}$ $\mathsf{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: $\Delta t\ =\ 10^{\,-4}$ and $\Delta 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 $\epsilon\,\leqslant\,0.01\cdot\Delta t^{2}$ 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$ $\mathsf{m}$) and right ($x\ =\ 0.1$ $\mathsf{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$ $\mathsf{h\,}$, due to the nonlinear behaviour of the material. The vapour pressure at $x\ =\ 0.1$ $\mathsf{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 $\mathcal{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 $\Delta t$, maintaining $\Delta x\ =\ 10^{\,-2}$. For each value of $\Delta t$, the $\mathcal{L}_{2}$ 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:

[TABLE]

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 $\mathcal{O}(\Delta t)$. The modification of the Crank–Nicolson scheme, in order to avoid the sub-iterations due to nonlinearities, loose the $\mathcal{O}(\Delta t^{2})$ accuracy.

Even with $\Delta 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 $\Delta 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$ $\mathsf{m\,}$, $\Delta 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 as:

[TABLE]

where $A$ and $\omega$ are the amplitude and the frequency of the signal. The error $\epsilon$ is given in Figure 12 as a function of $\Delta t$ for this case study. If an accuracy of $\epsilon\,\leqslant\,10^{\,-1}$ is required, a time discretisation lower than $\Delta t\,\leqslant\,2\cdot 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:

[TABLE]

The standard Euler implicit and Crank–Nicolson schemes require $N_{\,\mathrm{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 $\mathcal{O}(2\cdot N)$. Generally $\mathcal{O}(N_{\,\mathrm{NL}})\,\gg\,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 $\Delta 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 $\mathcal{O}(N_{\,\mathrm{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.

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.

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$ $\mathsf{m}\,$. The initial vapour pressure is $P_{\,v}^{\,i}\ =\ 1.16\cdot 10^{\,3}$ $\mathsf{Pa}\,$, equivalent to a relative humidity $\phi\ =\ 0.5$. To test the robustness of the scheme, with strong nonlinearities, the properties 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.

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$ $\mathsf{m}$ a moisture flux $g_{\,l}^{\,\mathrm{L}}\ =\ 3.4$ $\mathsf{kg/m^{2}/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 $0.2$, a frequency of $1$h and a mean of $0.5$, at $x\ =\ 0.1$ $\mathsf{m}$. The convective vapour coefficient equals $3\cdot 10^{-8}$ $\mathsf{s/m}$ and the final simulation time is $30$ $\mathsf{h}$.

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 ($\phi\ =\ 1$). In such case, the Robin type boundary condition Eq. (2.3) is modified in order to get a Dirichlet one:

[TABLE]

Here, $P_{\,v}^{\,\mathrm{L}}\ =\ 2.33\cdot 10^{\,3}$ $\mathsf{Pa}$. The vapour pressure of the ambient air is maintained constant at $P_{\,v}^{\,\mathrm{R}}\ =\ 1.16\cdot 10^{\,3}$ $\mathsf{Pa}$, at $x\ =\ 0.1$ $\mathsf{m}$, with a convective vapour coefficient set to $3\cdot 10^{-8}$ $\mathsf{s/m}$. The solution is computed for a final simulation time of $1$ $\mathsf{h}$.

5.4.3  Results and discussion

The discretisations used to compute the solution of both cases are $\Delta x\ =\ 0.01$ and $\Delta t\ =\ 10^{-5}$. The evolution of the vapour pressure at $x\ =\ 0$ $\mathsf{m}$ and $x\ =\ 0.1$ $\mathsf{m}$ is given in Figure 13(a) for the driving rain benchmark. At $x\ =\ 0$ $\mathsf{m}$, the vapour pressure increases due to the constant rain flux $g_{\,l}^{\,\mathrm{L}}$ imposed at the surface. At $x\ =\ 0.1$ $\mathsf{m}$, for $t\leqslant 5$ $\mathsf{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$ $\mathsf{h}$. By $t\ =\ 30$ $\mathsf{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$ $\mathsf{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$ $\mathsf{h}$. The dimensionless moisture transfer coefficient $d_{\,m}^{\,\star}$, represented in Figure 14, highlights the nonlinearities of the material properties. The coefficient has $\mathcal{O}(10^{3})$ orders of variation during the simulation.

The $\mathcal{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 $\mathcal{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 $\epsilon\,\leqslant\,10^{\,-2}\cdot\Delta t^{2}$, 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.

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 $\Delta t$. 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 $\Delta t$ and has the advantage of being unconditionally stable. The second one is also unconditionally stable and second-order accurate in space and time. The latter is second-order accurate in time and in space $\mathcal{O}(\Delta t^{2})$ 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 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 $\tau$ representing the amount of hyperbolicity added in the equation. An optimal choice seems to be $\tau\ =\ \Delta 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 $\Delta t\,$. In addition, for $\Delta t\leqslant 10^{-3}$, the error of the Dufort–Frankel scheme is first-order accurate in time, and for $\Delta t\geqslant 10^{-3}$ second-order accurate in time. The first conclusions revealed that for $\Delta t\leqslant\dfrac{\Delta x^{\,2}}{2\,\nu}$, it is preferable to use the hyperbolisation scheme, for its accuracy. For larger $\Delta 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 $\mathcal{O}(\Delta 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 $\Delta 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 $\Delta 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.

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

[TABLE]

Appendix A Stability analysis of the Dufort–Frankel scheme

The stability analysis is performed for the Initial Value Problem of the linear diffusion equation:

[TABLE]

The Dufort–Frankel numerical scheme yields to:

[TABLE]

Using the von Neumann’s stability analysis [6], we seek for discrete plane wave solutions of the form:

[TABLE]

where $k\ \Delta x$ is the wave length ($k\,\geqslant\,1$), $u_{\,0}$ is the initial wave amplitude ($t\ =\ 0$) and $\alpha$ is the amplification factor between two successive time layers.

In general $\alpha\ =\ \alpha(\,k\,\Delta x\,)$ and we shall compute this factor for the Dufort–Frankel scheme, knowing that:

[TABLE]

Thus, Eq. (A.2) becomes:

[TABLE]

where $\sigma\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}k\,\Delta x$ and $\mu\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\frac{2\,\nu\,\Delta t}{\Delta x^{2}}$. Eq. (A.4) admits two roots:

[TABLE]

In Figures 16(a) and 16(b), we show the dependence of the absolute value of the amplification factors \bigl{|}\,\alpha^{\,\pm}\,\bigr{|} on $\cos(\sigma)$ for two fixed values of the parameter $\mu$ — one below and one above the critical CFL value $\mu\ =\ 1\,$. It is important to notice that \bigl{|}\,\alpha^{\,\pm}\,\bigr{|}\,\leqslant\,1 for all values of the parameter $\sigma\,$. It can be shown (see [24]) that this property holds for all values of $\mu\,$. 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.

Appendix B Dimensionless parameters

Case from Section 4

Problem (2.6) is considered with $g_{\,l,\,\mathrm{L}}^{\,\star}\ =\ g_{\,l,\,\mathrm{R}}^{\,\star}\ =\ 0\,$. The dimensionless properties of the material are equal to $d_{\,m}^{\,\star}\ =\ 1$ and $c_{\,m}^{\,\star}\ =\ 8.6\,$. The final simulation time is fixed to $\tau^{\,\star}\ =\ 120$, while Biot numbers are $\mathrm{Bi}_{\,v}^{\,\mathrm{L}}\ =\ 101.5$ and $\mathrm{Bi}_{\,v}^{\,\mathrm{R}}\ =\ 15.2\,$. The boundary conditions are expressed as:

[TABLE]

Case from Section 5.3

The dimensionless properties of the materials are:

[TABLE]

The Biot numbers are also assumed to be $\mathrm{Bi}_{\,v}^{\,\mathrm{L}}\ =\ 101.5$ and $\mathrm{Bi}_{\,v}^{\,\mathrm{R}}\ =\ 15.2\,$. The ambient water vapour pressure at the boundaries are different from the previous case study. At the left boundary, $u^{\,\mathrm{L}}$ has a fast drop until the saturation state u^{\,\mathrm{L}}\ =\ 2,\;\forall t\in\bigr{[}\,10,\,40\,\bigl{]} and at the right boundary, $u^{\,\mathrm{R}}\ =\ 1\ +\ 0.8\;\sin\left(\,2\pi\,t^{\,\star}\,\right)\,.$

Case from Section 5.4

For both case-studies, an analytic expression of the material properties has been fitted, which dimensionless formulation is:

[TABLE]

The dimensionless parameters for the driving rain case (Section 5.4.1) are:

[TABLE]

while, for the capillary adsorption case (Section 5.4.2), they are:

[TABLE]