On the comparison of three numerical methods applied to building simulation: finite-differences, RC circuit approximation and a spectral method

TL;DR

This paper compares three numerical methods—finite-differences, RC circuit approximation, and spectral method—for building physics simulations, evaluating their accuracy and computational efficiency across different heat and moisture transfer scenarios.

Contribution

It provides a comparative analysis of the performance of three numerical methods in building physics, highlighting the spectral method's superior accuracy.

Findings

Spectral method is the most accurate among the three.

RC approach with few resistances is less accurate for temperature and vapor pressure.

Finite-difference method offers a balance between accuracy and complexity.

Abstract

Predictions of physical phenomena in buildings are carried out by using physical models formulated as a mathematical problem and solved by means of numerical methods, aiming at evaluating, for instance, the building thermal or hygrothermal performance by calculating distributions and fluxes of heat and moisture transfer. Therefore, the choice of the numerical method is crucial since it is a compromise among (i) the solution accuracy, (ii) the computational cost to obtain the solution and (iii) the complexity of the method implementation. An efficient numerical method enables to compute an accurate solution with a minimum computational run time (CPU). On that account, this article brings an investigation on the performance of three numerical methods. The first one is the standard and widely used finite-difference approach, while the second one is the so-called RC approach, which is a particular method brought to the building physics area by means of an analogy of electric circuits. The third numerical method is the spectral one, which has been recently proposed to solve nonlinear diffusive problems in building physics. The three methods are evaluated in terms of accuracy on the assessment of the dependent variable (temperature or vapor pressure) or of density of fluxes for three different cases: i) heat diffusion through a concrete slab, ii) moisture diffusion through an aerated concrete slab and iii) heat diffusion using measured temperatures as boundary conditions. Results highlight the spectral approach as the most accurate method. The RC based model with a few number of resistances does not provide accurate results for temperature and vapor pressure distributions neither to flux densities nor conduction loads.

Tables

Table 1. Table 1: Efficiency of the numerical models for linear heat diffusion.

Numerical Model	$\mathrm{scd}[-]$	${\epsilon }_{\infty }[-]$ for $T$	${\epsilon }_{\infty }[-]$ for $q$	$R_{\mathrm{cpu}}\left[\mathrm{𝗆𝗌}/𝗁\right]$
R$2$C	$1.5$	$0.02$	$0.47$	$0.03$
R$3$C	$1.8$	$0.02$	$0.45$	$0.06$
R$100$C	$3.6$	$8\cdot {10}^{-5}$	$6\cdot {10}^{-3}$	$12$
FDM	$4.2$	$7\cdot {10}^{-5}$	$1.7\cdot {10}^{-3}$	$15$
Spectral	$4.3$	$4\cdot {10}^{-5}$	$6\cdot {10}^{-4}$	$1.23$

Table 2. Table 2: Efficiency of the numerical models for nonlinear moisture diffusion.

Numerical Model	$\mathrm{scd}[-]$	${\epsilon }_{\infty }[-]$ for $P_v$	${\epsilon }_{\infty }[-]$ for $g$	$R_{\mathrm{cpu}}\left[\mathrm{𝗆𝗌}/𝗁\right]$
R$2$C	$0.69$	$0.20$	$10.7$	$0.02$
R$3$C	$0.72$	$0.12$	$10.6$	$0.11$
R$100$C	$3.5$	$2\cdot {10}^{-4}$	$5\cdot {10}^{-2}$	$38$
FDM	$3.4$	$3.7\cdot {10}^{-4}$	$8\cdot {10}^{-2}$	$37$
Spectral	$4.2$	$2\cdot {10}^{-5}$	$4\cdot {10}^{-3}$	$12.3$

Equations

$$\rho c\frac{\partial T}{\partial t}-k\frac{{\partial }^{2}T}{\partial x^2}=0,$$

$$T=T_L(t),$$

$$T=T_R(t),$$

$$\displaystyle T\ =\ T_{\,0}\,,\qquad t\ =\ 0\,,\qquad\forall\ x\ \in\ \bigl{[}\,0\,,\,L\,\bigr{]}\,.$$

$$q(t)$$

$$q_R(t)$$

$$E\stackrel{\mathrm{def}}{:=}{\int }_{t_1}^{t_2}q(\tau )\mathrm{d}\tau ,$$

$$\displaystyle\xi\,(P_{\,v})\,\frac{\partial P_{\,v}}{\partial t}\ -\ \frac{\partial}{\partial x}\,\biggl{(}\,\kappa\,(P_{\,v})\,\frac{\partial\,P_{\,v}}{\partial x}\,\biggr{)}\ =\ 0\,,$$

$$P_v=P_{v,L}(t),$$

$$P_v=P_{v,R}(t),$$

$$\displaystyle P_{\,v}\ =\ P_{\,v,\,0}\,,\qquad t\ =\ 0\,,\qquad\forall\ x\ \in\ \bigl{[}\,0\,,\,L\,\bigr{]}\,.$$

$$g(t)$$

$$u=\frac{T}{T_0},$$

$$x^{\star }=\frac{x}{L},$$

$$F{o}_T=\frac{k{t}^{\circ }}{\rho c{L}^2}.$$

$$\frac{\partial u}{\partial t^{\star }}=F{o}_T\frac{{\partial }^{2}u}{\partial x^{\star 2}},$$

$$u=u_R,$$

$$u=u_L,$$

$$u=1,$$

$$v=\frac{P_v}{P_{v,0}},$$

$$x^{\star }=\frac{x}{L},$$

$${\kappa }^{\star }=\frac{\kappa }{{\kappa }^{\circ }},$$

$$F{o}_m=\frac{{\kappa }^{\circ }{t}^{\circ }}{{\xi }^{\circ }{L}^2}.$$

