## 1 Introduction

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

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_2003, besides promoting severe disorders when reaching high levels berger_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_2017. 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_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_2001, ESP-r Clarke_2013, BSim Rode_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_2005; steeman_2009, the explicit Euler in tariku_2010; kalagasidis_2007 or the Crank–Nicolson in vangenuchten_1982. The second method is the so-called RC approach. This method was first used during the second world war Kirkpatrick_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_1953; Robertson_1958with 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_2002; Roels_2017; Naveros_2015. The third method is more advanced spectral method which was recently applied for the solution of diffusion problems through porous building elements Gasparin_2017a; Gasparin_2017b.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 composed of a single material in which the thermal conductivity is denoted as , the density as and the specific heat as . The problem can be formulated by the Fourier (or heat) equation, for and :

(1) |

where is the temperature within the wall at the distance and time .

For the sake of simplicity^{1}^{1}1Different 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:

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

One of the interesting outputs in the building physics framework is the heat flux density at , defined as:

(2) |

Particularly, we denote as the heat flux density computed at :

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:

(3) |

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 , with a single material of permeability and moisture capacity , both depending on the vapor pressure. The formulation of the problem, for and , yields to:

(4) |

where is the vapor pressure within the wall.

The boundary conditions at the extremity of the wall are:

A uniform vapor pressure is assumed as initial condition:

The vapor flux density is similarly computed according to:

### 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_1979. Moreover, the floating point numbers have the highest density within the interval 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 ) function in the Matlab environment. For a numerical model written using dimensionless equations, the accuracy of the floating points scales with . The potential commonly used in the physical model of heat and mass transfer is usually the temperature in , the vapor pressure or the capillary pressure both pressures in . 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 and , respectively. The range of the capillary pressure is between and . Therefore, the accuracy of the floating point can loose up to orders compared to a dimensionless numerical model. So, it is always better to manipulate numerically the quantities of the order of 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_2007, we define the following dimensionless quantities for the temperature:

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

where is a characteristic time. The Fourier dimensionless number is defined, characterizing the importance of the heat transfer through the material:

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

(5) |

together with the boundary conditions:

and the initial condition:

#### 2.3.2 Mass transfer

In a very similar way, we define the following dimensionless quantities related to the vapor pressure field Luikov_1966:

The unitless formulation of the time and space domains are:

The vapor pressure dependent moisture properties are transformed according to:

where and are reference property values. The Fourier dimensionless number is defined, characterizing here the importance of the moisture transfer through the material:

The Fourier number quantifies the first order of the diffusion transfer, while the dimensionless parameters and define the distortion or the nonlinearity of the phenomenon.

The unitless system of differential equation for moisture transfer is:

(6) |

with the boundary conditions:

and the initial condition:

## 3 Numerical methods

In order to describe the numerical schemes, let’s first consider a uniform discretisation for simplicity of the interval :

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

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:

(7) |

whose starting value is directly obtained from the initial condition:

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 for the cell illustrated in Figure 2(a):

(8) |

The average temperature of the cell is defined as

Thus, Eq. (8) becomes:

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

where is the thermal resistance defined as:

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

(9) |

being the thermal capacity defined as:

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 of resistances. Usually, such a model is denoted as with of the order of the unity . This hypothesis corresponds to compute the temperature using 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_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_1928:

(10) |

The algorithm is implemented in the Matlab™ environment. Interested readers are invited to consult Davies_2004; Fraisse_2002 for more details on this approach and Deconinck_2016; Reynders_2014; Jimenez_2009 for examples of applications in building physics.

*(a)*stencil and

*(b)*electrical analogy of the heat conduction transfer.

### 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 from Eq. (5) can be approximately represented as a finite sum Mendes_2017:

(11) |

Here, 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 approximation norm Gautschi_2004. The functions are the corresponding time-dependent spectral coefficients. The parameter

represents the number of degrees of freedom of the solution, also denoted as the order of the solution with

. The main advantage of the spectral method is that , where 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_2017a; Gasparin_2017b.The derivatives are written as follows:

(12a) | ||||

(12b) | ||||

(12c) |

where the dot denotes 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:

with,

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

(13) |

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

which is realised via the the Tau–Galerkin method, which requires Eq. (13) to be orthogonal to the Chebyshev basis functions . Here, the the scalar product is defined by :

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

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

. For linear problems, the system of ordinary differential equations is explicitly built:where , with constant coefficients,

