An improved explicit scheme for whole-building hygrothermal simulation

TL;DR

This paper introduces an improved explicit Dufort-Frankel scheme for whole-building hygrothermal simulation, offering unconditionally stable, second-order accuracy, and reduced computational costs compared to traditional implicit methods.

Contribution

The paper demonstrates the application of the Dufort-Frankel explicit scheme to coupled heat and moisture transfer, showing it outperforms classical explicit and implicit schemes in stability and efficiency for building simulations.

Findings

Dufort-Frankel scheme is unconditionally stable and second-order accurate.

It reduces computational costs by up to twenty times.

Enables perfect synchronism in whole-building co-simulation.

Abstract

Implicit schemes require important sub-iterations when dealing with highly nonlinear problems such as the combined heat and moisture transfer through porous building elements. The computational cost rises significantly when the whole-building is simulated, especially when there is important coupling among the building elements themselves with neighbouring zones and with HVAC (Heating Ventilation and Air Conditioning) systems. On the other hand, the classical Euler explicit scheme is generally not used because its stability condition imposes very fine time discretisation. Hence, this paper explores the use of an improved explicit approach - the Dufort-Frankel scheme - to overcome the disadvantage of the classical explicit one and to bring benefits that cannot be obtained by implicit methods. The Dufort-Frankel approach is first compared to the classical Euler implicit and explicit schemes to compute the solution of nonlinear heat and moisture transfer through porous materials. Then, the analysis of the Dufort-Frankel unconditionally stable explicit scheme is extended to the coupled heat and moisture balances on the scale of a one- and a two-zone building models. The Dufort-Frankel scheme has the benefits of being unconditionally stable, second-order accurate in time O(dt^2) and to compute explicitly the solution at each time step, avoiding costly sub-iterations. This approach may reduce the computational cost by twenty, as well as it may enable perfect synchronism for whole-building simulation and co-simulation.

Tables

Table 1. Table 1. Hygrothermal properties of the load bearing material.

Property	Value
Material density	${\rho }_{\mathrm{\hspace{0.17em}0}}=790\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Heat capacity	$c_{\mathrm{\hspace{0.17em}0}}=870\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
Water heat capacity	$c_w=4180\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
Sorption isotherm	$f=47.1{\left[1+{\left(-c_{\mathrm{\hspace{0.17em}1}}{P}_c\right)}^{1.65}\right]}^{-0.39}+109.9{\left[1+{\left(-c_{\mathrm{\hspace{0.17em}2}}{P}_c\right)}^{\mathrm{\hspace{0.17em}6}}\right]}^{-0.83}\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Permeability	${\delta }_v=\frac{1.88\cdot {10}^{-6}}{T}\cdot \frac{\left(1-\frac{f}{157}\right)}{0.503{\left(1-\frac{f}{157}\right)}^{\mathrm{\hspace{0.17em}2}}+0.497}\left[𝗌\right]$
Thermal conductivity	$\lambda =0.2+0.0045\cdot f\left[𝖶/\left(𝗆\cdot 𝖪\right)\right]$
Moisture transfer coeff.	$k_M=\exp [-0.07\cdot {\log }_{10}\left(2\frac{P_v}{P_s}\right)+1.44]\cdot 1.97\cdot {10}^{-10}+\exp {[-8\cdot \left(2\frac{P_v}{P_s}\right)-2]}^2\cdot 1.77\cdot {10}^{-7}\left[𝗌\right]$
Moisture storage coeff.	$c_M=-30.62\cdot \left(-c_1\frac{R_{v}T}{P_v}\right)\cdot {\left(-c_1{R}_{v}T{\log }_{10}\frac{P_v}{P_s}\right)}^{0.65}\cdot {\left[1+{\left(-c_1{R}_{v}T{\log }_{10}\frac{P_v}{P_s}\right)}^{1.65}\right]}^{-1.39}-549.5\cdot \left(-c_2\frac{R_{v}T}{P_v}\right)\cdot {\left(-c_2{R}_{v}T{\log }_{10}\frac{P_v}{P_s}\right)}^5\cdot {\left[1+{\left(-c_2{R}_{v}T{\log }_{10}\frac{P_v}{P_s}\right)}^6\right]}^{-1.83}\left[{𝗌}^{\mathrm{𝟤}}/{𝗆}^{\mathrm{𝟤}}\right]$
${}_{}{}^{\star }c_1^{}=1.25\cdot {10}^{-5}$ and $c_2=1.8\cdot {10}^{-5}$.

Table 2. Table 2. Parameter values for the one zone case.

Parameter	North Wall	South Wall	East and West walls
$c_M\left[{𝗌}^{\mathrm{\hspace{0.17em}2}}/{𝗆}^{\mathrm{\hspace{0.17em}2}}\right]$	$1.82\cdot {10}^{-2}$	$1.18\cdot {10}^{-1}$	$6.09\cdot {10}^{-2}$
$k_M\left[𝗌\right]$	$5.89\cdot {10}^{-9}$	$2.92\cdot {10}^{-8}$	$5.47\cdot {10}^{-9}$
$c_{TT}\left[𝖶.𝗌/\left(𝖪.{𝗆}^{\mathrm{\hspace{0.17em}3}}\right)\right]$	$7.7\cdot {10}^{\mathrm{\hspace{0.17em}5}}$	$1.28\cdot {10}^{\mathrm{\hspace{0.17em}6}}$	$8.61\cdot {10}^{\mathrm{\hspace{0.17em}5}}$
$k_{TT}\left[𝖶/\left(𝗆.𝖪\right)\right]$	$2.94\cdot {10}^{-1}$	$8.41\cdot {10}^{-1}$	$3.87\cdot {10}^{-1}$
$c_{TM}\left[𝖶.{𝗌}^{\mathrm{\hspace{0.17em}3}}/\left(\mathrm{𝗄𝗀}.{𝗆}^{\mathrm{\hspace{0.17em}2}}\right)\right]$	$1.52\cdot {10}^{\mathrm{\hspace{0.17em}3}}$	$9.88\cdot {10}^{\mathrm{\hspace{0.17em}3}}$	$5.09\cdot {10}^{\mathrm{\hspace{0.17em}3}}$
$k_{TM}\left[{𝗆}^{\mathrm{\hspace{0.17em}2}}/𝗌\right]$	$1.59\cdot {10}^{-2}$	$2.96\cdot {10}^{-3}$	$1.53\cdot {10}^{-2}$
$A\left[{𝗆}^{\mathrm{\hspace{0.17em}2}}\right]$	$18$	$18$	$9$
$h_T\left[𝖶/\left({𝗆}^{\mathrm{\hspace{0.17em}2}}.𝖪\right)\right]$	$5$	$25$	$12$
$h_M\left[𝗌/𝗆\right]$	$2\cdot {10}^{-7}$	$8\cdot {10}^{-7}$	$4\cdot {10}^{-7}$