$$\displaystyle\xi^{\,\star}\,(\,v\,)\,\frac{\partial v}{\partial t^{\,\star}}\ =\ \,Fo_{\,m}\,\frac{\partial}{\partial x}\,\biggl{(}\,\kappa^{\,\star}\,(\,v\,)\,\frac{\partial v}{\partial x^{\,\star\,2}}\,\biggr{)}\,,$$

$$v=v_R,$$

$$v=v_L,$$

$$v=1,$$

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

$$\frac{\mathrm{d}{u}_j}{\mathrm{d}t}=Fo\frac{u_{j-1}-2u_j+u_{j+1}}{\Delta x^2},j=1,\dots ,N-1,n⩾0,$$

$$u_j(0)=1.$$

$${\int }_{x_{j-1/2}}^{x_{j+1/2}}\rho c\frac{\partial T}{\partial t}\mathrm{d}x={\int }_{x_{j-1/2}}^{x_{j+1/2}}k\frac{{\partial }^{2}T}{\partial x^2}\mathrm{d}x.$$

$$T_j\stackrel{\mathrm{def}}{:=}\frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}T(x,t)\mathrm{d}x.$$

Full text

a Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LOCIE, 73000 Chambéry, France

b Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAMA, 73000 Chambéry, France

c Thermal Systems Laboratory, Mechanical Engineering Graduate Program,

Pontifical Catholic University of Paraná, Rua Imaculada Conceição, 1155, CEP : 80215-901, Curitiba - Paraná, Brazil

*∗*Corresponding author. E-mail address: julien.berger@univ-smb.fr

keywords: Heat Transfer; Moisture Transfer; Numerical Methods; Finite Differences; Thermal Circuit Model; Spectral Method.

1 Introduction

As the building sector represents almost $33$% of the world global energy consumption, current environmental issues lead to focus on energy efficiency of building envelopes Administration, (2015). Within this context, several tools have been developed since the $1970$’s for the accurate assessment of building energy performance. Many of them have been reported in the frame of the International Energy Agency Annex $41$ published by Woloszyn and Rode in Woloszyn and Rode, (2008) and more recently in Mendes et al. Mendes et al., (2016).

Among all the physical phenomena involved in building energy efficiency, energy losses associated to heat and moisture transfer through the building envelope are of major importance. They represent an important part of building energy consumption and moisture may considerably impact on conduction loads and on the size of HVAC systems Mendes et al., (2003), besides promoting severe disorders when reaching high levels Berger et al., (2015). Thus, it is of primary importance to have numerical models enabling to accurately represent the physical phenomena for the evaluation of heat losses and gains through building envelopes.

The numerical models are elaborated from the main governing equations representing the heat and/or moisture transfer in building porous materials detailed for instance in Mendes et al., (2016). The use of analytical solutions is often limited due to the nonlinearity of the material properties and to the non-periodicity of the boundary conditions. Thus, most of models referenced in literature are based on numerical approaches.

Consequently, the main challenge arises in elaborating efficient numerical models to perform the simulation. The word efficiency can designate several features. One decisive aspect is the accuracy of the computed solution. It is of capital importance for the design of energy efficient buildings. To predict reliable energy consumption, one must be certain that the numerical errors of the model are negligible. When comparing the predictions to experimental observations, such as performed in Yang and Becerik-Gerber, (2015), researchers often assume that the numerical errors are always lower than the uncertainties present in the measurements, in the inputs parameters and in the mathematical model that described the physical phenomena. Another important feature is the computational run time of the numerical model to compute the solution. Even with the increase of computer power in the recent decades, it is still a crucial issue. The numerical model needs to save the computational efforts to ease the work of building designers and engineers. It is also relevant in the research context of sensitivity analysis or parameter estimation problem, where a large number of direct model computations is required.

Surprisingly, despite the widespread use of models in research and practice, the efficiency of the numerical models have received little attention in the literature. Models are implemented in software such as EnergyPlus Crawley et al., (2001), ESP-r Clarke, (2013), BSim Rode and Grau, (2003), etc. The governing equations are well detailed but the description of the numerical methods is generally brief and no discussion on their efficiency are provided. Moreover, the validity of the assumption that the numerical error can be negligible has never been verified. Thus, this article proposes to overcome this issue by presenting a detailed evaluation of the efficiency of the three numerical methods.

The first one, is the standard finite-differences based approach, which is probably the most-used method to compute the solution of the diffusion problem. Different variations of this approach have been reported in the literature such as the implicit Euler in Mendes and Philippi, (2005); Steeman et al., (2009), the explicit Euler in Tariku et al., (2010); Kalagasidis et al., (2007) or the Crank–Nicolson in Van Genuchten, (1982). The second method is the so-called RC approach. This method was first used during the second world war nad J. Lee et al., (1943) where analogous electric networks were built to solve the solution of transient heat-flow problems. Since there were no computer devices, it was an ingenuous way to rapidly simulate the solution of the problem. Interesting details can be found in Lawson and McGuire, (1953); Robertson and Gross, (1958) with an investigation of the error devices as a function of the number of (physical) thermal resistances. Although the appearance of digital computers started in the fifties and the rapid and progressive hardware evolution since the seventies, this RC method is still used in many algorithms to solve the partial differential equation of heat transfer as for instance in Fraisse et al., (2002); Roels et al., (2017); Naveros and Ghiaus, (2015). The third method is more advanced spectral method which was recently applied for the solution of diffusion problems through porous building elements Gasparin et al., (2018, 2017).

The manuscript is organized as follows. The physical problem of heat and moisture transfer in building porous materials is recalled in Section 2. The three numerical methods are described in Section 3. Then, three case studies are analysed. The first one, in Section 4, considers a linear heat transfer in a concrete wall. Then, in Section 5, a nonlinear case of moisture diffusion is investigated. In Section 6, the efficiencies of the numerical methods are evaluated considering real measured temperatures as boundary conditions. Some conclusion and final remarks are outlined in Section 8.

2 Physical problem and mathematical formulation

2.1 Physical phenomenon of heat transfer

