1. Introduction
Energy demand of the buildings contributes significantly to the total global energy consumption and the emission of carbon dioxide. Due to the current actions of the human impact reduction on the environment, determination of the heat losses or the energy demand is an important step in the building designing process. Nowadays, advanced numerical simulation is commonly used in order to obtain accurate results. In such calculations, it is usually assumed that material properties and boundary conditions are strictly determined. However, numerous model parameters used for these calculations are uncertain. Those uncertainties can refer to geometry and construction (layer thickness, simplification of the geometry), material parameters (thermal conductivity, density, etc.), environmental data (temperature, wind speed, air flux, etc.) or model uncertainties (numerical errors, simplifications of the physical phenomena, etc.). Fadale and Emery [
1] distinguished 3 types of uncertainties for the material systems: intrinsic uncertainties, related to heterogeneity of the material or boundary conditions, information uncertainty, related to lack of data, and experimental uncertainty, related to the experimental errors. They may be also referred to as epistemic uncertainty, related to missing data, and aleatory uncertainty, caused by natural variations of certain parameters. Combination of the randomness of these parameters induces also considerable fluctuations of the total response function. Therefore, the analyses should be accompanied by the investigation of such uncertainties’ propagation.
The results of the laboratory measurement should contain information about the uncertainty of the measurement itself, i.e., EN 12524 [
2] declares 3% of uncertainty for the guarded hot plate and the heat flow meter hot plates, which are frequently used to determine thermal conductivity of the material at dry state. For the thermal conductivity tests, uncertainty for a wide range of materials has been defined, i.e., [
3], indicating that the uncertainty value depends on the thermal conductivity, with higher uncertainty for smaller conductivity values. Baldinelli et al. [
4] performed comparison of the measured thermal conductivity for aerogel, polystyrene, VIP and wood chips between different laboratories. The relative repeatability standard deviation and relative reproducibility standard deviation for polystyrene were equal to 0.9% and 2.0%, accordingly.
Thermal conductivity may vary strongly with reference to the value declared by the producer not only due to the experimental error and material heterogeneities. For majority of the materials, its value for a given production series remains approximately constant and only small variation may be usually observed. However, depending on the material delivery, workmanship errors or conditions of production and storage, differences in material properties may be observed for various production series. Furthermore, the values measured during the laboratory tests and the ones declared by the producer may differ due to the safety margin imposed. Additionally, thermal properties of the material may be changed due to the external loads present during their lifetime. Exact information about thermal conductivity may be missing at the design stage, what causes epistemic uncertainty. In [
5], a range of thermal conductivity for specific insulating materials has been defined in order to assess possible range of this parameter for a given material. Concluding, the real heat flux may diverge in reality from the results of deterministic calculations, which are based on the declared values of material parameters.
The typical reference year, which is usually assumed in the energy-related analyses to describe the external environment conditions, is based on the long-term observations of numerous parameters in the given location. However, due to the fluctuations in the weather, extension of the observation period taken into account during determination of the typical reference year and climatic changes, these values may also be considered as uncertain. Therefore, due to scatter of the results, such a set of data should be considered with attention and propagation of uncertainty of climatic conditions is essential in order to ensure thorough analyses of results.
Numerous researchers focused on the analyses of the impact of parameter uncertainty in the calculations of the heat transfer problems and the energy demand in the buildings. In majority of them, the Monte Carlo method or the quasi-Monte Carlo method has been applied. Randomness of parameters in the heat conduction [
6,
7,
8] and radiative transfer equation [
9,
10] have been taken into consideration by means of such approaches. Recently, non-sampling techniques are also often applied in the stochastic analyses. The polynomial chaos algorithm has been applied by Xiu and Karniadakis [
11] in the calculations of the transient heat conduction. The second-order perturbation method, which may be applied to the problems with relatively low coefficient of variation, has been used for the linear and nonlinear heat transfer problems [
12,
13,
14]. The generalized stochastic perturbation method, with a higher order of the Taylor series expansion, has been applied, i.e., by Kamiński and Strąkowski in the analyses of the steel tower structure exposed to fire [
15].
In this article, impact of the thermal conductivity and the ambient temperature uncertainties on the heat losses through the external wall have been studied. Distribution of thermal conductivity in the expanded polystyrene (EPS) is determined according to the available laboratory data. Two types of distributions are determined: when exact type of EPS is not defined (information uncertainty [
1], usually occurring at the early stages of the design process, which induces high range of possible material parameters), or when precise parameters are already assumed (intrinsic uncertainty [
1]). Furthermore, we analysed not only randomness of the material parameters, but also uncertainty of ambient temperature and took it into account by considering the Robin’s boundary conditions. The generalized perturbation method is presented.
The model is applied for the analyses of 1D heat flux through the two-layer external wall. Three cases have been considered, in which thermal conductivity or external temperature are assumed as uncertain. Expected value and variation of the heat flux on the internal side of the two-layers wall is calculated. The results are verified against data obtained using the Monte Carlo method for 10,000 individuals. Additionally, the Monte Carlo analysis has been repeated three times for 100, 1000 and 10,000 individuals in order to investigate possible scatter of results in the Monte Carlo method.
To the authors’ best knowledge, despite application of the perturbation method to the heat transfer problems in the past [
14], it has not been applied to analyses of the heat flux in the walls in realistic environmental conditions for the walls, based on the typical meteorological year, with the boundary conditions of the third type. This work should prove usefulness of the generalized perturbation method in such calculations and provide additional insight into propagation of uncertainties.
Determination of the distributions of thermal conductivity and ambient temperature is discussed in
Section 2. Mathematical model of the heat transfer and perturbation method is presented in
Section 3. Case study, in which propagation of the determined uncertainties in the heat flux calculations is considered, is presented and discussed in
Section 4.
2. Thermal Conductivity Uncertainty
In the stochastic analysis, thermal conductivity is usually assumed to be normally distributed [
16,
17,
18,
19,
20,
21]. Its value is influenced by numerous factors. Depending on the given material micro-structure, impact of three basic types of heat transfer (conduction, convection and radiation) can differ strongly. In the insulating materials, content of voids is very large. They may contain air or other gases introduced in order to decrease thermal conductivity. For materials with small density, thermal conductivity decreases with the density increase due to the long-wave radiation inside the pores. For heavier materials, thermal conductivity is directly proportional to density due to larger content of the solid matrix and its conductivity. In the materials with higher density, above approximately 50 kg/m
3, effect of radiation inside the pores is negligible [
5]. This effect may cause different distributions of thermal conductivities in materials of different density.
The time of aging, which is accounted for in the calculations of thermal conductivity, should be not smaller than half of the building’s lifetime, which is usually assumed as equal to 50 years. Effect of aging is of paramount importance to the materials in which foaming is performed with gases other than air, particularly extruded polystyrene and PIR. Interesting results concerning increase of thermal conductivity in such materials have been presented, i.e., by Katarzyna Firkowicz-Pogorzelska [
22], in which 22.2% and 31.1% increase of thermal conductivity has been observed for XPS and PIR, accordingly. Expanded polystyrene, which is considered in the presented work, has been investigated by Lee et al. [
23]. While resistance of XPS has decreases significantly, hardly any changes of thermal conductivity of EPS could have been noticed. Therefore, aging effect has not been considered in the presented work.
One of the most important parameters concerning thermal properties of the EPS is its moisture content. It has been vastly studied in the literature. Andreotti et al. [
24] proposed an in situ hot box which enables to determine hygrothermal performance of the building component and compared it with numerical simulations [
25] for analyses of energy refurbishment in historical buildings. Impact of moisture on thermal conductivity for EPS and XPS has been studied [
26,
27]. Even though for high moisture content resistance of the insulation decreases strongly, for low moisture content, below 0.3% by volume, only small increase of the EPS thermal conductivity of approximately 2% can be expected. In case of the walls with EPS on the external surface and appropriate coating, only small variations of relative humidity can be expected in this material and low relative humidity is usually maintained. Additionally, according to [
28] effect of moisture is independent. Therefore, it has not been considered in this research.
The database of thermal conductivity tests, performed in years 2016–2019, of insulating materials, provided by the Construction Control Authority in Poland [
6], has been used. The data obtained for expanded and extruded polystyrene has been considered. 1192 data points, with 4 values given for each material tested, have been used in the analyses.
Two histograms have been elaborated: in the first one, measured thermal conductivity is presented (
Figure 1a), while in the latter one, relative thermal conductivity (
λr), defined as relation of the measured thermal conductivity (
λm) to the value declared by the material producer (
λd),
λr =
λm/
λd, is given (
Figure 1b). For each histogram, parameters of the normal and lognormal distributions have been determined and the Pearson’s chi-square statistical test has been performed, indicating better fit for the thermal conductivity for normal distribution of parameters (0.0368, 0.0045
2) and for the relative thermal conductivity for lognormal distribution of parameters (0.0146, 0.0553
2). The thermal conductivity of the materials is in general agreement with the data presented in [
5] and covers a wide range of thermal conductivity values. The relative thermal conductivity indicated asymmetric shape, what contradicts assumption of normal distribution, often considered in such analyses. However, it should be noted that this curve is valid only for the case in which precise thermal conductivity has been defined based on the declaration of the producer. The shift in the graph is probably related to the fact that producers usually declare thermal conductivity slightly higher than the values observed in the laboratory as a safety margin. Finally, normal distribution of thermal conductivity will be assumed for the simulations with information uncertainty and lognormal distribution for the case with intrinsic uncertainty.
3. Mathematical Model
The transient heat flow can be given as:
where
is density,
—specific heat capacity,
—thermal conductivity and
—heat flux. Assuming that density, specific heat capacity and thermal conductivity are independent on temperature, one can write:
The initial conditions are written as:
In the article, boundary conditions of the third type, in which ambient temperature and heat transfer coefficients are defined, are applied:
where
α is the convection coefficient and
is ambient temperature.
The governing equation can be written as:
where
C is the capacity matrix, defined as:
K is the conductivity matrix:
and
Q is the loading vector:
The time derivative is discretized using truncated Taylor’s series expansion (finite difference scheme).
Denoting the random field by
b(
x) and its probability density function by
p(
b(
x)), first two probabilistic moments of the random variable can be written as:
In the perturbation method, variables and state functions,
f(
b(
x)), are interpolated using the Taylor expansion using small perturbation parameter
ε > 0 as:
where
. Value of
ε is usually assumed as equal to 1.
Substituting (11) into (9) one can obtain, for unsymmetrical distributions:
where
is the central probabilistic moment of the random variable
b of the
k-th order, defined as:
Analogically, substituting (11) into (10) one can define variance as:
4. Numerical Simulation
Numerical analyses have been performed for 1D heat flux through the two-layer external wall. Thickness of the insulation has been assumed in such a way, that maximum thermal transmittance for Poland for 2021 is fulfilled. The constructive layer is made of brick, while thermal insulation of expanded polystyrene (EPS). The material parameters are given in
Table 1. The scheme of the considered geometry and boundary conditions is presented in
Figure 2.
We will analyse how the uncertainty of thermal conductivity and ambient temperature impacts the probability distribution of temperature in a few points through the wall and the heat flux on the internal surface, which might be identified as the heat loss through the wall. As a result, we will obtain the expected value and variance of heat flux (on internal surface). The obtained statistical parameters allow the user to calculate the probability of any considered event, i.e., quantiles of heat flux.
In total, three cases have been analysed: Case 1, with normal distribution of thermal conductivity, assumed for lack of knowledge concerning exact declared value of thermal conductivity, Case 2, in which lognormal distribution of thermal conductivity around its declared value has been considered and Case 3, in which uncertainty of external temperature has been investigated.
The internal and external plasters were not taken into account due to their small thermal resistance and heat capacity. This simplification leads to relatively small errors, which might be neglected during the stochastic analysis. The heat transfer coefficients on the external and internal wall were taken from ISO 6946 [
29] as equal to 25 W/(m
2K) and 7.69 W/(m
2K), appropriately. All calculations were performed assuming constant internal temperature, equal to 20 °C. Ambient temperature changes according to the climatic data taken for the climatic reference year for Poland [
30]. For clarity of the results, the calculations have been performed for January only.
In the analyses, five points in the wall have been considered: P1 (internal surface), P2 (center of the brick), P3 (between EPS and brick), P4 (centre of EPS), P5 (external surface). In the analyses of impact of number of the Monte Carlo runs on the scatter of the results, simulations for 100, 1000 and 10,000 simulations have been repeated three times. In other analyses, results obtained for 10,000 individuals have been used.
4.1. Case 1: Uncertain Thermal Conductivity–Information Uncertainty
For the early design stage, exact type, producer and thus also thermal conductivity of the insulating materials is often not defined. Thermal conductivity has been taken for such a case as normally distributed with parameters N(0.037, 0.0045
2) (see
Figure 1a).
Comparison of the expected value and variance of the total heat loss through 1 m
2 of wall in January,
Qtot, obtained using the Monte Carlo method for different numbers of the Monte Carlo runs, is presented in
Figure 3. For 100 simulations, differences between three sets of results of appr. 53 Wh for expected value and 7 (kWh)
2 for variance have been observed, what is a high dispersion when compared with the values of appr. 2670 kWh and 95 kWh. Scatter of results for 1000 individuals was only slightly better. It can be noticed that even for 10,000 individuals, for which the calculation time cost was relatively long (approximately 3 h), quite a large scatter of the results has been obtained.
Comparison of the expected value and variance of temperature in the centre of the EPS (P4) and between brick and EPS (P3), obtained by means of the Monte Carlo method and perturbation method, is presented in
Figure 4a,b, while the expected value and variance of the heat flux on the internal side of the wall is presented in
Figure 4c,d. Only these points are presented in the graph for clarity of the results. It can be noticed that high similarity of variance distribution in the time in different points was obtained for the perturbation stochastic finite element method and the Monte Carlo method for all cases, when arbitrary material heat conductivity changes randomly. The validating calculations were performed for 48 h.
After validation, expected value and variance of temperature and heat flux on the internal surface of the envelope were calculated for January by means of the perturbation method. The expected value and variance of temperature in 5 points considered is shown in
Figure 5. Expected values of temperature correspond to the external environment temperature. It can be noticed that the highest variance is observed between the constructive material and the insulation during the entire month. The lowest value of variance is on the external surface due to the fact that the external surface resistance is low and temperature in this point will not vary strongly from the external temperature; thus, variations in the thermal conductivity of the insulation will also not induce high variations of temperature in this point. Similarly, variance of temperature on the internal surface is quite low, below 0.01 (°C)
2. Comparing P2 and P4, one can observe that variance in the brick, of high thermal mass, does not fluctuate strongly during the whole month. Meanwhile, variance in the middle of the EPS is strongly modified accordingly to the external temperature variations.
The expected value and variance of the heat flux on the internal side of the wall is presented in
Figure 5c,d. The expected value of heat flux varies approximately from 3 to 4.6 W/m
2, while its variance is in the range of 0.12 (W/m
2)
2 to 0.28 (W/m
2)
2. Peaks of variance are, obviously, strictly related to the expected value of heat flux and thus also to the external temperature.
4.2. Uncertain Thermal Conductivity–Intrinsic Uncertainty
The distribution of relative thermal conductivity, considered as relation of the actual thermal conductivity to the value declared by its producer, has been assumed as lognormal with parameters (−0.0146, 0.0553
2) (see
Figure 1b). Basic thermal conductivity equal to 0.037 W/mK has been assumed in the calculations. It should be noted that the thermal conductivity of EPS is distributed according to non-symmetric function, unlike in Case 1. Consequently, expected value and variance must be determined considering both even and odd statistic moments.
Expected value and variance of the total heat flux in January observed for different numbers of individuals in the Monte Carlo analysis is presented in
Figure 6. Similarly as in the previous case, when results obtained for sets of 100 or 1000 samples are compared, the results differ significantly even for the expected value of the heat losses–the highest difference of 30 Wh is obtained, what can be considered as a high scatter. For 10,000 individuals, relatively large scatter of the results is observed also for 10,000 Monte Carlo runs for variance of the heat flux.
Expected value and variance of temperature and heat flux as a function of time, determined using the Monte Carlo method and the perturbation method, are presented in
Figure 7a–d. Good accordance of the results for these two methods has been obtained.
Expected value and variance of temperature, determined using the perturbation method, is presented in
Figure 8a,b, accordingly, while of heat flux in
Figure 8c,d, respectively. Similarly to Case 1, the highest variation of temperature may be observed in P3, between EPS and brick. The variations of temperature and heat flux in Case 2 are approximately 5 times smaller than the ones observed in Case 1.
4.3. Uncertain Ambient Temperature
In the last section we analysed how the randomness of ambient temperature affects the statistical distribution of temperature and heat flux. In the analysis of uncertain ambient temperature, distribution of temperature
N(
Tref,year(
t), 3.0
2), with
Tref,year being the temperature according to the reference year data for polish time station [
30], has been assumed.
Comparison of the results obtained by means of the Monte Carlo method for different numbers of simulations is presented in
Figure 9. The results of the perturbation method are presented in
Figure 10, where the results of the Monte Carlo method for 10,000 individuals are also presented for comparison. Once again, very good resemblance between temperature variances calculated using the stochastic perturbation method and the Monte Carlo method can be observed. In all the cases analysed, small differences noticed in the variation of the heat flux may be contributed to scatter of the results obtained in the Monte Carlo method for 10,000 runs. Unlike in the previous cases, temperature variance calculated in arbitrary position do not change along with time, it remains constant in arbitrary nodes during entire analysed period.
Time evolution of the expected value and variance of temperature is presented in
Figure 11a,b, accordingly, while of heat flux in
Figure 11c,d, respectively. One can notice that variance of temperature is constant for given points, similarly as variance of heat flux. Highest variance of temperature is observed for the points on the external surface, while for internal surface it is negligible.
5. Discussion
In the first part of the manuscript, the statistical distribution of the EPS thermal conductivity was investigated. Based on the statistical testing it was concluded that this feature should be considered in different ways depending on the type of uncertainty. At the early design stage of the building project, when exact type of the insulation may be unknown, normal distribution, based on a wide range of thermal conductivities measured, can be assumed. Otherwise, if a value defined by the material producer is assumed, lognormal distribution accounting for randomness of the actual thermal conductivity with reference to the declared value should be considered. This finding was directly applied in the statistical analysis conducted by means of Stochastic Finite Element Method.
The impact of uncertainty of the material properties and boundary conditions on the heat flux and temperature distribution in the building envelope was presented. Two types of forward uncertainty analyses of the input parameters were performed. The Monte Carlo method is widely recognized as a very time-consuming technique due to necessity of performing multiple runs. Whilst giving very accurate results, the perturbation stochastic finite element method is much more efficient than the Monte Carlo method for transient heat transport in double-layer external envelope. Furthermore, the sampling size in the Monte Carlo method influences strongly the results, what has been presented in the manuscript. Therefore, limiting the population may strongly alter the results’ reliability.
The second-order perturbation method is usually applied to problems characterized by low variation coefficients. The generalized order perturbation method can be used for a wider range of problems. It involves a necessity to calculate derivatives of the response functions with relation to the uncertain parameters, what can be obtained by means of the direct differential method (DDM) or the response function (RF). The DDM, which has been used in this work, requires additional modifications in the source code. If the mathematical model is complicated, such implementations would be time consuming. The RF requires performing the simulations for a set of values of uncertain parameter. Even though the calculations have to be performed a few times, this number will be much smaller than number of simulations run in the Monte Carlo method.
After validating, the perturbation method has been applied to analyse propagation of the uncertainties of thermal conductivity and external temperature in the calculations of heat transfer in the building envelope. The largest values of temperature variance can be noted in the zone just in between two materials, when the material property uncertainties were investigated.
The information uncertainty, which is caused by lack of data concerning exact type of insulation material at the early design stage, causes much higher variance of results when compared to the intrinsic uncertainty—in the analysed case study it was approximately 7 times larger.
Since the internal temperature and heat transfer coefficient are assumed to be constant, variance of heat flux on the internal surface depends solely on the temperature variance on the internal surface. By analogy, high variance of temperature on the external side induces high disperse of the results of heat flux on the external side. For random thermal conductivity of EPS, high variance of heat flux on the internal surface and low on the external surface were identified.
The heat transfer uncertainty investigation allows not only to determine insignificant parameters and analyse robustness of the model, but it may determine errors that can be possibly done by assumption of strictly determined parameters. It has been presented that possible variance of the material parameters may cause large discrepancies between the model and reality. Considering the material and boundary conditions randomness, it allows to ensure greater assurance of the results obtained by the given mathematical model. The Monte Carlo method, which is usually applied for such a purpose, does not give precise results when number of the simulations is too small, while for a large number of individuals its calculation time cost is high even for simple problems considered. The perturbation method enables to perform such analyses in a quick manner with robust outcomes. It may allow to obtain probability distribution functions of the heat losses or quantiles of the energy demand, which can be further applied in the calculations of the energy demand or verification if a building fulfils certain national requirements. It also enables determining risk that the energy demand will be higher than certain threshold. The presented statistical analysis can be used as a quick approach to assess uncertainty impact on the buildings’ energy demand and may provide more information than the deterministic approach.