Table 3. Table 3. Computer run time required for the numerical schemes.

Wall model with nonlinear wall material properties
Numerical Scheme	CPU time ($𝗌$)	CPU time (%)	Average number of iterations
Dufort–Frankel	150	28	–
Euler explicit	164	31	–
Euler implicit	530	100	3
One-zone model with linear wall material properties
Numerical Scheme	CPU time ($𝗌$)	CPU time (%)	Average number of iterations
Dufort–Frankel	8.5	9	–
Euler explicit	$\infty$	–	–
Euler implicit	95	100	2.95
Two-zones model with nonlinear wall material properties
Numerical Scheme	CPU time ($𝗌$)	CPU time (%)	Average number of iterations
Dufort–Frankel	480	5.4	–
Euler implicit	8900	100	6

Equations

$$\frac{\partial {\rho }_w}{\partial t}$$

$$\displaystyle\bigl{(}\,\rho_{\,0}\ c_{\,0}\ +\ \rho_{\,w}\ c_{\,w}\,\bigr{)}\ \frac{\partial T}{\partial t}\ +\ c_{\,w}\ T\ \frac{\partial\rho_{\,w}}{\partial t}$$

$$P_c$$

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

$$\frac{\partial P_c}{\partial x}=\frac{\partial P_c}{\partial P_v}\frac{\partial P_v}{\partial x}+\frac{\partial P_c}{\partial T}\frac{\partial T}{\partial x}\simeq \frac{{\rho }_l{R}_{v}T}{P_v}\frac{\partial P_v}{\partial x}.$$

$$\frac{\partial {\rho }_w}{\partial t}=\frac{\partial {\rho }_w}{\partial \varphi }\frac{\partial \varphi }{\partial P_v}\frac{\partial P_v}{\partial t}+\frac{\partial {\rho }_w}{\partial T}\frac{\partial T}{\partial t}\simeq \frac{\partial {\rho }_w}{\partial \varphi }\frac{\partial \varphi }{\partial P_v}\frac{\partial P_v}{\partial t}.$$

$$\frac{\partial {\rho }_w}{\partial t}=\frac{f^{\prime }(\varphi )}{P_s}\frac{\partial P_v}{\partial t}.$$

$$k_M\stackrel{\mathrm{def}}{:=}{k}_l\frac{{\rho }_l{R}_{v}T}{P_v}+{\delta }_v$$

$$k_{TM}\stackrel{\mathrm{def}}{:=}{L}_v{\delta }_v$$

$$k_{TT}\stackrel{\mathrm{def}}{:=}\lambda$$

$$c_M\stackrel{\mathrm{def}}{:=}\frac{f^{\prime }(\varphi )}{P_s(T)}$$

$$c_{TT}\stackrel{\mathrm{def}}{:=}{\rho }_0{c}_0+f(\varphi )c_w$$

$$c_{TM}\stackrel{\mathrm{def}}{:=}{c}_{w}T\frac{f^{\prime }(\varphi )}{P_s(T)}$$

$$c_M\frac{\partial P_v}{\partial t}$$

$$c_{TT}\frac{\partial T}{\partial t}+c_{TM}\frac{\partial P_v}{\partial t}$$

$$k_M\frac{\partial P_v}{\partial x}$$

$$k_{TT}\frac{\partial T}{\partial x}+k_{TM}\frac{\partial P_v}{\partial x}$$

$$\displaystyle q_{\,\mathrm{rw}}\ =\ \sum_{w=1}^{m}\mathrm{s}\,\epsilon\,\sigma\,\Biggl{(}\,\biggl{(}\,T_{\,w}(\,x=0\,)\,\biggr{)}^{\,4}\ -\ \biggl{(}\,T(\,x=0\,)\,\biggr{)}^{\,4}\,\Biggr{)}\,,$$

$$P_v$$

$$T$$

$$c_M^{\star }\frac{\partial v}{\partial t^{\star }}$$

$$c_{TT}^{\star }\frac{\partial u}{\partial t^{\star }}+c_{TM}^{\star }\gamma \frac{\partial v}{\partial t^{\star }}$$

$$k_M^{\star }\frac{\partial v}{\partial x^{\star }}$$

$${\mathrm{Fo}}_{TT}{k}_{TT}^{\star }\frac{\partial u}{\partial x^{\star }}+k_{TM}^{\star }{\mathrm{Fo}}_{TM}\gamma \frac{\partial v}{\partial x^{\star }}$$

$$u\stackrel{\mathrm{def}}{:=}\frac{T}{T_i},$$

$$x^{\star }\stackrel{\mathrm{def}}{:=}\frac{x}{L},$$

$$c_{TM}^{\star }\stackrel{\mathrm{def}}{:=}\frac{c_{TM}}{c_{TM,0}},$$

$${\mathrm{Fo}}_M\stackrel{\mathrm{def}}{:=}\frac{t_0\cdot k_{M,0}}{L^2\cdot c_{M,0}},$$

$${\mathrm{Bi}}_M\stackrel{\mathrm{def}}{:=}\frac{h_M\cdot L}{k_{M,0}},$$

$$q_{\infty }^{\star }\stackrel{\mathrm{def}}{:=}\frac{L\cdot q_{\infty }}{T_i\cdot k_{TT,0}}.$$

$$\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{\}}\,.$$

$$\displaystyle\frac{\partial u}{\partial t}\ =\ \boldsymbol{\nabla}\boldsymbol{\cdot}\bigl{(}\,\nu\boldsymbol{\nabla}u\,\bigr{)}\,.$$

$$\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\,,$$