The physical problem involves heat conduction in a wall of thickness $L$ composed of a single material in which the thermal conductivity is denoted as $k$, the density as $\rho$ and the specific heat as $c\,$. The problem can be formulated by the Fourier (or heat) equation, for x\ \in\ \bigl{[}\,0\,,\,L\,\bigr{]} and t\ \in\ \bigl{[}\,0\,,\,\tau\,\bigr{]}\,:

[TABLE]

where $T(\,x\,,\,t\,)$ is the temperature within the wall at the distance x\ \in\ \bigl{[}\,0\,,\,L\,\bigr{]} and time $t\ \geqslant\ 0\,$.

For the sake of simplicity111Different approaches are reported in literature to represent the Neumann or Robin boundary conditions within the RC model framework. Thus, to limit the possible sources of error, the numerical investigation was performed considering Dirichlet boundary conditions. , Dirichlet boundary conditions are assumed at the extremity of the wall:

[TABLE]

At the initial state, the temperature of the wall is assumed to be uniform:

[TABLE]

One of the interesting outputs in the building physics framework is the heat flux density at x_{\,0}\ \in\ \bigl{[}\,0\,,\,L\,\bigr{]}, defined as:

[TABLE]

Particularly, we denote as $q_{\,R}$ the heat flux density computed at $x\ =\ L$:

[TABLE]

The conduction loads represent the heat fluxes at the building envelope internal surface and, in terms of the energy density, it can be evaluated as:

[TABLE]

which can be evaluated, for instance, over daily or monthly periods.

2.2 Physical phenomenon of moisture transfer

The moisture transfer occurs under isothermal conditions in a wall of thickness $L\,$, with a single material of permeability $\kappa$ and moisture capacity $\xi$, both depending on the vapor pressure. The formulation of the problem, for x\ \in\ \bigl{[}\,0\,,\,L\,\bigr{]} and t\ \in\ \bigl{[}\,0\,,\,\tau\,\bigr{]}\,, yields to:

[TABLE]

where $P_{\,v}$ is the vapor pressure within the wall.

The boundary conditions at the extremity of the wall are:

[TABLE]

A uniform vapor pressure is assumed as initial condition:

[TABLE]

The vapor flux density is similarly computed according to:

[TABLE]

2.3 Dimensionless formulation

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 such as the Biot and Fourier numbers. 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 Nayfeh, (2000). Finally, the floating point arithmetics is designed such as the rounding errors are minimal if computer manipulates the numbers of the same magnitude Kahan and Palmer, (1979). Moreover, the floating point numbers have the highest density within the interval $(\,0,\,1\,)$ and their density decays exponentially when we move further away from zero. Figure 1 shows the accuracy of the floating points in Matlab™ environment. It is generated using the $\mathrm{eps}(\,x\,$) function in the Matlab environment. For a numerical model written using dimensionless equations, the accuracy of the floating points scales with $\mathcal{O}(\,10^{\,-17}\,)\,$. The potential commonly used in the physical model of heat and mass transfer is usually the temperature $T$ in $\ [\,\mathsf{K}\,]$, the vapor pressure $P_{\,v}$ or the capillary pressure $P_{\,c}$ both pressures in $\ [\,\mathsf{Pa}\,]\,$. According to Figure 1, if the numerical model is written using the temperature or the vapor pressure with their physical dimension, the accuracy of the floating points scales between $10^{\,-14}$ and $10^{\,-13}\,$, respectively. The range of the capillary pressure is between $10^{\,3}\ [\,\mathsf{Pa}\,]$ and $10^{\,8}\ [\,\mathsf{Pa}\,]\,$. Therefore, the accuracy of the floating point can loose up to $8$ orders compared to a dimensionless numerical model. So, it is always better to manipulate numerically the quantities of the order of $\mathcal{O}(1)$ to avoid severe round-off errors and to likely improve the conditioning of the problem in hands.

2.3.1 Heat transfer

In this way, according to Incropera et al., (2007), we define the following dimensionless quantities for the temperature:

[TABLE]

The time and space domains are also modified through a unitless formulation:

[TABLE]

where $t^{\,\circ}$ is a characteristic time. The Fourier dimensionless number is defined, characterizing the importance of the heat transfer through the material:

[TABLE]

Therefore, the unitless system of the differential equation of heat transfer is formulated as:

[TABLE]

together with the boundary conditions:

[TABLE]

and the initial condition:

[TABLE]

2.3.2 Mass transfer

In a very similar way, we define the following dimensionless quantities related to the vapor pressure field Luikov, (1966):

[TABLE]

The unitless formulation of the time and space domains are:

[TABLE]

The vapor pressure dependent moisture properties are transformed according to:

[TABLE]

where $\kappa^{\,\circ}$ and $\xi^{\,\circ}$ are reference property values. The Fourier dimensionless number is defined, characterizing here the importance of the moisture transfer through the material:

[TABLE]

The Fourier number quantifies the first order of the diffusion transfer, while the dimensionless parameters $\kappa^{\,\star}$ and $\xi^{\,\star}$ define the distortion or the nonlinearity of the phenomenon.

The unitless system of differential equation for moisture transfer is:

[TABLE]

with the boundary conditions:

[TABLE]

and the initial condition:

[TABLE]

3 Numerical methods

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

[TABLE]

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

For the sake of simplicity, the standard finite-differences scheme and spectral approach will be described for the linear dimensionless heat diffusion equation (5). To our knowledge, the RC model is always described in the literature considering the physical dimension of the equation. In this way, the approach will be presented considering the heat diffusion equation (1).

3.1 The standard finite-differences method

The standard semi–discrete scheme based on central finite-differences can be written as:

[TABLE]

whose starting value is directly obtained from the initial condition:

[TABLE]

Many approaches can be used for the temporal discretisation of Eq. (7). Here the algorithm is implemented in Matlab™ environment using the function ode45 providing an efficient explicit Runge–Kutta scheme. In whole figures, this approach will be referenced as FDM.

