Stochastic Finite Element Method for Modelling Heat Transfer in the Building Envelope

Improving the energy efficiency of the buildings is one of the most effective and fastest ways of reduction of the carbon dioxide emission. However, in the assessment of the energy demand of the buildings, numerous factors are uncertain, i.e. layer thickness, material parameters, climatic conditions, etc. In the present study, mathematical model was developed for analyzing temperature distribution and heat flux in the wall with thermal conductivity of insulation as a random parameter. The results obtained employing a stochastic perturbation technique were compared against the results of the Monte Carlo simulation. Stochastic perturbation technique has been implemented using the tenth order Taylor series expansion. The direct differential method was used to determine the values of Taylor’s coefficient. The obtained results indicate good accordance of the stochastic perturbation technique with the Monte Carlo method. Afterwards, the expected value of the heat flux and its variance were studied for the reference year for a city in the Central Europe. Two cases of the external wall were investigated, in which the thermal insulation was localized either on the internal or the external side of the wall. Performed analyses serve as a good method for assessing the reliability of results obtained using standard, deterministic approach.


Introduction
In the calculations of the heat transfer in the building envelope, assumption of the strictly determined material parameters is usually made. However, both the material properties and initial or boundary conditions can be uncertain. Therefore, also the heat flux should be analysed concerning uncertainty propagation in the calculations and its possible impact on the results.
One of the properties which may differ from the values assumed in the mathematical models is the material thermal conductivity. Usually, its value is taken according the values defined in the standards or in the producers' declarations. Domínguez-Munoz et al. [1] investigated possible disperse of this parameter for, i.a., glass wool, rock wool, expanded polystyrene, extruded polystyrene, polyurethane, and reported high uncertainty factor for each material. Therefore, assumption of the thermal conductivity based only on the insulation type may result in very imprecise results. For the calculations which are performed in order to ensure that the heat losses are smaller than certain threshold value, assumption of upper limit of thermal conductivity for a given material seems acceptable. However, if the case is investigated in order to determine profitability of a given thermo-modernization, it may result in misleading conclusions. Also, the values obtained in the calculations may differ strongly from the actual energy demand in the buildings, what is often reported.
The values of the thermal conductivity declared by the producers may also diverge from the ones measured in reality. Except for influence of the material aging, temperature and relative humidity, which should be accounted for according to appropriate standard equations, uncertainty may arise from the measurement equipment or aleatory uncertainty, related to natural variations. Baldinelli et al. [2] reported variations of the thermal conductivity between different laboratories, in which aerogel, vacuum insulation panel, polystyrene and birch wood fibre insulation boards have been studied.
In majority of the works concerning propagation of the uncertainty in the calculations of the energy demand, Monte Carlo method is usually applied. In this method, one has to determine the probability distribution of the uncertain parameters, perform sampling procedure and execute multiple runs of the mathematical models in order to determine distribution of the results. Even though this method is easy to implement, it is usually time-consuming to perform numerous calculations in order to obtain satisfying convergence, even after application of more complex sampling methods, such as Latin hyper-cube sampling. Influence of the material properties uncertainty on the energy demand has been studied by means of this method by, i.a., Belazi et al. [3], Dolado et al. [4], Prada et al. [5] or Fagianelli et al. [6].
The second method applied in this article, namely perturbation method, belongs to non-sampling methods. It has been applied in the heat transfer calculations by, i.a., Hien and Kleiber [7,8], Kamiński and Hien [9]. In this work, the results of the temperature variance obtained by means of the perturbation method are compared against the ones received using the Monte Carlo method. Furthermore, expected value and variance of the heat flux are investigated for the typical reference year defined for a city in the Central Europe.
Normal distribution of the EPS thermal conductivity is usually assumed [10,11,12]. However, the standard deviation differs strongly between the analyses.

Mathematical model
The transient heat flow, with an assumption of density, specific heat capacity and thermal conductivity independent on temperature, may be written for the dry material as: where is density, is specific heat capacity, is thermal conductivity and is heat flux. Initial conditions may be given as: while boundary conditions can be defined as temperature on the surface, surface heat flux or ambient temperature with the heat transfer coefficient: , , , , , where is the prescribed heat flux, the surface temperature, the convection coefficient, the unit outward normal vector, is the ambient temperature.
Taking into account the Green's theorem the governing equation may be written in the compact matrix form as: , where is capacity matrix, the conductivity matrix and the loading vector defined as: .
(5c) The time derivative is discretized using finite difference scheme: , where is time step.
The governing equation for a given time step may be obtained in the following form: For the uncertain thermal conductivity, two first probabilistic moments of the random field b(x), expected value and variance, may be defined as: where p(b) is probabilistic density function of the random field.
In the nth order perturbation method, the state functions and variables functions are expanded using Taylor series with small perturbation parameter ε > 0 according to equation: (8c) with (9) for k = 1,..,10.
Substituting eq. (8c) into (8a), following equation for the expected value of the state function may be obtained: (10) where is the kth order central probabilistic moment of the random variable b, defined as: . Analogically, following equation can be derived for the variance:

Results
The model has been applied to investigate propagation of thermal conductivity uncertainty in 1D heat flow in external wall composed of 2 layers, brick and expanded polystyrene (EPS), thermal properties of which are defined in Table 1. Thermal conductivity of EPS was assumed as normally distributed with the parameters N(0.037, 0.0056). Due to small resistance of plaster, their impact has been neglected in the calculations.
In order to validate the model, comparison of the numerical results, obtained by means of the Monte Carlo and the stochastic finite element method, has been performed. Two separate cases, in which the insulation has been located either on the external or the internal side of the wall, have been studied in order to investigate influence of heat capacity of the constructive material. The temperature inside the room has been assumed as equal to 20°C. The external temperature has been determined according to the climatic reference year for Poland [13] -see Fig. 1. One can notice that the temperature varies strongly during the year. During winter, the temperature drops below -10°C°, while for summer months it reaches values higher than 30°C. High variation of the temperature in short periods of time can be seen. High gradient between the internal and external environment can be observed, especially during the winter months. The heat transfer coefficients have been assumed as equal to hsi=7.69 W/(m 2 K) and hse=25 W/(m 2 K) [14]. In the Monte Carlo uncertainty analysis, 10000 runs were performed for each case, for which no significant change of the results has been observed by further increasing the amount of runs. The results obtained using two methods after 0h, 10h and 20h of heat flux are presented in Fig. 2 and Fig. 3 for external and internal insulation, accordingly. Comparison of the results indicates good accordance of these methods. It can be noticed that variance of the temperature on the internal side of the wall is high, regardless localization of the insulation, what should result in high variance of the heat flux on the internal surface of the building envelope.    The comparison of the expected value of perturbation stochastic finite element method against the results of the Monte Carlo simulations are presented in Fig. 4 and Fig. 5 for external and internal insulation, respectively. Again, good resemblance of the those two methods has been obtained.  After validation of the SFEM, the heat flux on the internal surface of the wall has been investigated. The results of the heat flux are presented in Fig. 6 and Fig. 7 for the entire year and January, accordingly. It can be noticed that similar results have been obtained for the analyses performed for insulation on the internal and external side of the structural material. For external insulation, slight delay of the heat flux and heat flux variance can be observed. Obviously, heat flux depends on the temperature gradient between the internal and external side of the wall. Thus, the highest values (up to 9 W/m 2 ) can be observed for winter months, in which the temperature may decrease in Poland below 0°C.  The heat flux variance is presented in Fig. 8 and Fig. 9 for entire year and January, respectively. Evolution of the heat flux variance strictly corresponds to evolution of the heat flux and takes the highest values for the winter months, in which the temperature drops below 0°C, and drops to approximately 0 (W/m 2 ) 2 for summer months, in which the temperature fluctuates around 20°C, which is equal to the internal temperature. The peaks of the heat flux variance are slightly higher for the internal insulation. For the analyses performed for external insulation, due to heat capacity of brick, the heat flux variance is slightly delayed with respect to the temperature variation.    The relative value of variation coefficient, which has been calculated as the ration of the heat flux variation coefficient to the variation coefficient of thermal conductivity, is presented in Fig. 10 and Fig. 11 for the entire year and January, accordingly. This value fluctuates around 0.85 for the winter months, in which high temperature gradient between the wall surfaces occurs. For the summer months, during which the temperature gradient may drop to 0°, the relative value of variation coefficient changes rapidly. The probability density function of the annual heat flux density is presented in Fig. 12. One can notice similar distributions for the insulation placed on internal and external surface.

Conclusions
The influence of the thermal conductivity uncertainty on the heat transfer in the external wall during the reference year in a city in the Central Europe has been presented. The results of the temperature distribution and expected value of temperature in the wall, with internal and external insulation, indicate good accordance of the Monte Carlo method and the stochastic finite element method.
After validation, the model has been applied to study variance and expected value of the heat flux. For the cases with external insulation, evolution of heat flux and its variance is slightly delayed when compared against the results obtained for internal distribution. The variance changes can be strictly related to the heat flux variation. The heat flux variance is higher for the internal insulation.
The model can be applied to study propagation of the uncertainty of other material properties, such as density, boundary or initial conditions.