$$\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}}\ +\ \frac{1}{6}\,\Delta t^{\,2}\,\frac{\partial^{3}u}{\partial t^{3}}\\ \ -\ \frac{1}{12}\,\nu\,\Delta x^{\,2}\,\frac{\partial^{4}u}{\partial x^{4}}\ -\ \frac{1}{12}\,\nu\,\Delta t^{\,2}\,\Delta x\,\frac{\partial^{5}u}{\partial x^{3}\,\partial t^{2}}\ +\ \mathcal{O}\Bigl{(}\frac{\Delta t^{\,4}}{\Delta x^{\,2}}\Bigr{)}\,,$$

$$\displaystyle\mathcal{L}_{\,j}^{\,n}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \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}}\,.$$

$$\displaystyle\frac{u_{\,j}^{\,n+1}\ -\ u_{\,j}^{\,n-1}}{2\,\Delta t}\ =\ \frac{1}{\Delta x}\Biggl{[}\,\left(\,\nu\frac{\partial u}{\partial x}\,\right)_{\,j+{\textstyle{1\over 2}}}^{\,n}\ -\ \left(\,\nu\frac{\partial u}{\partial x}\,\right)_{\,j-{\textstyle{1\over 2}}}^{\,n}\,\Biggr{]}\,.$$

$$\frac{1}{\Delta x}\left[{\left(\nu \frac{\partial u}{\partial x}\right)}_{j+\frac{1}{2}}^n-{\left(\nu \frac{\partial u}{\partial x}\right)}_{j-\frac{1}{2}}^n\right]=\frac{1}{\Delta x^2}\left[{\nu }_{j+\frac{1}{2}}^n{u}_{j+1}^n+{\nu }_{j-\frac{1}{2}}^n{u}_{j-1}^n-\left({\nu }_{j+\frac{1}{2}}^n+{\nu }_{j-\frac{1}{2}}^n\right)u_j^n\right].$$

Full text

Key words and phrases: heat and moisture transfer; numerical methods; finite-differences; explicit schemes; Dufort–Frankel scheme; whole-building simulation

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Porous building element hygrothermal model

3. 3 Numerical schemes

4. 3.1 Improved explicit scheme: Dufort-Frankel method

5. 3.2 Validation of the numerical solution

6. 4 Numerical application: porous wall transfer

7. 5 Whole-building hygrothermal model

8. 5.1 Coupling the lumped air multizone model with the transfer through porous walls

9. 5.2 Implicit scheme for the whole-building energy simulation: problem statement

10. 5.3 Improved explicit schemes for the whole-building energy simulations

11. 6 Numerical application: multizone transfer

12. 6.1 One-zone case study

13. 6.2 Two-zones case study

14. 7 Conclusion

15. Acknowledgments

16. A The Dufort-Frankel scheme for weakly coupled equations

17. B Dimensionless numerical values

18. B.1 Whole-building model with linear material properties

19. Nomenclature

Introduction

Models for the combined heat and moisture transfer through porous building elements have been implemented in building simulation tools since the 1990’s in software such as Delphin [2], MATCH [22], MOIST [5], WUFI [9] and UMIDUS [17, 19, 7] among others. More recently, those models have 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 [28].

The Euler and Crank–Nicolson implicit schemes have been used in many studies and implemented in building simulation tools, as reported in the literature [18, 9, 2, 24, 23, 13, 14, 27, 16], due to their numerical property of unconditional stability. 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, as when material properties are moisture content or temperature dependent, one has to perform sub-iterations to linearise the system, increasing the total CPU time. In [13], thousands of sub-iterations are reported to converge to the solution of a mass diffusion problem. Another disadvantage of implicit schemes appears when coupling the wall model, representing the transfer phenomena in porous building elements, to the room air model. The wall and the room air models must iterate within one time step until reaching a given tolerance [12]. If it does not impose any limitation on the choice of the time discretisation, it induces sub-iterations that increase the computational time of the simulation of the whole-building model. Moreover, it is valuable to decrease this computational cost knowing that the hygrothermal and energy building simulation is generally carried out for time scale periods as long as one year, or even more. However, the phenomena and particularly the boundary conditions evolve a time scale of seconds.

Recently, in [10], the improved Dufort–Frankel explicit scheme was explored for the solution of moisture diffusion equation highlighting that the standard stability limitation can be overcome, which inspired to investigate the use of the Dufort–Frankel scheme for the solution of the combined heat and moisture transfer through porous building elements coupled with room air models.

In this way, this paper first describes in Section 2 the heat and moisture transfer model. In Section  3, basics of the Dufort–Frankel explicit scheme is detailed before exploring the features of the scheme applied to linear and nonlinear cases, presented in Section 4. Then, the benefits of using an unconditionally stable explicit scheme are investigated to perform a whole-building hygrothermal simulation based on coupling the porous element model to a room air multizone model.

Porous building element hygrothermal model

The physical problem considers one-dimensional heat and moisture transfer through a porous material defined by the spatial domain $\Omega_{\,x}\ =\ [\,0,\,L\,]\,$. The following convention is adopted: $x\ =\ 0$ corresponds to the surface in contact with the inside room and, $x\ =\ L\,$, corresponds to the outside surface. The moisture transfer occurs due to capillary migration and vapour diffusion. The heat transfer is governed by diffusion and latent mechanisms. The physical problem can be formulated as [13, 25]:

[TABLE]

where $\rho_{\,w}$ is the volumetric moisture content of the material, $\delta_{\,v}$ and $k_{\,l}\,$, the vapour and liquid permeabilities, $P_{\,v}\,$, the vapour pressure, $T\,$, the temperature, $R_{\,v}\,$, the water vapour gas constant, $P_{\,c}$ the capillary pressure, $c_{\,0}\,$, the material heat capacity, $\rho_{\,0}\,$, the material density, $c_{\,w}$ the water heat capacity, $\lambda$ the thermal conductivity, and, $L_{\,v}$ the latent heat of evaporation. Eq. (2.1a) 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}\,$, and the Clausius–Clapeyron equation:

[TABLE]

Neglecting the variation of the capillary pressure and the mass content with temperature [23], the partial derivative of $P_{\,c}$ can be written as:

[TABLE]

In addition, we have:

[TABLE]

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

[TABLE]

We denote by:

[TABLE]

Considering the previous notation, Eq. (2.1) can be rewritten as:

[TABLE]

The boundary conditions at the interface between the porous material and the air are expressed as:

[TABLE]