3.2 The RC model

Both sides of the heat equation (1) is integrated over $x$ for the cell illustrated in Figure 2(a):

[TABLE]

The average temperature of the cell is defined as

[TABLE]

Thus, Eq. (8) becomes:

[TABLE]

From the electric analogy of the heat conduction (Davies, 2004, Chap. 10), the heat flux is approximated by:

[TABLE]

where $R$ is the thermal resistance defined as:

[TABLE]

Finally, for each node $j$, the temperature is computed using:

[TABLE]

$C$ being the thermal capacity defined as:

[TABLE]

The electric analogy of the heat conduction equation is illustrated in Figure 2(b). It can be noted that Eq. (9) corresponds to the central finite-differences discretisation of the second space derivative. The RC model states that the temperature can be computed using a user-defined number r\,\in\bigl{\{}\,1\,,\,\ldots\,,\,N\,\bigr{\}} of resistances. Usually, such a model is denoted as $R_{\,r}\,C\,$ with $r$ of the order of the unity $r\,\simeq\,\mathcal{O}(\,1\,)\,$. This hypothesis corresponds to compute the temperature using $r\ -\ 1$ points of discretisation. It can also be seen as a low fidelity model to represent the physical phenomena of heat or moisture transfer in building porous material.

In this work, an Euler explicit approach is associated with the semi-discrete scheme (9), in agreement with Biddulph et al., (2014). In this case, it is important to note that the explicit scheme is only conditionally stable under the following Courant–Friedrichs–Lewy-type condition Courant et al., (1928):

[TABLE]

The algorithm is implemented in the Matlab™ environment. Interested readers are invited to consult Davies, (2004); Fraisse et al., (2002) for more details on this approach and Deconinck and Roels, (2016); Reynders et al., (2014); Jimenez et al., (2009) for examples of applications in building physics.

3.3 The advanced spectral method

The spectral method has different approach, other than the central differences and thus, the RC ones. It assumes that the unknown $u\,(\,x,\,t\,)$ from Eq. (5) can be approximately represented as a finite sum Mendes et al., (2016):

[TABLE]

Here, $\{\mathsf{T}_{\,i}\,(\,x\,)\}_{\,i\,=\,0}^{\,n}$ is a set of basis functions that remains constant in time. In this study, the Chebyshev polynomials are used as the basis functions since they are optimal in $\mathcal{L}_{\infty}$ approximation norm Gautschi, (2004). The functions $\{a_{\,i}\,(\,t\,)\}_{\,i\,=\,0}^{\,n}$ are the corresponding time-dependent spectral coefficients. The parameter $n$ represents the number of degrees of freedom of the solution, also denoted as the order of the solution with $n\,\simeq\,\mathcal{O}(\,10\,)\,\,$. The main advantage of the spectral method is that $n\,\ll\,p$, where $p$ is the number of degrees of freedom needed to solve problem (5) by means of conventional methods such as finite-differences, finite-volume or finite-element methods. For these reasons, the spectral method is also denoted as the spectral-Reduced Order Method (spectral-ROM) Gasparin et al., (2018, 2017).

The derivatives are written as follows:

[TABLE]

where the dot denotes $\dot{a}_{\,i}\,(\,t\,)\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{\mathrm{d}a\,(\,t\,)}{\mathrm{d}t}$ according to Newton notation. Using the properties of the Chebyshev polynomials, the space derivatives are re-expanded in the same Chebyshev basis function. The connection is explicitly given from the recurrence relation of the Chebyshev polynomial derivatives Peyret, (2002):

[TABLE]

with,

[TABLE]

Using the expression of the derivatives provided by Eqs. (12b) and (12c), the residual of the diffusion equation (5) is:

[TABLE]

which is considered a misfit of the approximate solution. The purpose is to minimize the residual:

[TABLE]

which is realised via the the Tau–Galerkin method, which requires Eq. (13) to be orthogonal to the Chebyshev basis functions $\langle\,R\,,\mathsf{T}_{\,i}\,\rangle\,=\,0$. Here, the the scalar product is defined by :

[TABLE]

Thus, it leads to the following relations among spectral coefficients:

[TABLE]

Finally, after the projection and expansion of the residual, the original partial differential equation (5) is reduced to a system of ordinary differential equations plus two algebraic expressions enabling to compute the time dependent coefficients $\{a_{\,i}\,(\,t\,)\}\,$. For linear problems, the system of ordinary differential equations is explicitly built:

[TABLE]

where $\mathcal{A}\in\mathrm{Mat}\,_{n\times n}(\mathds{R})\,$, with constant coefficients, $\mathrm{b\,}(\,t\,)\in\mathds{R}^{\,n}$ is a vector resulting from boundary conditions and $a_{\,0}$ is the vector of initial coefficients. Initial values of the coefficients $\{a_{\,i}\,(t\,=\,0)\}$ are calculated by the Galerkin projection of the initial condition Canuto et al., (2006):

[TABLE]

where $u_{\,0}\,(\,x\,)$, is the dimensionless initial condition. Interested readers may refer to Gasparin et al., (2018, 2017) for further details on the spectral method.

3.4 Extension of the methods for nonlinear problems

The extension of the three methods for nonlinear problem is now detailed. For the standard semi–scheme based on central finite-differences for Eq. (6), it is formulated as:

[TABLE]

where

[TABLE]

For the RC model, the extension for nonlinear problem in its physical dimension Eq. (4) is given by:

[TABLE]

where $R_{\,j}$ is the vapor resistance and $C_{\,j}$ being the moisture capacity, both defined respectively as:

[TABLE]

For the spectral method, Eq. (6) is rearranged as follows:

[TABLE]

where,

[TABLE]

As described in Section 3.3, the unknown $v\,(\,x\,,t\,)$ is approximated by the finite sum (11) with Chebyshev polynomials as basis functions. The derivatives are written as in the linear case, by Eqs. (12a), (12b) and (12c). Substituting them into Eq. (15), we get:

[TABLE]

Then, the Tau–Galerkin method is used to minimize the residual of the equation. The integrals of the nonlinear coefficients $\nu\,(\,v\,)$ and $\lambda\,(\,v\,)$ are computed using the Chebyshev–Gauß quadrature. At the end, it results in a system of Differential-Algebraic Equations (DAEs) with the following form:

[TABLE]

where, $\mathcal{M}$ is a diagonal and singular matrix containing the coefficients of the Chebyshev weighted orthogonal system, $\mathrm{b\,}\,(\,t\,)$ is a vector containing the boundary conditions and, $\mathcal{A}\cdot a_{\,n}\,(\,t\,)$ is composed by the right member of Eq. (3.4) projected on the Chebyshev basis functions. The initial condition is given by Eq. (14) and the DAE system is solved by ode15s or ode23t from Matlab™.

3.5 Methods implementation and metrics of their efficiency

All the numerical algorithms for the classic and advanced schemes are written using dimensionless variables, while for the RC model, the physical dimensional variables are used. The numerical solutions are computed using an adaptive time step $\Delta t$ using Matlab™ function ode45 Shampine and Reichelt, (1997) with an absolute and relative tolerances set to $10^{\,-4}$. The efficiencies of the method are evaluated in terms of three criteria: (i) the global error of the numerical solution, (ii) the significant digits of the solution and (iii) the computational run time to compute the solution.

To evaluate the error, a reference solution $u^{\,\mathrm{ref}}\,(\,x\,,t\,)$ is computed using a numerical pseudo–spectral approach obtained with the Matlab™ open source toolbox Chebfun Driscoll et al., (2014). Using the function pde23t, it permits to compute a numerical solution of a partial derivative equation using the Chebyshev functions. This useful package enables to compute reference solutions for one dimensional space-time problems. This tools is chosen since it can deal with more complex problems than analytical solution. Indeed, it can consider nonlinear coefficients or Robin-type time-dependent boundary conditions, which are more related to building physics application. The error between the solution, obtained by the numerical methods described above, and the reference one is computed as a function of $x$ by the following formula:

[TABLE]

where $N_{\,t}$ is the number of temporal steps. The global uniform error $\L_{\,\infty}$ is given by the maximum value of $\varepsilon_{\,2}\,(\,x\,)\,$:

[TABLE]

The significant correct digits of the solution are evaluated according to Soderling and Wang, (2006):

[TABLE]

As the RC approach computes directly the fields in their physical dimension, a scaling transformation is performed to compute the errors. The last criteria is the computational (CPU) run time required by the numerical model to compute the solution. It is measured using the Matlab™ environment with a computer equipped with Intel i$7$ CPU and $32$ GB of RAM. We define the ratio:

[TABLE]

where $t_{\,\mathrm{cpu}}\ \ [\,\mathsf{s}\,]$ is the CPU time and $\tau$ is the final physical time of the simulation.

4 Numerical investigations: linear heat diffusion

Since the three methods have been presented, a first linear case of heat diffusion is considered to evaluate their efficiency, considering the following input values for a concrete slab:

[TABLE]

The boundary conditions are defined as sinusoidal variations:

[TABLE]

The temperature is computed using three RC model approaches with the thermal resistances r\in\bigl{\{}\,2\,,\,3\,,\,100\,\bigr{\}}\,. In addition, the problem is solved using the standard finite-differences method and spectral approach with $N\ =\ 6$ modes. A spatial discretisation step $\Delta x\ =\ 10^{\,-3}\ \mathsf{m}$ is considered for both approaches. The temperature profiles at the last time of the simulation are shown in Figures 3 (a,b). It can be noted that the approaches with two or three resistances cannot compute an accurate temperature profile. A perfect agreement is observed between the reference and the solution computed using the R$100$C, the standard finite-differences and the spectral approaches. The time evolution of the temperature in the middle of the wall is shown in Figures 4(a) and 4(b). Apparently, it seems that each approach enables to represent the temperature evolution. However, as shown in Figure 4(c), the difference with the reference solution can reach $0.2\ [\,\mathsf{{}^{\circ}C}\,]$ for the approach with two resistances. For the spectral approach, the difference is of the order $\mathcal{O}(\,10^{\,-3}\,)\,$.

One can argue that the differences, between the reference temperature and the one computed with the RC approach using two or three resistances, are acceptable. However, the discrepancy increases drastically for the heat flux density – which is directly dependent on the temperature derivative at the boundary –, going up to $100\,\%$ as highlighted in Figure 5(a). For the RC model with $100$ resistances, the heat flux density is computed with a very satisfying accuracy. As expected, similar results are observed for the standard finite-differences approach. The heat flux density is however more accurate when computed with the spectral approach. Figures 6(a) and 6(b) shows the error $\varepsilon_{\,2}$, computed using the dimensionless fields. It confirms that the temperature, computed with an RC model approach with only two or three resistances, lacks of accuracy to represent the physical phenomenon of heat diffusion. It should be noted that the error of the spectral method scales with $\mathcal{O}(\,10^{\,-5}\,)$ with only $N\ =\ 6$ modes. For the flux, the best accuracy reaches $\mathcal{O}(\,10^{\,-2}\,)\,$, which is obtained with the spectral approach. The error of the RC model with two or three resistances is completely unacceptable. It can be noted that the accuracy is more sensible to the heat flux density. This can also be understood by analyzing the propagation of the numerical errors. According to the definition in Eq. (2), the finite-difference approximation of the heat flux density is given by

[TABLE]

If we consider that the temperature is computed with a numerical perturbation $T\simeq T\ +\ \delta\,T$, with \bigl{|}\bigl{|}\,\delta\,T\,\bigr{|}\bigr{|}\ \ll\ 1\,, the approximation of the heat flux is:

[TABLE]