is a vector resulting from boundary conditions and

is the vector of initial coefficients. Initial values of the coefficients are calculated by the Galerkin projection of the initial condition canuto_spectral_2006:(14) |

where , is the dimensionless initial condition. Interested readers may refer to Gasparin_2017a; Gasparin_2017b 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:

where

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

where is the vapor resistance and being the moisture capacity, both defined respectively as:

As described in Section 3.3, the unknown 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:

(16) |

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

where, is a diagonal and singular matrix containing the coefficients of the Chebyshev weighted orthogonal system, is a vector containing the boundary conditions and, 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 using Matlab™ function ode45 Shampine_1997 with an absolute and relative tolerances set to . 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 is computed using a numerical pseudo–spectral approach obtained with the Matlab™ open source toolbox Chebfun Chebfun_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 by the following formula:

where is the number of temporal steps. The global uniform error is given by the maximum value of :

The significant correct digits of the solution are evaluated according to Soderling_2003:

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 CPU and GB of RAM. We define the ratio:

where is the CPU time and 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:

The boundary conditions are defined as sinusoidal variations:

The temperature is computed using three RC model approaches with the thermal resistances . In addition, the problem is solved using the standard finite-differences method and spectral approach with modes. A spatial discretisation step 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 RC, 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 for the approach with two resistances. For the spectral approach, the difference is of the order .

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 as highlighted in Figure 5(a). For the RC model with 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 , 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 with only modes. For the flux, the best accuracy reaches , 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

If we consider that the temperature is computed with a numerical perturbation , with , the approximation of the heat flux is:

Since the perturbations are uncorrelated, and the approximation of the heat flux becomes

Since , 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 for the standard finite-difference approach. With the spatial discretisation , the term is of the order . 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 resistances computes the field with more than twice digits accuracy than the RC one. The spectral method presents the highest accuracy. For the computational time, as expected, the numerical models with and resistances are the fastest approach. Indeed they require a very few computations at each time step ( and , 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 resistances. A good compromise between speed of computation and accuracy is the Spectral approach.

*(a,b)*Temperature evolution at .

*(c,d)*Temperature difference with respect to the reference solution.

*(a,b)*Evolution of the heat flux density at the right boundary.

*(c,d)*Heat flux density difference to the reference solution.

*(a)*and on the heat flux density

*(b)*.

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 , as shown in Figure 7. If the field of interest is the temperature, it can be noted that a number of resistances is sufficient. The error remains stable at 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 to reach an error lower than .

Numerical Model | for | for | ||
---|---|---|---|---|

RC | ||||

RC | ||||

RC | ||||

FDM | ||||

Spectral |

## 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 and the simulation horizon to . The material properties are inspired from Bednar_2005 with a constant moisture capacity and a vapor pressure dependent moisture permeability:

Sinusoidal variations are imposed as boundary conditions:

where is the saturation pressure at . These conditions correspond to variations around the relative humidity to the dry and almost-saturated states.

The solution to problem (4) is computed using the RC approaches for , the standard finite-differences one and the spectral one with modes. For the last two methods, the spatial discretisation is . The profile of vapor pressure at 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 ). 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 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 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 . 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 RC and RC 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 in the previous case, here the vapor pressure scales with . 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 RC and RC 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 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 and resistances are necessary to represent accurately the vapor pressure or moisture flux evolution, respectively.

*(a,b)*Vapor pressure evolution at .

*(c,d)*Vapor pressure difference to the reference solution.

*(a,b)*Evolution of the moisture flux density at the right boundary.

*(c,d)*absolute and

*(e,f)*relative differences of the moisture flux density to the reference solution

Numerical Model | for | for | ||
---|---|---|---|---|

RC | ||||

RC | ||||

RC | ||||

FDM | ||||

Spectral |

## 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 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––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 during one year. Interested readers are invited to consult Cantin_2010; Berger_2016 for complementary information.

In Berger_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 and , corresponding to the summer season.

The wall is considered as homogeneous using the following thermal properties: for the thermal conductivity, for the volumetric heat capacity. Considering the length of the wall and a reference time of , the Fourier dimensionless number equals .

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 .

*(a)*Illustration of the wall and the sensor positions.

*(b)*Variation of the temperature boundary conditions.

### 6.2 Results and discussion

The temperature is computed within the wall using the RC with

Comments

There are no comments yet.