where $P_{\,v,\,\infty}$ and $T_{\,\infty}$ stand for the vapour pressure and temperature of the air and $h_{\,M}$ and $h_{\,T}$ are the convective transfer coefficients. If the bounding surface is in contact with the outside building air, $g_{\,\infty}$ is the liquid flux from wind driven rain and $q_{\,\infty}$ is the total heat flux from radiation and the heat contribution from the inward liquid water penetration. If the bounding surface is in contact with the inside building air, $g_{\,\infty}\ =\ 0$ and $q_{\,\infty}$ is the distributed short-wave radiative heat transfer rate $q_{\,\mathrm{rw}}$ in the enclosure and long-wave radiative heat exchanged among the room surfaces:

[TABLE]

where $s$ is the view factor between two surfaces, $\sigma$ is the Stefan–Boltzmann constant, $\epsilon$ is the emissivity of the wall surface, $w$ represents the $m$ bounding walls. We consider a uniform vapour pressure and temperature distributions as initial conditions:

[TABLE]

The governing equations can be written in a dimensionless form as:

[TABLE]

and the boundary condition as:

[TABLE]

where the dimensionless quantities are defined as:

[TABLE]

The dimensionless formulation enables to determine important scaling parameters (Biot and Fourier 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 [20]. 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.

Numerical schemes

Let us consider a uniform discretisation of the interval $\Omega_{\,x}\ \rightsquigarrow\ \Omega_{\,L_{\,x}}\,$:

[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 the 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 generality, simple diffusion equation is considered:

[TABLE]

First, the numerical schemes are explained considering the linear case. Then, the extension to the nonlinear case is described.

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}\ from the explicit scheme is replaced by $u_{\,j}^{\,n-1}\ +\ u_{\,j}^{\,n+1}\,$. The Scheme (3.2) has the stencil depicted in Figure 1. At a first glance, the scheme (3.2) looks like an implicit scheme, however, it is not truly the case. Eq. (3.2) 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 shows that the Dufort–Frankel scheme is unconditionally stable [10, 21, 26]. The consistency error analysis of the Scheme (3.2) shows the following result:

[TABLE]

where:

[TABLE]

So, from the asymptotic expansion for $\mathcal{L}_{\,j}^{\,n}$ one can see 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}$ .

In the nonlinear case, the  numerical scheme can be derived as follows:

[TABLE]

The right-hand side term can be expressed as:

[TABLE]

Using the Dufort–Frankel stencil (see Figure 1), the term $u_{\,j}^{\,n}$ is replaced by $\dfrac{u_{\,j}^{\,n+1}\ +\ u_{\,j}^{\,n-1}}{2}\,$.

Thus, considering Eq. (3.4), 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 layer $n\,$, the material properties $\nu_{\,j+{\textstyle{1\over 2}}}\,$, $\nu_{\,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]

Validation of the numerical solution

All the numerical results in this paper were computed using Matlab. One possibility for comparing the computed solutions is by computing the $\mathcal{L}_{\,\infty}$ error between the solution $u_{\,\mathrm{num}}$ and a reference solution $u_{\,\mathrm{ref}}\,$:

[TABLE]

The reference solution is computed using the Matlab open source package Chebfun [8]. Using pde functions, it enables to compute a numerical solution of a partial derivative equation with the Chebyshev polynomials adaptive spectral methods.

The $\mathcal{L}_{\,\infty}$ error can be computed along the space or time domains, according to:

[TABLE]

Numerical application: porous wall transfer

Heat and moisture transfer are strongly nonlinear due to the variation of the material properties with the field. For this reason, this case study is devoted to investigated these effects by means of the Dufort–Frankel numerical scheme. A physical application is performed in a $10-\mathsf{cm}$ slab of bearing material. To test the robustness of the scheme, the properties of the material are gathered from [13, 11, 10] and given in Table 1. At $t\ =\ 0\ \mathsf{h}\,$, the material is considered with uniform fields, with a temperature of $T_{\,i}\ =\ 20\ \mathsf{{}^{\,\circ}C}$ and relative humidity of $\phi_{\,i}\ =\ 50\ \mathsf{\%}\,$. The boundary conditions, represented by the relative humidity $\phi$ and temperature $T\,$, are given in Figure 2. The convective mass and heat transfer coefficients are set to $h_{\,M}\ =\ 2\cdot 10^{\,-7}\ \nicefrac{{\mathsf{s}}}{{\mathsf{m}}}\,$, $h_{\,T}\ =\ 25\ \nicefrac{{\mathsf{W}}}{{\mathsf{(m^{\,2}\,.\,K)}}}$ and to $h_{\,M}\ =\ 3\cdot 10^{\,-8}\ \nicefrac{{\mathsf{s}}}{{\mathsf{m}}}\,$, $h_{\,T}\ =\ 8\ \nicefrac{{\mathsf{W}}}{{\mathsf{(m^{\,2}\,.\,K)}}}\,$, for the left and right boundary conditions, respectively. The dimensionless numerical values of this case are provided in B.

The numerical solution was computed using three different schemes, considering a spatial discretisation parameter of $\Delta x^{\,\star}\ =\ 10^{\,-2}\,$. For the time domain, the discretisation parameter equals $\Delta t^{\,\star}\ =\ 10^{\,-3}$ for the Dufort–Frankel and Euler implicit scheme. A tolerance $\eta\leqslant 10^{\,-2}\ \Delta t^{\,\star}$ has been used for the convergence of the sub-iterations of the implicit scheme, using a fixed-point algorithm. For the Euler explicit scheme, a time discretisation $\Delta t^{\,\star}\ =\ 10^{\,-5}$ has been used, to respect the CFL stability condition. The latter has been computed for the heat transfer Equation (2.5a):

[TABLE]

and for the mass transfer Equation (2.5b):

[TABLE]

which give the numerical values of $\Delta t^{\,\star}\ \leqslant\ 3\cdot 10^{\,-4}$ and $\Delta t^{\,\star}\ \leqslant\ 7.4\cdot 10^{\,-5}\,$, corresponding to physical values of $\Delta t^{\,\star}\ \leqslant\ 1.08\ \mathsf{s}$ and $\Delta t^{\,\star}\ \leqslant\ 0.2\ \mathsf{s}\,$, respectively. It can be noted that the stability of the coupled equation is given by the moisture equation stability.