Since the perturbations are uncorrelated, $\delta\,T_{\,j}\ -\ \delta\,T_{\,j-1}\simeq 2\,\delta T$ and the approximation of the heat flux becomes

[TABLE]

Since $\Delta x\ \ll\ 1\,$, the error on the heat flux density is higher than the one on the temperature. Considering the numerical values read in Figure 6(a), the numerical perturbation scales with $\delta T\ =\ \mathcal{O}(\,10^{\,-4}\,)$ for the standard finite-difference approach. With the spatial discretisation $\Delta x\ =\ 10^{\,-2}\,$, the term $\displaystyle\frac{2\,\delta T}{\Delta x}$ is of the order $\mathcal{O}(\,10^{\,-2}\,)\,$. Thus, the error on the flux cannot be lower than this value. This value is in consistent with the results observed in Figure 6(b).

In terms of digits accuracy, the results are reported in Table 1. It is noticed the RC approach with $100$ resistances computes the field with more than twice digits accuracy than the R$2$C one. The spectral method presents the highest accuracy. For the computational time, as expected, the numerical models with $2$ and $3$ resistances are the fastest approach. Indeed they require a very few computations at each time step ($3$ and $4\,$, respectively). However, the speed of computation decreases with the accuracy of the solution. The computational effort of the standard finite-difference approach scales with the RC one for $100$ resistances. A good compromise between speed of computation and accuracy is the Spectral approach.

These results highlight that the RC model with a small number of resistance cannot provide an accurate solution. A natural question rises: how many resistances are required to accurately compute the temperature within the wall? The answer to this question strongly depends on the numerical values of the application. For this case study, the error with the reference solution has been computed as a function of the resistance number $r$, as shown in Figure 7. If the field of interest is the temperature, it can be noted that a number $r\ =\ 10$ of resistances is sufficient. The error remains stable at $\mathcal{O}(\,10^{\,-4}\,)$ corresponding to the tolerance set in the ode45 solver of the Matlab™ environment. Although, if one is interested in the heat flux density, the minimum number is $r\ =\ 30\,$ to reach an error lower than $10^{\,-2}\,$.

5 Numerical investigations: nonlinear moisture diffusion

The previous section considered a linear model of diffusion. It is important to evaluate the performance of the methods for nonlinear problems since the accuracy of the solution can be deteriorated. For this, a nonlinear moisture diffusion problem is considered. The length of the wall is set to $L\ =\ 0.1\ [\,\mathsf{m}\,]$ and the simulation horizon to $t\ =\ 72\ [\,\mathsf{h}\,]\,$. The material properties are inspired from Bednar and Hagentoft, (2005) with a constant moisture capacity and a vapor pressure dependent moisture permeability:

[TABLE]

Sinusoidal variations are imposed as boundary conditions:

[TABLE]

where $P_{\,\mathrm{sat}}\,(\,T^{\,\circ}\,)$ is the saturation pressure at $T^{\,\circ}\ =\ 25\ [\,\mathsf{{}^{\,\circ}C}\,]\,$. These conditions correspond to variations around the relative humidity $0.5$ to the dry and almost-saturated states.

The solution to problem (4) is computed using the RC approaches for r\ =\ \bigl{\{}\,2\,,\,3\,,\,100\,\bigr{\}}, the standard finite-differences one and the spectral one with $N\ =\ 10$ modes. For the last two methods, the spatial discretisation is $\Delta x\ =\ 10^{\,-3}\ \mathsf{m}\,$. The profile of vapor pressure at $t\ =\ 72\ [\,\mathsf{h}\,]$ is shown in Figure 8 (a,b). The time evolution of the vapor pressure in the middle of the layer is given in Figures 9(a) and 9(b). The standard finite-differences and spectral approaches represent accurately the field evolution. On the contrary, an important discrepancy is noted for the approach using two or three resistances (at the order of $\ \simeq\ 200\ [\,\mathsf{Pa}\,]$). The differences become more important when looking at the moisture flux density, illustrated in Figures  10(a) and 10(b). The flux is underestimated by the RC approach based on a few resistances. If we look at Figures 11(a) and 11(b), the absolute differences of the flux reach $\ \simeq 10^{\,-2}$ for the RC approach. However, as shown in Figure 11(c) and 11(d), for these approaches the relative difference on the flux scales with $\ \simeq\ 100\%\,$ for a field of the same order of magnitude. The global error for each approach is given in Figure 12. The most accurate approach is the spectral one, with an error at the order $\mathcal{O}(\,10^{\,-2}\,)\,$. Due to the nonlinearity of the problem, the error of the RC approach, with two or three resistances, have increased when compared to the previous case. As mentioned in the previous case, the error is more important for the R$2$C and R$3$C approaches since the space derivative of the field is computed with a low accuracy as illustrated in Figure 13.

The accuracy digits of each method are reported in Table 2. The digit accuracy of the RC approach with a few number of resistance is much lower compared to the previous case. This lack of accuracy may also be due to the computation of the solution in its physical dimension. Where the temperature scales with $\mathcal{O}(\,10\,)$ in the previous case, here the vapor pressure scales with $\mathcal{O}(\,10^{\,3}\,)\,$. It may introduce additional computational rounding errors. For the spectral approach, it remains stable with four digits of accuracy in the computed vapor pressure. In terms of computational time, the R$2$C and R$3$C have low computational ratio. However, these methods lack of accuracy to compute the fields, particularly the moisture flux density. The Spectral method has the best efficiency providing an accurate solution at a reasonable computational cost.

A parametric study on the number of resistances have been carried out and the error $\varepsilon_{\,\infty}$ is shown in Figure 14. Compared to the previous case of linear diffusion, more resistances are required to reach a satisfying accuracy. A minimal number of $r\ =\ 20$ and $r\ =\ 90$ resistances are necessary to represent accurately the vapor pressure or moisture flux evolution, respectively.

6 Evaluating the methods efficiencies for a real case study