The time evolution of the fields in the middle of the material is presented in Figure 3, showing a good agreement between the solution of each numerical scheme. The coupling effect between heat and mass transfer can be specially noted in Figure 3(a). The temperature at $x\ =\ 0.05\ \mathsf{m}$ varies according to both frequencies of the left side air temperature and relative humidity, illustrated in Figures 2(c) and 2(a). The relative humidity rises until $80\ \mathsf{\%}$ at $t\ =\ 110\ \mathsf{h}\,$, demonstrating that the material wont be saturated in the middle during the simulation. Figure 4 represents the field profiles at different time. We can observe that this material is non-hygroscopic but fairly permeable. It diffuses rapidly but do not retain the moisture.

The errors between the reference solution and the ones computed with the different numerical schemes are given in Figures 5(b) and 5(a). It confirms that all the numerical schemes enable to compute an accurate solution, at the order of $10^{\,-4}$ for both fields and considered spatial and temporal discretisations.

For the Euler implicit and the Dufort–Frankel explicit schemes, the CPU time has been calculated, using Matlab platform on a computer with Intel i$7$ CPU and $32$GB of RAM, and shown in Table 3. The implicit scheme requires around $9$ sub-iterations per time step to treat the nonlinearities of the problem. The Dufort–Frankel approach computes directly the solution and therefore has a reduced computational costs, around $15\%$ of the Euler implicit scheme based algorithm.

A numerical analysis of the behaviour of the three numerical schemes has been carried out for different values of the temporal discretisation $\Delta t^{\,\star}\,$. The spatial discretisation is maintained to $\Delta x^{\,\star}\ =\ 10^{\,-2}\,$. Results of the $\mathcal{L}_{\,\infty}$ error $\varepsilon$ are shown in Figure 6. As expected, the Euler explicit scheme was not able to compute a solution when the stability CFL condition is not respected (around $\Delta t^{\,\star}\ \leqslant\ 8\cdot 10^{\,-5}$). The values computed from Eq.  (4.1) and (4.2) are in accordance with the results from the convergence study. It also confirms that the Dufort–Frankel scheme is unconditionally stable, as it computes a solution for any discretisation parameter $\Delta t^{\,\star}\,$. The error of the numerical scheme is second-order accurate in time $\mathcal{O}\,(\Delta t^{\,2})\,$.

The choice of the discretisation parameters have strong influences on the solution computed. When analysing the error of the Dufort–Frankel scheme as a function of the time discretisation $\Delta t^{\,\star}\,$, three regions have been highlighted. The first one corresponds to small discretisation parameter $\Delta t^{\,\star}\,\leqslant\,10^{\,-3}$, where the solution obtained with the Dufort–Frankel scheme reaches the constant error value which is lower than the Euler implicit one. In this region, all schemes provide an accurate solution. The second region is where the error is proportional to $\mathcal{O}\,(\Delta t^{\star\,2})$ but higher than the Euler implicit one. The last region, when $\Delta t^{\,\star}\,\geqslant\,5\cdot 10^{\,-2}\,$, includes the large time step discretisation. With this in mind, the first conclusion we can make is that until the second region the Dufort–Frankel scheme is more recommended, while in the third region, for bigger values of $\Delta t^{\,\star}$ the Euler implicit scheme is suggested. However, when the time discretisation $\Delta t^{\,\star}$ is too large, both schemes do not succeed in representing the physical phenomena. As mentioned in [10], the time step has to be carefully chosen in accordance with the characteristic time of the physical phenomena. For the Dufort–Frankel scheme and the case studied, this condition is reached for time discretisation up to $180$-$\mathsf{s}\,$.

Whole-building hygrothermal model

Coupling the lumped air multizone model with the transfer through porous walls

The lumped multizone model divides the whole-building into $N_{\,z}$ perfectly mixed air zones. For each zone, the evolutions of the air temperature $T_{\,a}$ and the humidity ratio $w_{\,a}$ are given by the equations [25, 3, 4]:

[TABLE]

Besides the dependence of the air heat capacity on the humidity ratio, the right-hand terms of the room air energy balance equation — Eq. (5.1a) — can also be strongly dependent on the humidity ratio, which may lead to a highly nonlinear problem.

Equation (5.1b) can be expressed using the vapour pressure of the air zone $P_{\,v,\,a}\,$, by means of the following relation:

[TABLE]

where $P_{\,v}^{\,\circ}\ =\ 1.61\cdot 10^{\,5}\ \mathsf{Pa}\,$. In addition, according to the notation used before, we have:

[TABLE]

The terms $G$ and $Q$ are associated to occupants and their activities, denoted with the subscript $o\,$, as well as air flow from ventilation systems or infiltration, denoted with the subscript $v\,$. The sources due to interzone airflow are designated by $Q_{\,\mathrm{inz}}$ and $G_{\,\mathrm{inz}}\,$. Interested readers may find a detailed list of these sources and their physical description in [7]. The first can be expressed as:

[TABLE]

where $g_{\,\mathrm{o}}$ is the time variant vapour production. The sources due to ventilation system are written as:

[TABLE]

where $g_{\,\mathrm{v}}$ stands for the air flow rate due to ventilation and infiltration, $T_{\,\infty}$ and $w_{\,\infty}$ are the outside temperature and humidity ratio, and $c_{\,w}$ the water heat capacity depending on the air temperature and humidity ratio. In a similar way, the sources due to an airflow $g_{\,\mathrm{inz}}$ between zones $i$ and $j$ are given by:

[TABLE]

Here, a mean value of the latent heat of evaporation $L_{\,v}\,$, evaluated at the considered temperatures, is used for numerical application. The terms $Q_{\,w}$ and $G_{\,w}$ represent the mass and heat quantities exchanged between the air room and the the $N_{\,w}$ bounding walls:

[TABLE]

where $A$ is the wall surface area. The terms $T_{\,s}$ and $P_{\,v,\,s}$ are the wall surface temperature and vapour pressure. Considering Eq. (2.2), they correspond to $T_{\,s}\ =\ T\,(x\ =\ 0)$ and $P_{\,v,\,s}\ =\ P_{\,v}\,(x\ =\ 0)$ (with the adopted convention). Thus, the coupling between the wall model Eq. (2.2) and the lumped multizone model Eq. (5.1) is operated by the sources $G_{\,w}$ and $Q_{\,w}$ from Eq. (5.2).

The temperature and vapour pressure in the zone are initially at $T\ =\ T_{\,i}$ and $P_{\,v}\ =\ P_{\,v,\,i}\,$. To write a dimensionless formulation of Eq. (5.1), the following quantities are defined:

[TABLE]

Thus, the dimensionless formulation of Eq. (5.1) can be written as:

[TABLE]

The quantities $u_{\,\infty}\,$, $v_{\,\infty}$ and $u_{\,a,\,2}\,$, $v_{\,a,\,2}$ come from building outside (provided by weather data) and from the adjacent zone. The coupling between the wall model Eq. (2.5) and the lumped multizone model Eq. (5.3) is operated through the wall source terms. The dimensionless parameters $\mathrm{Bi}_{\,TT}\,$, $\mathrm{Bi}_{\,TM}$ and $\mathrm{Bi}_{\,M}$ qualify the penetration of the heat and moisture through the wall according to the physical mechanism. The parameters $\theta_{\,T,\,i}$ and $\theta_{\,M,\,i}$ depend on the wall surface on the room air, providing the weighted contribution of the wall to the energy and moisture balances.

Implicit scheme for the whole-building energy simulation: problem statement

For the sake of simplicity and without loosing the generality, the coupling procedure is explained considering only the linear heat diffusion equation Eq. (3.1) for one wall model with the following boundary conditions for the surface in contact with the outside ($x\ =\ 0$) and inside($x\ =\ 1$) air of the building:

[TABLE]

The room air energy conservation equation for the multizone model is expressed as:

[TABLE]

Many whole-building models reported in the literature, such as MATCH [22], MOIST [5] and DOMUS [17, 19, 6], use implicit (Euler or Crank–Nicolson) approaches to solve these equations, mainly due to the unconditional stability property. Thus, using an Euler implicit approach, the discretisations of Eqs. (3.1), (5.4a), (5.4b) yield to:

[TABLE]

Therefore, at iteration $t^{\,n}\,$, solution $u^{\,n\,+\,1}\,$, of the wall model, is computed using the resolvent operator $\mathcal{R}_{\,\mathrm{imp}}^{\,\mathrm{w}}$ written as:

[TABLE]

In the same way, the discretisation of Eq. (5.5) yields to:

[TABLE]

and, at iteration $t^{\,n}\,$, solution $u_{\,a}^{\,n\,+\,1}\,$, of the multizone model, is computed using the resolvent operator $\mathcal{R}_{\,\mathrm{imp}}^{\,\mathrm{a}}$ written as:

[TABLE]

By coupling the wall and the multizone models to perform a whole-building energy simulation, ones must solve at each time iteration the following system of equations:

[TABLE]

As the system Eq. (5.6) is nonlinear, it is not possible to solve it directly. Thus, \bigl{(}\,u^{\,n\,+\,1},\,u_{\,a}^{\,n+1}\,\bigr{)} are computed using, for instance, a fixed point algorithm, until reaching a prescribed tolerance $\eta\,$, as illustrated in Alg. 1. The number of subiterations should increase with the number of walls considered in the whole-building model, when dealing with nonlinear models and when considering interzone airflows. Interested readers may report to [6] for more details on this problem statement. Indeed, authors investigate the influence of the time step on the computational time of a whole-building simulation model, considering only heat transfer, for different numerical schemes and analytical solution.

Improved explicit schemes for the whole-building energy simulations

Using the improved Dufort–Frankel explicit scheme, results from Section 3.1 have shown that the discretisation of Eqs. (3.1), (5.4a), (5.4b) is given by:

[TABLE]

Moreover, to solve Eq. (5.5), we get the scheme:

[TABLE]

Therefore, by coupling the wall and the multizone models, at each time iteration, the solution of the equations is directly computed:

[TABLE]

Thus, the use of improved explicit schemes, such as the Dufort–Frankel one, is particularly advantageous to avoid subiterations at each time step to compute the dependent variable fields. In addition, the numerical property of being unconditionally stable avoids any limitation on the values of the time step $\Delta t\,$, which is chosen only according to the characteristic time of the physical phenomena [10].

Numerical application: multizone transfer

One-zone case study

A single zone, surrounded by four walls is considered, as illustrated in Figure 7. All the walls have $0.1\ \mathsf{m}$ of length constituted by the same material, the load bearing material, whose properties are given in Table 1. Each wall is consider to have linear transfer and their coefficients are computed as follows: the north wall with $23^{\,\circ}\mathsf{C}$ and $15\ \%\,$, the south wall with $23^{\,\circ}\mathsf{C}$ and $90\ \%$ and the east and west walls with $23^{\,\circ}\mathsf{C}$ and $50\ \%\,$. The values of the linear coefficients are shown in Table 2. Their dimensionless values are provided in B. The convective mass and heat transfer coefficients are set to $h_{\,M}\ =\ 3\cdot 10^{\,-8}\ \nicefrac{{\mathsf{s}}}{{\mathsf{m}}}$ and $h_{\,T}\ =\ 8\ \nicefrac{{\mathsf{W}}}{{\mathsf{(m^{\,2}\,.\,K)}}}$ between the inside air zone and the walls. It is assumed that no transfer occurs through the ceiling and the floor. The zone has a floor area of $18\ \mathsf{m^{\,2}}$ and a volume of $54\ \mathsf{m^{\,3}}\,$. The room is subjected to a constant moisture load of $25\ \mathsf{g/h}$ with an increase of $400\ \mathsf{g/h}$ from $6\ \mathsf{h}$ to $9\ \mathsf{h}$ each day. External air enters the room through ventilation and infiltrations at a constant rate of $g_{\,v}\ =\ 0.5\ \mathsf{h^{\,-1}}\,$. There are no radiative heat exchange among the walls. The initial temperature and relative humidity in the zone are $T_{\,i}\ =\ 20\ \mathsf{{}^{\,\circ}C}$ and $\phi_{\,i}\ =\ 50\ \mathsf{\%}\,$, respectively. The outside boundary conditions $T_{\,\infty}$ and $\phi_{\,\infty}$ are given in Figures 8(a) and 8(b). The convective transfer coefficients between the walls and the outside conditions are also provided in Table 2. For this case, neither radiation or rain sources are considered for the boundary conditions.

The whole simulation is performed for $80\ \mathsf{h}\,$. The walls and zone fields are computed using the Euler implicit and Dufort–Frankel explicit schemes. The Euler explicit scheme has a stability restrictions imposing a very small time discretisation ($\Delta t^{\,\star}\ \leqslant\ 10^{\,-5}$ for the given space discretisation) and therefore an important computer run time. Moreover, to our knowledge, this scheme is not commonly used in building simulation programs. Thus, it was not used in this case study. The space and time discretisation parameters are $\Delta x^{\,\star}\ =\ 10^{\,-2}$ and $\Delta t^{\,\star}\ =\ 10^{\,-3}\,$, corresponding, from a physical point of view, to $\Delta x\ =\ 10^{\,-3}\ \mathsf{m}$ and $\Delta t\ =\ 3.6\ \mathsf{s}\,$. Each solution is compared to a reference solution computed using the Chebfun open-source package for Matlab as explained in Section 3.2.

The time evolution of the temperature and relative humidity for different walls and for the air zone are given in Figures 9 and 10. The solutions computed with the Euler and Dufort–Frankel schemes are in very good agreement with the reference. As for the material properties of the wall, differences between the field evolution can be observed. Particularly, in Figures 9(a) and 9(c), the increase and decrease of temperature at $x\ =\ 0\ \mathsf{m}$ are higher for the Eastern wall than for the Northern one. It is due to a higher Biot number for the latter: $\mathrm{Bi}_{\,TT}\ =\ 3.1$ against $\mathrm{Bi}_{\,TT}\ =\ 1.7\,$. In Figures 10(a) and 10(b) the increase and decrease of relative humidity are different at $x\ =\ 0\ \mathsf{m}$ between the Northern and Southern walls. In Figure 10(d), a slight increase of relative humidity due to moisture generation in the zone from $6$ to $9\ \mathsf{h}$ can be observed.

The $\mathcal{L}_{\,\infty}$ error with respect to space $x$ has been computed for the fields in the four walls, as illustrated in Figure 11(a). Moreover, the $\mathcal{L}_{\,\infty}$ error with respect to $t$ for the fields in the zone is given in Figure 11(b). The error is approximately $\mathcal{O}\,(10^{\,-4})$ for both models, proving the accuracy of the solution computed with the Euler implicit and Dufort–Frankel explicit schemes.

For each scheme, the CPU time has been calculated, using the Matlab platform on a computer with Intel i$7$ CPU and $32$GB of RAM, which is shown in Table 3. The Dufort–Frankel explicit scheme requires less than 10% of the time needed for the Euler implicit one, to compute the solution. This difference is due to the sub-iterations needed to solve the nonlinear system composed by equations of the wall and zone models. In this case, a fixed-point algorithm with tolerance parameter $\eta\ =\ 0.01\cdot\Delta t^{\,\star}\,$, as the one illustrated in Algorithm. 1, has been used. As emphasized in Figure 12, this algorithm required almost three sub-iterations to compute the solution at each time iteration even in the linear case.

A convergence study of the whole-building model has been carried out by fixing the space discretisation to $\Delta x^{\,\star}\,=\,10^{\,-2}$ and varying the time discretisation. Results are reported in Figures 13(a) and 13(b), for each field. The error with the reference solution has been computed for each field (temperature and relative humidity) and for each model (wall and zone). It can be noted that the wall model reach a constant absolute accuracy lower than the one for the zone. As discussed in the previous case, the Dufort–Frankel scheme is preferable to use for $\Delta t^{\,\star}\,\leqslant\,5\cdot 10^{\,-2}$ the equivalent to $180$-$\mathsf{s}\,$. For bigger time step values the implicit scheme would be more advisable at the strict condition that it respects the characteristic time of the physical phenomena [10]. Moreover, in building simulation, the recommended time step is about $100$-$\mathsf{s}$ [6] which makes the Dufort–Frankel scheme suitable for this kind of application.

Two-zones case study

Previous case study considered a single zone building hygrothermal simulation with linear wall material properties. It enabled to enhance the sub-iterations required using Euler implicit scheme when coupling the wall and zone models. The present case focus on a bi-zone building, as shown in Figure 14, and take into account the variation of the wall material properties with temperature and relative humidity, such as the ones used in Section 4. Furthermore, the radiative heat exchanged among the room air surfaces are included. All the inside walls have an emissivity $\epsilon\ =\ 0.5$ except the wall $7\,$, which have a higher value $\epsilon\ =\ 0.9$. The view factor is set to $s\ =\ 0.2\,$. Incident radiation is considered for the outside boundary condition, shown in Figure 15. An airflow of $g_{\,\mathrm{inz}}\ =\ 0.3\ \mathsf{h^{\,-1}}$ occurs between both zones. Only zone $1$ is subjected to moisture sources due to occupants and to air ventilation. Zone $2$ receives a heat source $q_{\,h}\ =\ 500\ \mathsf{W}$ from $3$ to $4\ \mathsf{h}\,$. The other parameters have the same numerical values as in the previous case study. This case study was designed in order to enforce the nonlinearities of the whole-building model, to increase the sub-iterations of the Euler implicit approach and thus enhance the efficiency of the Dufort–Frankel explicit scheme. The solution is computed using both schemes, with a time and space discretisations $\Delta t^{\,\star}\ =\ 10^{\,-3}$ and $\Delta x^{\,\star}\ =\ 10^{\,-2}\,$, equivalent to $\Delta x\ =\ 10^{\,-3}\ \mathsf{m}$ and $\Delta t\ =\ 3.6\ \mathsf{s}\,$.

Figures 16(a) and 16(b) present the time evolution of the dependent variable fields for the wall $1$, showing a very good agreement among the three solutions. By comparing Figures 9(a), 10(a), 16(a) and 16(b), the effect of considering nonlinear material properties can be observed. Similar observations can be done by analysing Figures 17(a), 17(b), 17(c) and 17(d). The moisture production in zone $1$ is enhanced by an increase of the relative humidity from $6$ to $9\ \mathsf{h}\,$. As the moisture generation only occurs in this zone, there is no increase of the relative humidity in zone $2$ due to this phenomena. Moreover, the heat production in zone $2$ is highlighted from $2$ to $6\ \mathsf{h}\,$. The temperature reach $24^{\circ}\,\mathsf{C}\,$. A slight increase is observable in zone $1$ during this period due to the heat generation in zone $2$ and to the airflow between both zones. The error with the reference solution is of the order of $\mathcal{O}\,(10^{\,-4})$ for the wall and zone fields, as shown in Figures 18(a) and 18(b).