6.1 Description

The purpose is to evaluate now the efficiencies of the numerical methods considering a more realistic case study of heat transfer. The fields were assessed on a building built in Bayonne, France, at the end of the $19^{\,\mathrm{th}}$ Century. With three basements, a west-oriented wall, located at the first floor in the living room, was monitored. Two calibrated monitoring sensors HOBO TMC–$6$–HA were placed at the surface of the wall, as shown in Figure 15(a). A thermal conductive paste was added on the sensor to reduce the contact resistance. An insulated protection has also been placed to avoid incident radiation on the sensors. The location of the sensors was chosen on a part of the wall where heat transfer can be supposed as uni-dimensional between inside and outside ambient. The data measurements were stored with a period of $1\ [\,\mathsf{h}\,]$ during one year. Interested readers are invited to consult Cantin et al., (2010); Berger et al., (2016) for complementary information.

In Berger et al., (2016), complementary measures were used to estimate the thermal conductivity of the wall. Here, the purpose is to use the surface one year measurements to provide the boundary conditions. As shown in Figure 15(b), there are distinguished daily variations with a raising between $150\ [\,\mathsf{days}\,]$ and $300\ [\,\mathsf{days}\,]\,$, corresponding to the summer season.

The wall is considered as homogeneous using the following thermal properties: $k\ =\ 2.48\ [\,\mathsf{W/m/K}\,]$ for the thermal conductivity, $\rho\cdot c\ =\ 2.8\cdot 10^{\,6}\ [\,\mathsf{J/(m^{\,3}.K)}\,]$ for the volumetric heat capacity. Considering the length of the wall $L\ =\ 0.5\ [\,\mathsf{m}\,]$ and a reference time of $1\ [\,\mathsf{h}\,]$, the Fourier dimensionless number equals $Fo\ =\ 1.3\cdot 10^{\,-2}\,$.

For this real and occupied building the initial condition is not known. In addition, the installation of the sensors may cause perturbations on the thermal behavior of the wall. Thus, the first seven days of measurements are discarded and a linear reconstruction between the measured temperatures at given locations was considered as the initial condition at $t\ =\ 0\,$.

6.2 Results and discussion

The temperature is computed within the wall using the RC with r\ =\ \bigl{\{}\,2\,,\,3\,,\,100\,\bigr{\}} resistances, the standard finite-differences and the spectral approaches. For the last two methods, the spatial discretisation is $\Delta x\ =\ 10^{\,-3}\ \mathsf{m}\,$. For the last two methods, the spatial discretisation is $\Delta x\ =\ 10^{\,-2}\ \mathsf{m}\,$. No reference solution is computed for this case. Thus, as the spectral method was shown as the most precise method in the previous cases, it was used as a reference to compute the absolute temperature difference for the standard finite-differences  and R$2$C approaches, illustrated in Figures 16(a) and 16(c). The standard finite-differences  approach has less discrepancies with the spectral solution than the R$2$C one. The maximum difference for the R$2$C reaches $3.1\ [\,\mathsf{{}^{\,\circ}C}\,]\,$. As shown in Figures 16(b) and 16(d), there is good agreement between the R$100$C, standard finite-differences and spectral solutions. As expected, the RC approach with a sufficient high number of resistances enables to accurately represent the physical phenomenon. The heat flux at the right boundary is shown in Figures 17(b) and 17(d). More important discrepancies are noted between the RC approaches with two or three resistances than for the other solutions. In addition, the absolute difference with the spectral solution reaches $\simeq\ 10^{\,2}\ [\,\mathsf{W/m^{\,2}}\,]$ for the R$2$C approach. For the standard finite-differences solution, the difference remains lower than $10^{\,-1}\ [\,\mathsf{W/m^{\,2}}\,]\,$.

In terms of computational cost, the ratio of CPU time $R_{\,\mathrm{cpu}}$ is $0.17\ \mathsf{s/h}$ for the R$2$C, $0.48\ \mathsf{s/h}$ for the R$3$C, $1.9\ \mathsf{s/h}$ for the R$100$C, $1.7\ \mathsf{s/h}$ by using the finite-difference approach and $1.0\ \mathsf{s/h}$ by the spectral model. The latter presents the best efficiency to compute the solution of the problem.

It can be noted that the daily and monthly averaged temperatures are well represented by all methods, as shown in Figures 18(a) and 18(b). There is a perfect agreement among all solutions. The conduction loads are presented in Figures 19(a) and 19(b). For the analysis on the monthly period, the conduction loads are negative for a ten-month period corresponding to winter time. During the summer, the loads are positive, indicating an inward heat flux density. Globally, all the numerical methods enable to estimate the conduction loads with the same order of accuracy. No important discrepancies are noted. This is due to the definition of the heat loads from Eq. (3). The integration is performed for a monthly period. As noticed in Figure 17(b), the R$2$C approximates well the mean heat flux, compared to the R$100$C approach. Therefore, when integrating during a monthly period, the error on the conduction loads is smoothed. On the contrary, when looking at the daily loads of a winter period in Figure 19(b), the discrepancies among the three methods are higher. Particularly, the R$2$C and R$3$C solutions lack of accuracy.

These results lead to the conclusion that the choice of the methods depends on the characteristic timescale of the selected output. If the issue is to accurately analyze the physical phenomenon at a short time scale, then the methods (i) standard finite-differences, (ii) spectral and (iii) RC with a sufficient enough number of resistances, provide a sufficient accuracy. If the output evolves on a relatively large time scale, the RC approach with a few number of resistances may provide a satisfactory approximation of the physical phenomenon.

7 Synopsis of the discussions