By considering nonlinear wall material properties, the Euler implicit scheme require more sub-iterations at each time step as illustrated in Figure 19. In the previous linear case study, the algorithm required around $3$ iterations whereas for the present case, it needs at least $6$ to achieve the same accuracy. Therefore, the CPU time of the Euler implicit scheme to compute the numerical solution increases, as reported in Table 3. As no sub-iterations are necessary for the Dufort–Frankel explicit scheme, the computation gain rises compared to the previous case study. It needs only $5\%$ of the CPU time of the Euler implicit scheme. These gains might considerably increase when considering highly nonlinear phenomena such as driving rain and iteration with HVAC systems [1].

The purpose of this study was essentially to compare the numerical methods and the innovative Dufort–Frankel explicit scheme on a nonlinear case of whole building hygrothermal simulation. The computational time was measured for different simulation time horizon:

[TABLE]

It can be noted that the increase of the CPU time is linear for both approaches. Therefore for a simulation time horizon of $1$ year, the CPU ratio between the Dufort–Frankel scheme and the Euler implicit approach would be conserved. It is important to compare the ratio between the different approaches and not only the absolute value of CPU time with the classical building simulation programs. Indeed, the algorithms developed in this study are not optimised and the computing platforms may be different.

Conclusion

Implicit methods are extensively used in building simulation codes due to their stability conditions. However, they may require important extra computation when dealing with highly nonlinear problems, such as the combined heat and moisture transfer through porous building elements or when the whole-building is simulated, demanding perfect synchronism. In this way, this study aimed at exploring the use of the improved Dufort–Frankel explicit scheme. The advantages of this approach are studied here in whole-building hygrothermal simulations, with one- and two-zones building models.

The solution computed with the Dufort–Frankel scheme was compared to three other different solutions. Results have shown that the Dufort–Frankel scheme enables to compute an accurate solution, which is second-order accurate in time — $\mathcal{O}\,(\Delta t^{\,2})\,$. It also enhanced that the Dufort–Frankel scheme can overcome some disadvantages of the Euler implicit and explicit approaches.

First, contrarily to the Euler explicit scheme, it is unconditionally stable, enabling to compute the solution without satisfying any stability condition, although, it does not mean that any value of $\Delta t$ can be used. When analysing the error of the Dufort–Frankel scheme as a function of the discretisation in the time domain — $\Delta t\,$, three regions can be highlighted. The first one corresponds to small values of $\Delta t\,$, where the solution obtained with the Dufort–Frankel scheme reaches a constant error value, which is lower than the Euler implicit one. The second region is where the error is proportional to $\mathcal{O}\,(\Delta t^{\,2})$ but higher than the Euler implicit. The third and last region includes large values of $\Delta t\,$ where the Euler implicit scheme seems to be more accurate. However, both schemes do not succeed in representing the physical phenomena. As mentioned in [10], the time step has to be chosen in accordance with the characteristic time of the problem. For the case studied, this condition is reached for a time discretisation up to $180-\mathsf{s}\,$, which is in accordance with the precision given in [6] for whole building simulation.

Then, when dealing with nonlinearities, the scheme does not require any sub-iteration at each time step. The solution is directly computed, reducing consequently the computational cost of the algorithm. For the first case study, the Euler implicit scheme required around $3$ sub-iterations, making it three times more costly than the Dufort–Frankel scheme. When coupling the wall and zone models using implicit schemes, a nonlinear system of equations has to be solved. By using the Dufort–Frankel explicit scheme, the system of equations becomes linear and no sub-iterations are necessary. Therefore, within the Dufort–Frankel approach, the algorithm requires only $9\%$ of the CPU time of the implicit approach. When considering the nonlinearities of the wall material properties and long-wave radiative heat transfer among the room surfaces, the computational savings rise to $95\%\,$.

These results encourage to apply the Dufort–Frankel approach in building simulation tools. The computational gains should increase with the number of rooms, walls, partitions, and furniture. The proposed method enable perfect synchronism for simulation and co-simulation, which can reduce even more the computation efforts.

Acknowledgments

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.

Appendix A The Dufort-Frankel scheme for weakly coupled equations

This Appendix details the application of the Dufort–Frankel approach for two weakly coupled equations, as the ones describing the heat and moisture transfer in porous materials (Section 2). For this purpose, we consider a uniform discretisation as described in Section 3 and the following linear coupled diffusion equation:

[TABLE]

The boundary conditions associated to Eq. (A.1) are written as:

[TABLE]

Considering Eq. (A.1a) and the straightforward application of the Dufort–Frankel scheme described in Section 3.1, we get:

[TABLE]

where:

[TABLE]

For Eq. (A.1b), the numerical scheme is expressed as:

[TABLE]

With Eq. (A.3) and rearranging the terms of Eq. (A), the numerical scheme is derived as follows for the field $u\,$:

[TABLE]

where:

[TABLE]

For the boundary conditions, the application of the Dufort–Frankel scheme to Eq. (A.2a) gives:

[TABLE]

Here, the node $j\ =\ 0$ is a ghost one located a distance $\Delta x$ from the node $j\ =\ 1$. From Eq. (A.6), we can deduce $v_{\,0}\,$:

[TABLE]

In a similar way for Eq. (A.2b), we get:

[TABLE]

Using Eqs. (A.7) and (A), it is possible to compute Eqs (A.3) and (A) for $j\ =\ 1\,$. A similar approach is adopted for node $j\ =\ N\,$.

Appendix B Dimensionless numerical values

This Appendix provides the dimensionless values of the linear case study considered in this work.

Whole-building model with linear material properties

The dimensionless properties of the wall materials, considered in Section 6.1, are:

[TABLE]

and for all the walls:

[TABLE]

At $x\ =\ 0$, the Biot numbers are:

[TABLE]

and at $x\ =\ 1$

[TABLE]

For the zone model, the properties are $\kappa_{\,TT,\,0}\ =\ 1$, $\kappa_{\,TT,\,1}\ =\ 1.43\cdot 10^{\,-2}$ and the coupling parameter $\theta$:

[TABLE]

The source term due to ventilation system equals:

[TABLE]

The source term due to moisture load equals:

[TABLE]

[TABLE]

The outside boundary conditions are expressed as:

[TABLE]

The final simulation time is fixed to $\tau^{\,\star}\ =\ 80\,$.

Nomenclature

[TABLE]

[TABLE]

[TABLE]

[TABLE]