It is of major importance to have efficient numerical methods to investigate physical phenomena associated to heat or moisture diffusion through building envelopes. Building simulation programs detail the governing equations but the efficiency of the numerical methods is rarely discussed. This article proposed to investigate the efficiency of three methods by evaluating their accuracy and calculation speed. The first one is the standard one based on central finite-differences and on the Runge–Kutta explicit approach. The second method is the so-called RC model. This approach was initially used to build analogous electric devices to simulate the heat conduction in the $1940$’s–$1950$’s. It was then implemented in the computer algorithm and is still used for educational, research and regulation purposes in building physics. The last method is the spectral reduced order method. It has been recently proposed for building simulation purposes in Gasparin et al., (2018, 2017).

Three case studies have been presented in Section 4, a case of linear heat diffusion was considered. The solutions were computed using the three approaches and compared to a reference solution. The results highlighted that accuracy is satisfactory for the standard finite-differences, the R$100$C and the spectral approaches. Among the three methods, the spectral is the most accurate with the higher number of digits accuracy. It actually gives a good compromise between accuracy and calculation speed. When using a low number of thermal resistances, the accuracy of the RC approach decreases with an absolute difference of the order of $\mathcal{O}(\,10^{\,-1}\,)\ [\,\mathsf{{}^{\,\circ}C}\,]$ with the reference solution. However, the accuracy is lower when looking at the flux density. Indeed, the R$2$C and R$3$C approximate the flux with an absolute difference of the order $\mathcal{O}(\,10^{\,2}\,)\ [\,\mathsf{W/m^{\,2}}\,]\,$. A parametric study was performed to evaluate the minimal number of thermal resistances to reach a satisfactory precision. For this case study, the minimal number is around $10$ and $30$ resistances for the temperature and the heat flux, respectively.

A similar numerical investigation was then accomplished for a case of nonlinear moisture diffusion, in Section 5. The results are very similar and highlight the spectral method is again the most accurate approach with a reasonable computational run time. The speed of computation of the R$2$C and R$3$C models is very low. However, the speed of computation decreases with the accuracy of the solution. It was noted that the number of significant digits of the R$2$C and R$3$C solutions was much lower than for the previous linear case study due mainly to two reasons. First, this case considered a non-linear problem inducing more errors for coarse approaches. Then, the RC approach computes the solution directly in its physical dimension where the others consider a dimensionless formulation of the problem. The vapor pressure scales with $\mathcal{O}(\,10^{\,3}\,)$ where the temperature with $\mathcal{O}(\,10\,)\,$. Therefore, additional computational rounding errors may be introduced in the second case study since the errors of the floating point arithmetic of the computers are minimal for quantities of the magnitude $\mathcal{O}(\,1\,)\,$. The minimal number of resistances to reach a sufficient accuracy with the RC approach was $10$ and $90$ for the temperature and the heat flux outputs of interest, respectively.

The third case study was more realistic case with real measured temperature as boundary condition for the study of heat transfer in a $19$ Century building wall. The simulation was performed for $1\ [\,\mathsf{year}\,]$. No reference solution was available for this case. The standard finite-differences, R$100$C and spectral methods had similar tendencies to represent the temperature evolution in the wall. The R$2$C approach showed a maximum absolute difference of $3.1\ [\,\mathsf{{}^{\,\circ}C}\,]$ with the spectral solution. For the heat flux, the absolute differences were of the order of $\mathcal{O}(\,10^{\,2}\ [\,\mathsf{W/m^{\,2}}\,]\,)$ between both approaches. Within the context of energy efficiency, a relevant output is the time-integrated conduction loads. The latter were calculated for both monthly and daily periods. All the methods have very similar tendencies to estimate the monthly conduction loads. Nevertheless, the R$2$C and R$3$C approaches reveal significant discrepancies with the others when looking at the daily loads. It is due to the integration periods and the fact that the RC approaches with a few temperature approximates well the time mean of the heat flux.

8 Conclusion

With the issue of elaborating reliable models to predict some physical phenomena involved in building energy efficiency, the selection of the numerical methods has to be done considering the characteristic time scale of the aimed output, among other criteria. For the analysis of physical phenomena in a short time scale it is better to use the methods (i) standard finite differences, (ii) spectral or (iii) RC with a sufficient number of resistances. The RC model with a few number of resistances decreases the order of the model. With three resistances for instance, the number of degrees of freedom is only two, reducing consequently the computational effort to solve the problem. However, the fidelity of the model to represent the physical phenomena of heat or moisture transfer is strongly impacted and lacks of accuracy. These results are corroborated with the ones presented in Kircher and Zhang, (2015). With the issue of reducing the computational effort, the spectral method is the most efficient among the three approaches. Indeed, a model with less than $10$ degrees of freedom enables to compute the most accurate solution among the other approaches. Another important point concerns the computation of the surface density of fluxes at the boundaries. The numerical approximation to compute the fluxes introduces errors so that the accuracy is generally decreased compared to the one of the governing fields (temperature or vapor pressure). Since the fluxes are of capital importance on the design of energy efficient buildings, the accuracy of the numerical method to compute them should always been carefully verified.

All these results were obtained on rather a simple case study. The physical phenomena are more complex, involving coupled heat and mass transfer and Robin–type boundary conditions with radiation heat transfer and wind-driven rain. Therefore, before analyzing the physical phenomena in building physics, it is of capital importance to compare the accuracy of the numerical model, including the mathematical formulation of the physical problem together with the numerical methods and the discretisation mesh, with a reference solution. This procedure has to be carried out keeping in mind the output of interest.

Nomenclature

[TABLE]

[TABLE]

Acknowledgments

The authors acknowledge the Brazilian Agencies CAPES of the Ministry of Education, the CNPQ of the Ministry of Science, Technology and Innovation, for the financial support for the project CAPES-COFECUB Ref. 774/2013 as well as the support of CNRS/INSIS (Cellule énergie) under the program “Projets Exploratoires — 2017”. The authors also acknowledge the Junior Chair Research program “Building performance assessment, evaluation and enhancement” from the University of Savoie Mont Blanc in collaboration with The French Atomic and Alternative Energy Center (CEA) and Scientific and Technical Center for Buildings (CSTB).