Next Article in Journal
Robustness of Wireless Power Transfer Systems with Parity-Time Symmetry and Asymmetry
Previous Article in Journal
Application of Phase-Shifting Transformer with Longitudinal and Transverse Voltage Regulation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Comparison of the Performance of New and Traditional Numerical Methods for Long-Term Simulations of Heat Transfer in Walls with Thermal Bridges

1
Department of Fluid and Heat Engineering, University of Miskolc, 3515 Miskolc, Hungary
2
Institute of Physics and Electrical Engineering, University of Miskolc, 3515 Miskolc, Hungary
3
Mechanical Power Engineering Department, Al-Baath University, Homs 77, Syria
4
Mechanical Engineering Department, University of Technology—Iraq, Baghdad 10066, Iraq
*
Author to whom correspondence should be addressed.
Submission received: 14 April 2023 / Revised: 3 June 2023 / Accepted: 6 June 2023 / Published: 8 June 2023
(This article belongs to the Section G: Energy and Buildings)

Abstract

:
Several previous experiments showed that the leapfrog–hopscotch and the adapted Dufort–Frankel methods are the most efficient among the explicit and stable numerical methods to solve heat transfer problems in building walls. In this paper, we extensively measure the running times of the most successful methods and compare them to the performance of other available solvers, for example, ANSYS transient thermal analysis and the built-in routines of MATLAB, where three different mesh resolutions are used. We show that the running time of our methods changes linearly with mesh size, unlike in the case of other methods. After that, we make a long-term simulation (one full winter month) of two-dimensional space systems to test the two best versions of the methods. The real-life engineering problem we solve is the examination of thermal bridges with different shapes in buildings to increase energy efficiency.

1. Introduction

Significant progress toward a sustainable economy may be made via improved building energy efficiency, where the buildings account for 40% of the primary energy use and 24% of the generation of greenhouse gases [1]. When compared to other industries, the construction industry has the ability to significantly reduce energy usage and, by extension, greenhouse gas emissions [2]. Heat transfer calculations are widely used in buildings, including determining how much energy is lost or gained through the building’s envelope (heat conduction), conducting environmental analyses inside the building, and troubleshooting issues with specific materials or structural elements [3,4].
Exterior wall conduction, interior mass conduction, heat gain/loss conversion to cooling and heating load, and ground heat loss from the slab-on-grade floor and basement walls are all issues related to conduction heat transfer that arise in the context of buildings, so the rate at which heat moves through walls can be changed by changing the thickness of the walls or by making the layer of insulation thicker [5]. Transitory wall conduction heat transfer reacts to weather conditions such as temperature swings, sunshine, air movement, etc. Increasing the building envelope’s thermal insulation and decreasing the heat loss via walls is one of the most efficient strategies to raise a structure’s energy efficiency and decrease its energy consumption [6]. Walls account for around 35% of the heat loss in a house, with the vast majority of the heat entering and leaving them [7]. Since the walls of a house are in direct contact with the colder air outside, heat is often lost by conduction or physical contact.
EN ISO 13789 [8] is a commonly used standard for calculating transmission heat transfer through a building envelope; which in turn mentions two standards for thermal bridges; EN ISO 10211 [9] provides a structure for thorough calculations of thermal bridges in building construction and EN ISO 14683 [10] provides a simplified approach with default linear thermal transmittance values.
Because thermal bridges have a big effect on heat loss and reduce building energy efficiency, durability, and air quality, requiring integrated thermal and structural design [11], and the thermal bridge models that were looked into should be especially interesting to architects, civil engineers, and people who work in the insulation materials industry. The examination of the effects of thermal bridges is meant to illuminate the feasibility of implementing energy renovation techniques in existing structures [12]. Thermal bridge effects are analyzed using a stationary linear thermal transmittance method for a wide variety of foundation types and construction parameters. Validation tests on over a thousand real-world examples corroborated the validity of those correlations, allowing for a straightforward and practical analysis of thermal bridges in existing structures and providing practitioners with a resource that accounts for the vast majority of scenarios [13].
A thermal analysis may reveal the heat distribution patterns of a system or its individual parts. Most research into thermal quantities focuses on temperature distributions, thermal fluxes, and heat capacities. Since diverse heat transfer applications within engineering fields involve many thermal models, the analysis of transient heat transfer is an essential problem that is often solved using numerical rather than analytical methods. Analytic techniques provide precise results, but they can only be used for isotropic, homogeneous situations with straightforward geometries and boundary conditions [14].
Heat transfer across a layer is found to be proportional to the Δ T temperature difference across the layer and the heat transfer area, but inversely proportional to the layer thickness. It is called “transient conduction” when the mode of thermal energy transfer occurs throughout a time period in which temperatures vary at any location inside an item. Time-dependent temperature fields are referred to as “non-steady-state” conduction. Conduction of heat through the composite wall may be modeled if required. Assuming this is the case, we simply require boundary conditions at the outside wall surface, with the same conditions applied to the inside wall surface. Solving heat conduction equations across the wall layers yields their temperatures.
The heat conduction equation, often known as the “heat equation”, is a partial differential equation (PDE) used to explain the conduction of heat through a solid. Innovative analytical approaches [15] exist for spatially homogeneous systems, and mathematicians create and evaluate most numerical techniques for homogeneous cases. However, some analytical solutions for steady-state and transient situations are also available for the one-dimensional inhomogeneous system, for example, multilayer systems [16]. Most of the time, these systems are used to determine how much heat is gained or lost through the outside envelope and how much heat is stored inside buildings [3]. Density, thermal conductivity, and specific heat are the key examples of material qualities that may change greatly in a building’s system. Weather conditions also change over time. These imply that most heat conduction issues are multi-dimensional and transient [15]. As a result, we need to use numerical computer simulation.
Our research group works on transient heat transfer calculations using fundamental physical laws (ab initio approach). Therefore, it is expected that these results are much more accurate than those based on the usual (ISO) standards, which are steady-state calculations without solving the transient PDE and therefore cannot properly take into account, for example, the heat capacity of the envelope. Our long-term goal is to revolutionize these simulations (at this stage by the numerical methodology) in order to make transient simulations more available due to reduced computational cost and programming difficulty, as follows.
Heat transfer and other similar problems are routinely solved by lots of numerical methods. Most of them fall into one of two categories: explicit method or implicit method.
Compared to one another, each type has some key benefits and at least one serious drawback. The widely used explicit algorithms, such as the classical FTCS (forward time central space), execute time steps in a small amount of time and are quite straightforward to parallelize. The problem is that the so-called Courant-Friedrichs-Lewy (CFL) limit constrains their stability domain and causes them to become unstable if the time step size is above this limit [16]. Moreover, since an insulated wall consists of very different materials, the coefficients can change sharply from mesh point to mesh point; thus, this limit is rather tiny and the stiffness ratio is large.
Since the stability of the implicit algorithms is substantially greater, many scholars see them as superior and frequently employ them to solve these and other similar PDEs [17,18,19]. The drawback comes from the fact that in each time step, a system of algebraic equations must be solved, a procedure that is difficult to parallelize. In addition, the large amounts of the required RAM can also make calculations extremely slow. These problems are remarkable typically when the physical space is two or three-dimensional.
It has been shown that explicit approaches can be more effective even when tiny time step sizes are required [20]. Moreover, various ingenious combinations of explicit and implicit methods, such as semi-explicit or semi-implicit methods, have also been presented [21,22,23,24]. However, they do not really solve the problems and the drawbacks of the explicit and implicit approaches discussed above.
In light of the aforementioned data, it is reasonable to assume that explicit algorithms, particularly if they have improved stability qualities (see, for example, [25,26,27,28,29,30,31,32]), will have a growing comparative advantage over time. We started working on new explicit schemes a couple of years ago for determining heat conduction in any number of spatial dimensions. In our original articles [33,34,35,36,37,38,39,40], the novel algorithms were theoretically and experimentally investigated. They were proven to be stable for the linear diffusion- or diffusion-reaction equation, and their theoretical order of convergence is also stated. The algorithms were tested using analytical as well as numerical reference solutions. It was demonstrated that they deliver fairly accurate results at a remarkably greater speed than the widely used methods, for example, the built-in ODE solvers of MATLAB and that they can be used in a wide variety of situations involving random discontinuous parameters and initial conditions.
In our paper [38], we tested 12 explicit and stable numerical algorithms to simulate heat conduction in building’s walls, with and without insulation, using equidistant and non-equidistant meshes. We obtained that the original odd–even hopscotch (OOEH), the leapfrog–hopscotch (LH), and the Dufort–Frankel (DF) are the most effective methods. Then, in [37] we adopted some of the above-mentioned 14 methods to include not only heat conduction but also convection and radiation. According to the results, the LH and the non-standard version of the OOEH are the most accurate if the system is not really stiff. However, the OOEH becomes less accurate when the stiffness is larger, which is the case if the mesh is non-equidistant and/or there are materials in homogeneity. In these cases, the LH takes the lead, but the DF, as well as the shifted-hopscotch (SH) and asymmetric hopscotch (ASH) methods, also perform well. Those methods which are unconditionally stable for the simple conduction case can be used without stability problems with fairly large time step sizes, so they outperform the conventional explicit time integrators. Since usually the LH was the best method, we devoted a whole paper [37] on the question of how to implement the convection and radiation term in an optimal way. However, in those papers [37,38], no running-time measurements were made, and the performance was evaluated only in terms of accuracy versus time step size.
In our current work, we continue the above-mentioned investigations. We use ANSYS’ thermal analysis solutions that help engineers solve the most complex thermal challenges, and predict how their designs will perform with temperature changes. However, because simulation by this kind of software takes a long time and requires serious computer resources, we compare our methods with ANSYS to investigate runtime, stability, and other features [41]. Now the goal is to systematically evaluate how the performance of the various solvers (including MATLAB routines and ANSYS) depends on the mesh settings to see which one is optimal for certain accuracy requirements. Therefore, the paper is organized as follows: Section 2 presents the examined system with the equations, followed by the numerical methods and implementations of the convection and radiation terms, as well as the methods used for comparison purposes. Section 3.1, Section 3.2, Section 3.3 display the preliminary conditions for the simulation of the wall: materials, mesh construction, initial, and boundary conditions.
In Section 3.4, we start by verifying our codes in a 2D system for three different cases and comparing the results against simple analytical solutions on a grid that is equidistant. Section 4 displays the results of the numerical tests performed to compare differences in errors and running times. Section 5 shows the second part of this study, where we chose our best methods to perform a long-term simulation of one month for a real wall in Miskolc city with and without insulation. Two kinds of the thermal bridges are included, and it is calculated how much energy is lost due to the thermal bridges. Section 6 finally concludes with a summary of our findings.

2. The Studied Cases

2.1. The Equation and Its Discretization

Using a linear parabolic PDE, we can characterize the phenomenon of the simplest Fourier-type heat conduction inside a homogeneous medium with a heat source as follows:
u t = α 2 u + q ,
where u = u r , t is the temperature, q is the heat generation or heat source, and α is the thermal diffusivity, which is given as α = k / ( c ρ ) , where c = c r , t , k = k r , t , ρ = ρ r , t are the specific heat, heat conductivity, and the density of the material, respectively.
Newton’s law of cooling suggests that a term K u a u can describe convective heat transfer [42], where the ambient temperature u a (measured in Kelvin) can be considered as independent from the unknown variable u. Consequently, it makes sense to include the K u a term in the expression q, which represents the heat source. Stefan-Boltzmann’s law [43] says that a term σ u 4 can be used to describe the surface’s radiative heat loss, where now the proportionality constant σ is the product of the Stefan-Boltzmann constant and the surface area, which are all positive quantities. We can incorporate the incoming radiation, which may include direct sunshine, into the source term q similarly to the K u a term above. The following is an extension of the heat conduction Equation (1) to include the terms for convection, radiation, and the source of heat:
u t = α 2 u + q K u σ u 4 .
Note that all terms in Equation (2) are local, except the conduction term. In the case of Equation (1) in one space dimension, we apply to the α 2 u term the most common central difference Equation
2 x 2 u ( x i ) u ( x i + 1 ) u ( x i ) Δ x + u ( x i 1 ) u ( x i ) Δ x Δ x = u i 1 2 u i + u i + 1 Δ x 2
which is second order in Δ x , where i = 1 , , N and N is the overall number of nodes. By doing this, we are able to derive the spatially discretized version of the heat transfer Equation (2) in one space dimension as follows:
d u i d t = α u i 1 2 u i + u i + 1 Δ x 2 + q K u i σ u i 4 .
Let us now present the discretization of the heat transfer equation assuming that the quantities describing the properties of materials, namely α, k, c, and ρ, are functions of the space, rather than a fixed value. Now in one space dimension, instead of the α 2 u term, we have to deal with the term
1 c x ρ x x k x u x
In this case, the heat conduction equation can be discretized as follows
  c ( x i ) ρ ( x i ) u t x i = 1 Δ x   k x i + Δ x 2 u ( x i + Δ x ) u ( x i ) Δ x + k x i Δ x 2 u ( x i Δ x ) u ( x i ) Δ x .
Section 5 of the book [44] presents more details about this procedure for the case of underground reservoirs. The nodes are surrounded by cells, and k i , i + 1 is the heat conductivity between node i and its right neighbor. The discretized equation attains the following form
d u i d t = 1 c i ρ i Δ x k i , i + 1 u i + 1 u i Δ x + k i 1 , i u i 1 u i Δ x   +   q K u i σ u i 4 .
The dimensions of a cell, measured along its length and across its (typical) cross-section, are represented as Δ x and S . Since the volume of the cell can be given as V = S Δ x , the heat capacity is estimated as C i = c i ρ i V . The thermal resistance between the two neighboring nodes can be determined as R i , i + 1 Δ x / k i , i + 1 S . Based on these new quantities, one can obtain the following expression for the time derivative of each cell variable:
d u i d t = u i 1 u i R i 1 , i C i + u i + 1 u i R i + 1 , i C i + q K u i σ u i 4
It is not hard to generalize Equation (6) even more for the case of arbitrary number of neighbors to obtain the following spatially discretized version of Equation (2):
d u i d t = j i u j u i R i , j C i + q K u i σ u i 4 .
This ODE system is applicable to arbitrary (i.e., unstructured) grids with cells that may have different forms and characteristics. Naturally, uneven discretization might reduce spatial accuracy, but in this research, only cells with a rectangular form are employed.

2.2. The Leapfrog–Hopscotch Structure

It is necessary to discretize the space domain using a special, so-called bipartite mesh before applying the leapfrog–hopscotch or any other odd–even hopscotch algorithm.
Therefore, the mesh is split into two separate subgroups. The nodes or cells belonging to the first and second subgroups are designated as odd and even, respectively. The primary criterion is that, like on a checkerboard, all the near neighbors of the odd cells must be even and vice versa. We explain it using a 1D interval as an example x x 0 , x N , L = x N x 0 on which coordinates are used to create an equidistant grid x 0 , x 1 , , x N of nodes, so x j = x j 1 + Δ x , j = 1 , , N , Δ x = L / N . The time domain is t t 0 , t fin , and it is discretized as usual: t j = t 0 + j h , j = 1 , , T , h T = t fin t 0 , where h is the time step size.
The leapfrog–hopscotch (LH) space-time structure [33] is presented in Figure 1. The calculation consists of two half and several full-time steps. The calculation starts by taking a half-sized time step for the odd nodes based on the initial values. This “zeroth” stage is symbolized by a blue half-circle in the figure. After that, full-sized time steps (red circles in the figure) are taken strictly alternately for the even and odd nodes up until the end. At the end, the last time step (purple half-circle) should be cut in half for odd nodes to reach the same final time point as the even nodes. It is essential that when a stage-calculation is executed for node i, the most recently obtained values of the neighbors i − 1 and i + 1 are used in the equations.

2.3. The Mathematical Equations We Use

In this paper, we give only the formulas applied for Equation (7). Note that the discretization and the formulas for a 1D equidistant mesh can be found in the given references.
It is well known that the theta-method for the ODE y ' = f t , y has the Equation
y n + 1 = y n + 1 + Δ t θ f t n , y n + 1 θ f t n + 1 , y n + 1 ,
where θ 0 , 1 . If we adapt it for Equation (7), we obtain
u i n + 1 = u i n + θ A i n r i u i n Δ t K u i n Δ t σ u i n 4 + ( 1 θ ) A i n + 1 r i u i n + 1 Δ t K u i n + 1 Δ t σ u i n + 1 4 ,
where
A i n = Δ t j i u j n R i , j C i + Δ t q i ,   and   r i = Δ t j i 1 R i , j C i , i = 1 , , N , n = 0 , , T .
The first quantity mediates the accumulated effect of the temperatures of the cell i’s neighbors. The second quantity is the generalization of the often-used mesh ratio r = α Δ t Δ x 2 , more specifically for Equation (4), the r i = 2 r relation holds.
One should note that the sum in Equation (7) is separated into two terms: the first contains the actual cell variable, whereas the other, A i n , contains only its neighbors. The theta method given in Equation (9) is an implicit method if θ < 1 . It can be made explicit by the so-called pseudo-implicit trick: the neighbors in the second term at the r. h. s. of (9) must be taken at the n-th time level. Furthermore, 3 of the four powers of u i n + 1 in the nonlinear term can be replaced by u i n . With this, we obtain:
u i n + 1 = u i n + θ A i n r i u i n Δ t K u i n Δ t σ u i n 4 + 1 θ A i n r i u i n + 1 Δ t K u i n + 1 Δ t σ u i n + 1 u i n 4 ,
which can be rearranged as
u i n + 1 = u i n + A i n θ r i u i n + Δ t K u i n + Δ t σ u i n 4 1 + 1 θ r i + Δ t K + Δ t σ u i n 3 .
The value θ = 1 yields the explicit Euler method, which has a low CFL limit. If θ = 0 and q is non-negative, Equation (12) preserves the positivity of the temperature for arbitrary time step sizes. In this case, it can be considered as an adaptation of the so-called UPFD (unconditionally positive finite difference) scheme invented a decade ago [45]. Generally, smaller values of θ mean better stability, but it can imply worse accuracy as well. Now we fill the LH structure with this generalized theta-formula as follows.
In the pure conduction case (where K = 0 and σ = 0 ), the following theta values were obtained [39] during optimization: Stage 0: θ = 0 , all other stages: θ = 1 2 , which will be used everywhere in this work for the term A i n r i u i n . However, there is no reason to believe that for the convection and radiation terms, the optimal theta values are the same. In fact, in our last work [43], we examined the different treatments of these two terms, and according to our experience, only a few versions are competitive. For the convection term, θ = 1 2 is always the best choice, since our analytical calculations showed that it preserves second order convergence and unconditional stability at the same time. We currently do not have analytical proofs in the presence of the nonlinear term, but we found that the three theta values are worth examining, namely θ = 0 , 1 2 , 1 . We exemplify these by presenting the “zeroth” stage equations as follows.
1.
Pseudo-implicit treatment: θ = 0 for the radiation term, which yields:
u i 1 2 = u i 0 + A i 0 / 2 K Δ t u i 0 / 4 1 + r i + K Δ t / 4 + σ Δ t u i 0 3 / 2
2.
“Inside” treatment: θ = 1 for the radiation term, which means that it is taken into account explicitly, which yields
u i 1 2 = u i 0 + A i 0 / 2 K Δ t u i 0 / 4 σ Δ t u i 0 4 / 2 1 + r i + K Δ t / 4
3.
Mixed treatment with an equal share of the previous two treatments, that is, θ = 1 2 for the radiation term, which yields
u i 1 2 = u i 0 + A i 0 / 2 K Δ t u i 0 / 4 σ Δ t u i 0 4 / 4 1 + r i + K Δ t / 4 + σ Δ t u i 0 3 / 4

2.4. Other Explicit and Stable Methods

One classic example of a scheme that meets both of the criteria of explicitness and unconditional stability is the one of Dufort and Frankel (DF) [46] (p. 313). As this algorithm does not start itself, u i 1 has to be derived from u i 0 using some other approach. Here, we use the UPFD calculation method [43] for this purpose:
u i 1 = u i 0 + A i 0 1 + r i + Δ t K + Δ t σ u i 0 3
After this first step, the procedure is described by a simple formula for heat conduction, but the terms of convection and radiation can be handled in several ways. We present here only the most promising treatments.
4.
DF-D: the only place the Sigma and K terms appear is in the denominator:
u i n + 1 = 1 r i u i n 1 + 2 A i 1 + r i + 2 Δ t K + 2 Δ t σ u i n 3
5.
DF-M: the Sigma and K terms in a mixed way:
u i n + 1 = 1 r i u i n 1 + 2 A i Δ t K u i n Δ t σ u i n 4 1 + r i + Δ t K + Δ t σ u i n 3
6.
DF-KD: the K term appears only in the denominator, and the Sigma term in a mixed way:
u i n + 1 = 1 r i u i n 1 + 2 A i Δ t σ u i n 4 1 + r i + 2 Δ t K + Δ t σ u i n 3
7.
DF-SD: the Sigma term appears only in the denominator, and the K term in a mixed way:
u i n + 1 = 1 r i u i n 1 + 2 A i Δ t K u i n 1 + r i + Δ t K + 2 Δ t σ u i n 3
8.
Over half a century has passed since the discovery of the original odd–even hopscotch (OOEH) algorithm [47]. Its temporal and spatial organization has been described in [45]. After the first step by the FTCS formula (which is based on explicit Euler time discretization, so θ = 1 ) for the odd cells, the BTCS formula (which is based on implicit Euler time discretization) is used for the even cells. The labels odd and even are interchanged after each time step. We modify this method here to include the convection component, which is always considered at the new time level for enhanced stability. The radiation term is handled first explicitly and then implicitly [29]. These are the equations that are being used:
First   stage :                       u i n + 1 = 1 r i u i n + A i Δ t σ u i n 4 1 + Δ t K
Second   stage :                       u i n + 1 = u i n + A i new 1 + r i + Δ t K + Δ t σ + u i n 3 .
9.
NS-OOEH algorithm:
In addition, we try to alter the first (Explicit Euler) stage of the OOEH to make the terms for convection and radiation appear in the denominator to improve stability. The first stage, instead of (21), is as follows:
u i n + 1 = 1 r i u i n + A i 1 + Δ t K + Δ t σ ( u i n ) 3 .
The second stage is the same as in Equation (22).

2.5. Professional Solvers Used for Comparison Purposes

We compare our results with the direct, iterative, and programmable solvers available in student version of ANSYS 2023 R1 Academic software based on the finite element method (FEM).
(1)
The direct solver gives the exact solution to the system of equations defining the FE model. For the system represented by the equation [ K ] x = b , the exact solution is x = [ K ] 1 b , where [ K ] 1 is the inverse of the matrix [ K ] . The high computational cost of finding the inverse of a matrix means that direct solvers do not usually calculate the inverse, but use LU decomposition to solve the equation;
(2)
To provide an approximate solution within a certain convergence tolerance, iterative solvers assume an initial solution and iterate until they converge. Therefore, if the convergence tolerance is 0.01%, the solver will repeat until the difference between the current and past estimates of the solution is less than 0.01%.
The direct solutions tend to be more robust for complex systems and low-quality grids but require a relatively lot of memory. Iterative solvers are generally more efficient, because they use much less memory. Sometimes, however, they fail to converge even for a well-built model. Therefore, the requirements and capabilities should guide the choice of solver type. It is preferable to use an iterative solver to solve the model because that is ideal. The use of a direct solver is recommended if there is trouble with convergence. If the size of the model is quite large and the computer has limited RAM, it will need to use an iterative solution.
(3)
The programmable solver is controlled by a program and employes a mixture of direct and iterative solvers.
MATLAB solvers have been used for comparison purposes, namely ode15s, ode23t, ode23tb, ode23, ode45, and ode113. While implicit solutions are used for the other odes, it is known that odes 45, 23, and 113 employ explicit methods. Since the time step sizes cannot be calculated explicitly for the MATLAB solvers, we instead define the tolerances, beginning with a large value, such as Tol = 10 1 until an extremely small minimum value, which is Tol = 10 12 . We note that an additional solver, namely ode23s, is also available in MATLAB, but in our previous work [40,48], it is proved to be really slow for these problems; thus, we omitted it.
All the running times are measured on a desktop computer with an Intel Core i7-11700F (16 CPUs), and 64 GB RAM is used, whereas the program we used is MATLAB R2020b. The built-in tic-toc timer of that program is used to keep track of how long the algorithms have been running. We present a list of the names and a short description of the methods used in Table 1.

3. Numerical Simulation

3.1. Geometry and Material Properties

A piece of wall is considered with dimensions 1   m in the x and z directions and 0.02 m in the y direction, as can be seen in Figure 2. Two geometries are taken into consideration:
(a)
One-layer of brick is examined only for verification case;
(b)
Two-layers consisting of brick and rigid polyurethane foam insulation with a straight thermal bridge steel beam for running time measurements.
In the current work, the real material properties listed in Table 2 are taken into account. Note that these coefficients are constants inside a material, that is, they do not change with time, space, or temperature, but they have a discontinuity at the border of the materials.

3.2. Mesh Construction

We considered a dimension of wall 1 m × 1 m × 0.02 m in the calculations of Section 4. The y-direction is orthogonal to the surface of Figure 2A,B, and we use the approximation where no physical quantities change in that y-direction; as a result, this dimension from the mathematical point of view is irrelevant. In mathematical terms, this means that we only need to handle a two-dimensional problem (the cross-section) and can use Δ y i = 0.02 m. So, we have constructed three versions of meshes of size 1 m2, which means x , z 0 , 1 × 0 , 1 . The cells have a square or rectangular shape as shown in Figure 3.
We apply an equidistant grid to discretize the space variables in all cases. The computations were performed for three different uniform grids, namely, coarse = 40 × 40 = 1600, medium = 80 × 80 = 6400, and fine = 120 × 120 = 14,400. We have to note that the temperature is calculated at the nodes to be comparable to the ANSYS program, which is based on the finite element method, thus, we have a grid with these three grids.
The heat capacity of the cells around the nodes that are not at the boundary may be expressed in the case of a homogeneous material as C i = c i ρ i Δ x i Δ z i Δ y i . If the material properties are different, we can write it as follows C i = c i ρ i Δ x Δ z Δ y 2 + c i + 1 ρ i + 1 Δ x Δ z Δ y 2 . It is calculated at the upper, lower, left, and right borders by multiplying the heat capacity by a factor of 1/2 and at the corners by multiplying it by a factor of 1/4.
While the thermal resistance in the x-direction between two nodes has the approximate equation R x i Δ x k i S x i , where S x i is the surface element orthogonal to x, which can be given as S x = Δ y Δ z if the nodes are not at the boundary. The horizontal and vertical resistances can be calculated in the case of a homogeneous material as
R x Δ x k Δ z Δ y   and   R z Δ z k Δ x Δ y ,
but it is calculated at the upper, lower, left, and right borders as
R x Δ x k ( Δ z / 2 ) Δ y   and   R z Δ z k ( Δ x / 2 ) Δ y ,
respectively. If the material properties are different, we can write for the resistance between two nodes i and i + 1 that
R x i , i + 1 Δ x k i + k i + 1 2 Δ z Δ y ,
and since the node i + Nx is below the node i, we have
R z i , i + N x Δ z k i + k i + N x 2 Δ x Δ y ,
for the vertical resistance.
Also, it is calculated at the borders as
R x i , i + 1 Δ x k i + k i + 1 2 Δ z 2 Δ y   and   R z i , i + N x Δ z k i + k i + N x 2 Δ x 2 Δ y

3.3. The Initial and the Boundary Conditions

For running time measurements, zero Neumann boundary conditions have been applied to all borders, preventing the passage of conductive heat at the boundaries. To do this, we set zero for the matrix components describing heat conduction across the border by setting the relevant resistances to infinity.
On the left and right sides of the wall, convective and radiative heat transfer have been assumed. The interior components cannot lose or absorb heat via radiation or convection, and there is no heat source in those elements. Table 3 shows that elements on the left and right sides may transfer heat through radiation and convection in the x direction. The final time t fin = 20 , 000 s indicates the end of the time interval that is examined; additionally, the time step size is specified in seconds.
As can be seen in Table 3, we have used numbers from publication [42] for the convection heat transfer coefficient hc. The Stefan–Boltzmann constant is a universal constant for radiation: 5.67 10 8 W m 2 K 4 . Since the surface is not a perfect black body, we have to multiply the Stefan-Boltzmann constant by an emissivity constant in order to get results that are more in line with realistic values for σ . The value of q is estimated for the definition of “heat source”, which includes solar radiation in the table below. It is assumed that the ambient air temperature on the left side is 17     ° C 290     K .
Because of the nonzero temperature of the air u a , the term q also includes heat gain from convection; this allows us to calculate the value of q as follows. The values of the coefficients in our Equations (2) and (3) are calculated using the values we obtain:
K = h c c ρ Δ x 2 , σ = σ c ρ Δ x 2 , q = q c ρ Δ x 2 + h c c ρ Δ x 2 u a
We assumed that the left and right elements have the following heat sources as follows:
In   terms   of   the   left - hand   side :   q = 1 c ρ Δ x 2 × 360.95 W m 2 + h c c ρ Δ x 2 × 290     K
In   terms   of   the   left - hand   side :   q = 1 c ρ Δ x 2 × 435.39 W m 2 + h c c ρ Δ x 2 × 313     K 2 u v .
The constant initial temperature is prescribed: u x , z , t = 0 = 290     K .

3.4. Verification Using Analytical Solutions

For verification, we used three different cases. In the first one, we considered a one-layer wall made of homogeneous material (see Figure 2A) and applied a sinusoidal initial condition with a zero Dirichlet boundary condition to examine the conduction term only. Then in the second and third cases, we considered spatially homogeneous initial temperature to exclude conduction, which means we have only one temperature in the system, that is, the behavior can be described by an ODE. In the second case, we took into account only convection, which is perpendicular to the surface and happens in the y direction, and in the third case, we considered only radiation perpendicular to that surface. After that, we compared our results in MATLAB and ANSYS with the analytical solution.
  • First verification (for conduction):
1.
Two sine functions are multiplied to provide the initial condition:
u ( x , z , t = 0 ) = sin ( k x π x ) sin ( k z π z )
2.
Zero Dirichlet boundary conditions are used:
u ( x = 0 , z , t ) = u ( x = 1 , z , t ) = u ( x , z = 0 , t ) = u ( x , z = 1 , t ) = 0
3.
One can quickly verify that the problem’s analytical solution is
u ( x , z , t ) = sin ( k x π x ) sin ( k z π y ) e α ( k x 2 + k z 2 ) π 2 t
where the wave numbers are fixed to k x = 1 ,   k z = 1 , and substituting the physical properties of the brick, we obtained α. We used the analytical solution Equation (23) at t fin = 2000 s as the reference solution.
  • Second verification (for convection and heat generation)
Since the temperature is spatially homogeneous, we have the simple ODE, valid for each node:
d u d t = K u + Q
Its analytical solution is:
u ( x , z , t ) = Q K + u 0 Q K e K t
where the initial condition here is a constant temperature that equals to 290 K.
  • Third verification (for radiation)
We have the ODE:
d u d t = σ u 4
with the analytical solution:
u ( x , z , t ) = u 0 3 + 3 σ T 1 3
The error is defined as the largest difference in absolute terms between the reference temperature u i ref , which is the analytical solution and the temperature u i num obtained by the studied numerical method at t fin = 2000   s :
M a x E r r o r = max 1 i N u i ref ( t fin ) u i num ( t fin ) .
For all systems, the obtained results in the case of MATLAB and ANSYS are very similar to the analytical solution, and the error is below 10−5. This means that the MATLAB code using the ode15s solver, as well as the ANSYS solver, has been verified. The temperature contour is presented in the Supplementary Materials (see Figure S1) as a function of the x- and z-coordinates in the case of ANSYS and the analytical solution. The temperature as a function of the x-coordinate for z = 0.5 m was plotted in the case of ode15s, ANSYS, and the analytical solution for the coarse grid. These lines are so close to one another that they are indistinguishable by the naked eye, (see Figure S2 in the Supplementary Materials). For other system sizes, the same behaviour is experienced. In the case of the second and third verification, the maximum error is around 10−9 between ode15s and the analytical solution, but it is slightly larger than 10−5 between ANSYS and the analytical solution, as shown in Table 4. Since the MATLAB ode15s solver is proven to be more accurate in all cases, we choose this as a reference solution when we measure and compare the running times in the next section.

4. Results: Comparison of Performances by Measuring the Running Times

4.1. Comparison with MATLAB Methods and ANSYS Solvers for the Coarse Mesh System

The two-layer wall (brick + insulation) with the straight thermal bridge (see Figure 2B) is simulated using the coarse mesh and the conditions in Section 3.3, and values from Table 2. The maximum errors are shown in Figure 4 as a function of the time step size. The comparison of the running times of the tested methods and three solvers of ANSYS are shown in Figure 5, and the temperature distributions are shown in Figure 6. We can see that for the coarse mesh, the LH pseudo-implicit treatment of radiation was the most accurate, but two of the DF versions and the three solvers of ANSYS are not really accurate. Regarding speed, we can say that the explicit and stable schemes are clearly faster than ANSYS, but they are more efficient than the MATLAB routines only if small or moderate accuracy is required.

4.2. Comparison with MATLAB Methods and ANSYS Solvers for the Moderate Mesh System

Now a moderately fine resolution mesh is applied, and the same conditions and values are used as for the previous coarse grid. Most of the proposed methods, especially the LH-PI is again better than ANSYS solvers and MATLAB routines for both accuracy and speed. Figure 7 compares the running times of the three ANSYS solvers with the tested methods in MATLAB. (We note that the errors as a function of the time step size can be seen in Figure S3 in the Supplementary Materials).

4.3. Comparison with MATLAB Methods and ANSYS Solvers for the Fine Mesh System

For the fine mesh, the errors are plotted in Figure 8 as a function of the running time and in Figure S4 (Supplementary Materials) as a function of the time step size. We notice that the accuracy of the ANSYS solvers is positively affected by the smoothing of the mesh, the errors are decreasing with the time step size. In spite of this, some of the proposed explicit methods coded in MATLAB remain the best in both accuracy and speed. The solvers of ANSYS are very slow if we have a big system, but the built-in routines of MATLAB are getting slower at a larger rate with increasing system size, so we think that for even finer mesh, they would be the slowest. The leapfrog–hopscotch with the pseudo-implicit treatment of the radiation term (LH-PseudoImp) and the Dufort–Frankel schemes with the pseudo-implicit treatment of both the convection and the radiation term (DF-D) are the two best methods. Their advantages are in fact, increasing with increasing system size. That is why, in the next section, we choose these two methods to make a long-term simulation.
For mesh-independence, a horizontal line is considered in the middle of the thermal bridge (z = 0.225 m), and along this line, the temperature of the LH pseudo-Imp method is plotted for all three meshes. The result is presented in the Supplementary Materials (see Figure S5). As one can see, the values of the medium mesh are close to those of the fine mesh, where the maximum difference was quite small 0.097 in Kelvin units, so to reduce the computational cost, we can adapt the medium mesh in long-term calculations.
We also compared three of the solvers, the iterative ANSYS, the LH-PseudoImp, and ode15s using the medium mesh. The maximum difference between LH-PseudoImp and ode15s is equal to 7.3 × 10−7 in Kelvin units, but ANSYS has a little noticeable difference from them, which is equal to 0.2 K. The figure is presented in the Supplementary Materials (see Figure S6).
Finally, we examined how the running times depend on the system size. We obtained that the running time of our methods increases linearly with the number of nodes, but it increases faster than linearly for other solvers. Figure S7 in the Supplementary Materials shows the concrete data for one of our best methods with one of the MATLAB routines ode23t, and ANSYS solvers.
Nevertheless, we absolutely do not intend to suggest that ANSYS and similar solvers are redundant since there are faster numerical algorithms than those they employ. This simulation software automatically generates the mesh and presents graphically the results, whereas if someone chooses to write their own code, it is tedious work, especially in the case of complicated geometries. Moreover, the explicit and stable methods have not been applied to cases involving the simulation of the motion of the air, for example, for the problem of air leakage in building components, whereas the simulation software routinely handles these problems. However, we think that in the already well-established cases, such as those types studied in this paper, the simulation software could incorporate the proposed algorithms as an option can be chosen by the users or by an artificial intelligence.

5. Long-Term Simulations

5.1. Geometry and Mesh

The simulations were performed for a standard residential wall in Miskolc, Hungary, and the outdoor temperature, convection coefficient, and solar radiation values for Miskolc were used in the simulations where the data was taken from the website every 3 h. Then we used linear interpolation to calculate the values for every 100 s.
The month of January is the coldest. The highest temperature measured in January 2022 was 11 °C, whereas the lowest was −6 °C. Depending on the topography, the predominant wind might blow in various different directions. The maximum measured wind speed was 10.8 m/s [49].
We applied the best methods, namely, the LH-PI and DF-D, to calculate temperatures and heat losses across the wall, with a fixed time step size of Δ t = 100 s, where the total time was T = 31 24 3600 = 2 , 678 , 400 s. We chose this time step size because from Figure 4, it can be seen that the error of the best methods is around 0.01 °C, so sufficient accuracy is reached. We stress again that the CFL limit for explicit Runge-Kutta methods is orders of magnitude smaller than this time step size. The difference between the DF-D and LH-pseudoImp was continuously monitored during the simulation and found to be very small, so the presented results are for the LH-pseudoImp method. The running times for long-time simulation was 145.5 s for the case of two-layers with the straight thermal bridge case on the same computer which was used in the previous section.
Figure 9 illustrates the different models of the wall. First, one layer has been used, which is only brick; the dimension of the wall height, width, and thickness is (0.45, 1, 1 m) as shown in Figure 9A. In the second model, there are two layers: the brick wall has the same dimension as in the first model, and in addition to it, there is an insulation layer with a thickness of 0.15 m, as shown in Figure 9B. The third model is the same as the second one with a straight thermal bridge. The width of the thermal bridge is the same as that of the insulator, as shown in Figure 9C. The horizontal position of the thermal bridge is between x = 0.45 m and x = 0.6 m. The vertical position of the top of the thermal bridge is z = 0.75 m (25 cm from the bottom of the wall), and the thickness of the bridge is 5 cm in the z direction.
The fourth model contains the same two layers as the second one, but the bar thermal bridge has a curved shape. The thermal bridge is made up of three straight bars, two horizontal and one vertical. The first horizontal bar is connected from the outside at the same vertical position as the straight bridge above. It goes horizontally with a length 0.05 m. It is connected to the vertical bar, which has a length 0.3375 m in the z-direction. It then connects to the other horizontal bar, which has the same dimensions as the first horizontal bar, whereas the vertical position of the bottom of this piece is 0.5125 m as shown in Figure 9D. The medium mesh is applied for all cases, but in the first case, it means Δ x = 0.0057 , whereas in the other cases, the widths of the cells are larger, Δ x = 0.0076 , since the wall is thicker.
We are going to examine which thermal bridge yields a higher rate of heat transfer and more extra heating cost. It is a nontrivial question since, in the case of the straight bridge, the way of heat is shorter, but in the case of the bent bridge, the average conductivity of the wall is higher.

5.2. Initial and Boundary Conditions

Zero Neumann boundary conditions have been applied to all borders and wall cases, and convective as well as radiative heat transfer are assumed on the right and left sides. As in Section 4, the interior parts cannot lose or absorb heat via convection, radiation, and heat source. As it is indicated in Table 5, elements on the right and left sides can transfer heat through convection and radiation in the x direction. The left elements have constant values for the convection, radiation, and heat source parameters, whereas the right elements have changing values depending on the external conditions. In this part, the final time is tfin = 2,678,400 s, which indicates the simulation of a full winter month.
We obtain the values of the coefficients in our equations as follows: [42],
K = h c c ρ Δ x , σ = σ c ρ Δ x , q = q c ρ Δ x + h c c ρ Δ x u a .
We also supposed that the left and right elements have the heat sources as follows:
For   the   interior   side   elements : q l = 1 c ρ × q l + h c l c ρ Δ x × 295     K
For   the   external   side   elements : q r ( t ) = 1 c ρ × q r t + h c r t c ρ Δ x × u r t
And   K t = h c r t c ρ Δ x , σ = σ r c ρ Δ x , q t = q r t c ρ Δ x + h c r t c ρ Δ x u r t
where q l = ε l σ l ( 295 ) 4 and q r t = α sun G c r t + α Low ε r σ r u r t 4 [50].
The convection heat transfer coefficient for outside elements as a function of air velocity is estimated as follows [51]:
h c r t = 0.6 + 6.64 v t
v t : The air velocity is taken for each 100 s in January [m/s].
u r t : The outside air temperature for each 100 s in January month [°C].
u l : The inside air temperature [°C] on the left side.
G c r t : The solar radiation was taken for each 100 s in January [W/m2].
α sun : The absorptivity of brick surface to solar radiation.
α Low : The absorptivity of brick surface to low-temperature thermal radiation.
The environment air temperature is taken to be 22     ° C 295     K inside, and changing depending on weather conditions outside.
We calculated the initial temperature inside the wall using the assumption that before the simulation time, a stationary heat flow with constant flux evolved between the given boundary values of the internal and external air temperatures. In the case of one-layer, it yields a linear function of the x variable for the initial condition:
u x , z , t = 0 = ( u r , initial u l ) x / L x + u l
where u r , initial = 278     K .
For the remaining three cases, the assumption of stationary heat conduction with the initial values at the boundaries implies that we have to use two linear functions of the x variable for the initial condition:
For the brick part: u x , z , t = 0 = ( u mid u l ) x / L b + u l
For the insulator part:
u x , z , t = 0 = ( ( u r , initial u mid ) x / L ins ) ( ( u r , initial u mid ) L b / L i n s ) + u mid where   u mid = u l ( q flux L b / k b ) and   q flux = ( u l u r , initial ) / ( ( L b / k b ) + ( L ins / k ins ) )

5.3. Result for the Simulation of the Wall

The results will be displayed for specific points depicted in Figure 10, and the subfigures contain the following information.
(A)
There are 3 points denoted for the one-layer wall: on the interior surface, at the middle of the wall, and on the exterior surface;
(B)
Two-layer wall (brick + insulator) with five points on the interior surface, the middle of the brick, the border of the two layers, the middle of the insulation, and the exterior surface;
(C)
Two-layers with a straight thermal bridge with 6 points: on the interior surface, the middle of the brick, the border of the two layers, the exterior surface of the straight thermal bridge, the middle of the insulation, and the exterior surface of the insulation;
(D)
Two-layers with the bent thermal bridge with six points: on the interior surface, middle of the brick, the border of the two layers with bent thermal bridge, the exterior surface of the bent thermal bridge, the middle of the insulation, and the exterior surface of the insulation.

5.3.1. One Layer (Brick)

The temperatures as a function of the time are presented in Figure 11. Note that without the presence of the insulator, the temperature of the inner surface of the wall is initially equal to the temperature of the internal air (22 °C), then decreases because of the cooling effect of the wall and the outside cold weather. As for the temperature of the middle of the wall, it is low at first, and then it rises due to heat transfer by conduction from the inside. The external surface temperature is greater than the external air temperature due to the heat flowing to the surface through conduction.

5.3.2. Two Layers (Brick and Insulation) without Thermal Bridge

The temperature-time functions are presented in Figure 12. The temperature in the brick part is changing very slowly, but in the insulator part it is changing very sharply. We note that in the presence of an insulator, the temperature of any point of the brick follows the temperature of external air only very slightly. As for the temperature of the middle of the insulator, it rises and falls due to the effects of external conditions, and the minimum values of the outer surface temperatures of the wall are slightly less than the inside air temperature because the insulator limits the flow of heat from inside to outside; for higher values, it is greater owing to the effects of solar radiation. The comparison between the one-layer and two-layer wall cases in terms of the final temperature is presented in the Supplementary Materials (see Figure S8).

5.3.3. Two Layers with Straight Thermal Bridging

Figure 13 clearly shows that the temperature of the middle of the wall and the temperature of the interior surface of the insulator decreased slightly compared to the previous case due to the heat loss caused by the straight thermal bridge. Moreover, the external temperature of the thermal bridge is higher than the temperature of the exterior surface of the insulator due to the flow of heat from inside to outside as well as the physical properties of the material.

5.3.4. Two Layers with Bent Thermal Bridging

In this case, which includes the thermal bridge created by the bent thermal bridge, the mid-insulator temperature rises (see Figure S9 in the Supplementary Materials). Figure 14 shows the contours of both types of thermal bridge cases, where we note the way of the heat in the bent bridge is longer than in the straight bridge.
In Figure 15, we compared the temperatures of the external points of the two thermal bridges with the external points of the one- and two-layer cases. We concluded that the straight thermal bridge allows the heat to flow faster than the bent thermal bridge because the passage of heat in the bent bridge is longer than in the straight one.
Figure 16 shows a comparison between the losses in energy. One can see the largest thermal loss in the case of a one-layer wall, and these losses are less with the presence of the insulator. When a thermal bridge is present, these losses are larger than without a thermal bridge, and the losses are slightly larger in the case of a straight thermal bridge than for the bent bar.

5.3.5. Calculation of Heat Loss

The total energy requirement Q in kWh units was calculated from Equation (25)
Q = t = 1 T i m e Δ t i = 1 N Q cell , i ( t ) ,
where Q cell = q * cell S and q * is the current density of the heat loss, which can be calculated as q * cell = ( T i + 1 T i ) K cond Δ x .
Δ t is the time step size, and S is the area of the wall in the direction of heat loss.
Now we can calculate the cost of the energy used to make up the energy loss by the multiplication of the current price of electricity (per kWh) used by the amount of heat losses. Table 6 below shows the heat loss and costs of energy consumption.

5.4. Simulation of the Coldest Day of the Month

For the design of the heating load, we adopt the hardest weather conditions, and the design is based on them. Therefore, we will review the coldest day of the month, which is 25 January, where the lowest recorded temperature is −6 °C. We also applied the previous best methods to calculate temperatures and heat losses across the wall, but with a new time step size of Δ t = 60 s in this simulation, to track and record the temperature change every minute because the total simulation time here is short, only T = 24 × 3600 = 86,400 s. The temperature as a function of the hours of the day for the exterior points in all cases of the wall, and a comparison between the losses in energy on the coldest day of the month are presented in the Supplementary Materials (see Figures S10 and S11). One can see the largest losses in the case of one layer, and Table 7 contains the cost of energy consumption on this day.

6. Conclusions

The purpose of this article was to develop and test a simulation methodology based on fundamental physical principles and laws to study transient heat transfer in the wall. First, we numerically investigated transient heat transfer in a two-dimensional wall without an insulator, with an insulator, and two types of thermal bridges. For our purposes, we applied and tested nine explicit algorithms coded by us in MATLAB and three solvers included in the commercial software called ANSYS. Then the running time was measured for our methods and compared with the ANSYS solvers and six built-in MATLAB routines. We used three mesh sizes: 40 × 40, 80 × 80, and 120 × 120, which we applied to walls with 1600, 6400, and 14,400 cells, respectively. Three simple analytical solutions of the heat equation were used with an equidistant mesh for verification in the case of homogeneous material properties (one brick layer). All the methods used and the ANSYS solvers are confirmed to be convergent. These experiments suggest that our methods are better than all ANSYS solvers and MATLAB routines, whereas ANSYS was less accurate and slower, and it was observed that the best performance was achieved by the leapfrog–hopscotch and the Dufort–Frankel algorithms with the pseudo-implicit treatment of the nonlinear radiation term. Therefore, these two methods were applied to real problems, and a long-term simulation of four cases was performed. The temperature distribution and total heat losses of all cases were calculated. We found the straight thermal bridge to be energetically worse than others, and the total heat loss during the month (one-layer, two-layer, two-layer with a straight thermal bridge, and two-layer with a bent thermal bridge) was, respectively, 19.14, 1.99, 5.29, and 5.01 kWh for a 1 m2 wall surface.
We can conclude that the numerical simulation methodology is established in this paper. In our next works, we will use these solvers and approaches to study real problems like real thermal bridges in blocks of flats, and we will use physical experiments as well as standards such as ISO 10211:2017 and ISO 14683:2017 to validate them.

Supplementary Materials

The following supporting information can be downloaded at: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/en16124604/s1, Figure S1: The temperature distribution contour for the sinusoidal initial temperature (left) for ANSYS and analytical solution on (right) in Kelvin units; Figure S2: The temperature comparison for the sinusoidal initial temperature in the case of the coarse mesh; Figure S3: The maximum errors as a function of the time step size for the tested methods in the case of the medium mesh; Figure S4: The maximum errors as a function of the time step size for the tested methods in the case of the fine mesh; Figure S5: Comparison between three types of mesh for the LH-PseudoImp method temperature; Figure S6: Comparison between three solvers for the medium mesh; Figure S7: Running time as a function of the total number of cells for the examined method and other solvers, where the left axis refers to LH-PseudoImp, and the right one refers to ode23t solver of MATLAB and ANSYS; Figure S8: The temperature in Celsius units for the long-term simulation, for one layer and two layers at the end of the last day; Figure S9: The temperature in Celsius units as a function of time in days for the long-term simulation of the wall with the bent thermal bridging; Figure S10: The temperature in Celsius units as a function of the hours of the day for the one day simulation of the wall on the exterior point of the thermal bridge (straight and bent) compared with the one-layer and two-layer cases; Figure S11: The rate of the total heat loss in W units with hours of the day for all cases of the one-day wall simulation.

Author Contributions

Conceptualization and resources: E.K.; methodology and supervision: E.K. and B.B.; software and visualization: I.O. and A.H.A.; investigation: I.O.; writing—original draft preparation, I.O.; writing—review and editing: I.O., E.K. and B.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data available on request from the first author.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

SymbolsGreek Symbols
cSpecific heat [J/(kgK)]αThermal diffusivity [m2/s]
Δ t Time step size [s]ΔDifference
hcHeat transfer coefficient [W/(m2K)]ΡMass density [kg/m3]
KConvection coefficient [1/s]ΣCoefficient of the radiation term [s−1K−3]
kThermal conductivity [W/(m·K)] σ realistic values of the non-black body [W/(m2·K4)]
QHeat transfer rate [W]Subscripts
q heat generation [W/m2]AAmbient air
qHeat source rate [K/s]LLeft side
ttime [s]RRight side
uTemperature [K]b, insBrick and Insulation
LThickness [m]cconvection

References

  1. Berggren, B.; Wall, M. State of Knowledge of Thermal Bridges—A Follow up in Sweden and a Review of Recent Research. Buildings 2018, 8, 154. [Google Scholar] [CrossRef] [Green Version]
  2. Liu, C.; Sharples, S. A Review of Building Energy Retrofit Measures, Passive Design Strategies and Building Regulation for the Low Carbon Development of Existing Dwellings in the Hot Summer—Cold Winter Region of China. Energies 2023, 16, 4115. [Google Scholar] [CrossRef]
  3. Kusuda, T. Fundamentals of Building Heat Transfer. J. Res. Natl. Bur. Stand. 1977, 82, 97. [Google Scholar] [CrossRef] [PubMed]
  4. Dai, B.; Qi, H.; Liu, S.; Zhong, Z.; Li, H.; Song, M.; Ma, M.; Sun, Z. Environmental and Economical Analyses of Transcritical CO2 Heat Pump Combined with Direct Dedicated Mechanical Subcooling (DMS) for Space Heating in China. Energy Convers. Manag. 2019, 198, 111317. [Google Scholar] [CrossRef]
  5. Cao, X.; Dai, X.; Liu, J. Building Energy-Consumption Status Worldwide and the State-of-the-Art Technologies for Zero-Energy Buildings during the Past Decade. Energy Build. 2016, 128, 198–213. [Google Scholar] [CrossRef]
  6. Li, T.; Xia, J.; Chin, C.S.; Song, P. Investigation of the Thermal Performance of Lightweight Assembled Exterior Wall Panel (LAEWP) with Stud Connections. Buildings 2022, 12, 473. [Google Scholar] [CrossRef]
  7. Lima, S.A.; Kamrujjaman, M.; Islam, M.S. Numerical Solution of Convection--Diffusion--Reaction Equations by a Finite Element Method with Error Correlation. AIP Adv. 2021, 11, 085225. [Google Scholar] [CrossRef]
  8. EN ISO 13789:2017; Swedish Standards Institute Thermal Performance of Buildings—Transmission and Ventilation Heat Transfer Coefficients—Calculation Method. International Organization for Standardization: Geneva, Switzerland, 2017. Available online: https://www.iso.org/standard/65713.html (accessed on 20 April 2023).
  9. EN ISO 10211:2017; Swedish Standards Institute Thermal Bridges in Building Construction—Heat Flows and Surface Temperatures—Detailed Calculations. International Organization for Standardization: Geneva, Switzerland, 2017. Available online: https://www.iso.org/%0Astandard/65710.html (accessed on 20 April 2023).
  10. EN ISO 14683:2017; Swedish Standards Institute Thermal Bridges in Building Construction—Linear Thermal Transmittance—Simplified Methods and Default Values. International Organization for Standardization: Geneva, Switzerland, 2017. Available online: https://www.iso.org/standard/65706.html (accessed on 20 April 2023).
  11. Alhawari, A.; Mukhopadhyaya, P. Thermal Bridges in Building Envelopes—An Overview of Impacts and Solutions. Int. Rev. Appl. Sci. Eng. 2018, 9, 31–40. [Google Scholar] [CrossRef]
  12. Borelli, D.; Cavalletti, P.; Marchitto, A.; Schenone, C. A Comprehensive Study Devoted to Determine Linear Thermal Bridges Transmittance in Existing Buildings. Energy Build. 2020, 224, 110136. [Google Scholar] [CrossRef]
  13. Karabulut, K.; Buyruk, E.; Fertelli, A. Numerical Investigation of the Effect of Insulation on Heat Transfer of Thermal Bridges with Different Types. Therm. Sci. 2016, 20, 185–195. [Google Scholar] [CrossRef]
  14. Mátyás, L.; Barna, I.F. General Self-Similar Solutions of Diffusion Equation and Related Constructions. arXiv 2021, arXiv:2104.09128v1. [Google Scholar]
  15. Amiri, S.; Mazaheri, M.; Bavandpouri Gilan, N. Introducing a New Method for Calculating the Spatial and Temporal Distribution of Pollutants in Rivers. Int. J. Environ. Sci. Technol. 2021, 18, 3777–3794. [Google Scholar] [CrossRef]
  16. Jejeniwa, O.A.; Gidey, H.H.; Appadu, A.R. Numerical Modeling of Pollutant Transport: Results and Optimal Parameters. Symmetry 2022, 14, 2616. [Google Scholar] [CrossRef]
  17. Yadav, V.S.; Singh, A.; Maurya, V.; Rajpoot, M.K. New RK Type Time-Integration Methods for Stiff Convection–Diffusion–Reaction Systems. Comput. Fluids 2023, 257, 105865. [Google Scholar] [CrossRef]
  18. Kumar, V.; Chandan, K.; Nagaraja, K.V.; Reddy, M. V Heat Conduction with Krylov Subspace Method Using FEniCSx. Energies 2022, 15, 8077. [Google Scholar] [CrossRef]
  19. Mbroh, N.A.; Munyakazi, J.B. A Robust Numerical Scheme for Singularly Perturbed Parabolic Reaction-Diffusion Problems via the Method of Lines. Int. J. Comput. Math. 2022, 99, 1139–1158. [Google Scholar] [CrossRef]
  20. Essongue, S.; Ledoux, Y.; Ballu, A. Speeding up Mesoscale Thermal Simulations of Powder Bed Additive Manufacturing Thanks to the Forward Euler Time-Integration Scheme: A Critical Assessment. Finite Elem. Anal. Des. 2022, 211, 103825. [Google Scholar] [CrossRef]
  21. Beuken, L.; Cheffert, O.; Tutueva, A.; Butusov, D.; Legat, V. Numerical Stability and Performance of Semi-Explicit and Semi-Implicit Predictor--Corrector Methods. Mathematics 2022, 10, 2015. [Google Scholar] [CrossRef]
  22. Fedoseev, P.; Pesterev, D.; Karimov, A.; Butusov, D. New Step Size Control Algorithm for Semi-Implicit Composition ODE Solvers. Algorithms 2022, 15, 275. [Google Scholar] [CrossRef]
  23. Ndou, N.; Dlamini, P.; Jacobs, B.A. Enhanced Unconditionally Positive Finite Difference Method for Advection--Diffusion--Reaction Equations. Mathematics 2022, 10, 2639. [Google Scholar] [CrossRef]
  24. Ji, Y.; Xing, Y. Highly Accurate and Efficient Time Integration Methods with Unconditional Stability and Flexible Numerical Dissipation. Mathematics 2023, 11, 593. [Google Scholar] [CrossRef]
  25. Appadu, A.R. Performance of UPFD Scheme under Some Different Regimes of Advection, Diffusion and Reaction. Int. J. Numer. Methods Heat Fluid Flow 2017, 27, 1412–1429. [Google Scholar] [CrossRef] [Green Version]
  26. Karahan, H. Unconditional Stable Explicit Finite Difference Technique for the Advection—Diffusion Equation Using Spreadsheets. Adv. Eng. Softw. 2007, 38, 80–86. [Google Scholar] [CrossRef]
  27. Sanjaya, F.; Mungkasi, S. A Simple but Accurate Explicit Finite Difference Method for the Advection-Diffusion Equation. Proc. J. Phys. Conf. Ser. 2017, 909, 12038. [Google Scholar] [CrossRef]
  28. Pourghanbar, S.; Manafian, J.; Ranjbar, M.; Aliyeva, A.; Gasimov, Y.S. An Efficient Alternating Direction Explicit Method for Solving a Nonlinear Partial Differential Equation. Math. Probl. Eng. 2020, 2020, 9647416. [Google Scholar] [CrossRef]
  29. Harley, C. Hopscotch Method: The Numerical Solution of the Frank-Kamenetskii Partial Differential Equation. Appl. Math. Comput. 2010, 217, 4065–4075. [Google Scholar] [CrossRef]
  30. Al-Bayati, A.Y.; Manaa, S.A.; Al-Rozbayani, A.M. Comparison of Finite Difference Solution Methods for Reaction Diffusion System in Two Dimensions. AL-Rafidain J. Comput. Sci. Math. 2011, 8, 21–36. [Google Scholar] [CrossRef] [Green Version]
  31. Nwaigwe, C.; An Unconditionally Stable Scheme for Two-Dimensional Convection-Diffusion-Reaction Equations. No. January 2022. Available online: https://www.researchgate.net/publication/357606287_An_Unconditionally_Stable_Scheme_for_Two-Dimensional_Convection-Diffusion-Reaction_Equations (accessed on 5 June 2023).
  32. Savović, S.; Drljača, B.; Djordjevich, A. A Comparative Study of Two Different Finite Difference Methods for Solving Advection--Diffusion Reaction Equation for Modeling Exponential Traveling Wave in Heat and Mass Transfer Processes. Ric. Mat. 2021, 71, 245–252. [Google Scholar] [CrossRef]
  33. Nagy, Á.; Omle, I.; Kareem, H.; Kovács, E.; Barna, I.F.; Bognar, G. Stable, Explicit, Leapfrog-Hopscotch Algorithms for the Diffusion Equation. Computation 2021, 9, 92. [Google Scholar] [CrossRef]
  34. Kovács, E.; Nagy, Á.; Saleh, M. A Set of New Stable, Explicit, Second Order Schemes for the Non-Stationary Heat Conduction Equation. Mathematics 2021, 9, 2284. [Google Scholar] [CrossRef]
  35. Jalghaf, H.K.; Kovács, E.; Majár, J.; Nagy, Á.; Askar, A.H. Explicit Stable Finite Difference Methods for Diffusion-Reaction Type Equations. Mathematics 2021, 9, 3308. [Google Scholar] [CrossRef]
  36. Kareem Jalghaf, H.; Omle, I.; Kovács, E. A Comparative Study of Explicit and Stable Time Integration Schemes for Heat Conduction in an Insulated Wall. Buildings 2022, 12, 824. [Google Scholar] [CrossRef]
  37. Pirasaci, T. Investigation of Phase State and Heat Storage Form of the Phase Change Material (PCM) Layer Integrated into the Exterior Walls of the Residential-Apartment during Heating Season. Energy 2020, 207, 118176. [Google Scholar] [CrossRef]
  38. Nagy, Á.; Saleh, M.; Omle, I.; Kareem, H.; Kovács, E. New Stable, Explicit, Shifted-Hopscotch Algorithms for the Heat Equation. Math. Comput. Appl. 2021, 26, 61. [Google Scholar] [CrossRef]
  39. Saleh, M.; Nagy, Á.; Kovács, E. Part 1: Construction and Investigation of New Numerical Algorithms for the Heat Equation. Multidiszcip. Tud. 2020, 10, 323–338. [Google Scholar] [CrossRef]
  40. Mahmoud, S.; Ádám, N.; Endre, K. Construction and Investigation of New Numerical Algorithms for the Heat Equation: Part III. Multidiszcip. Tud. 2020, 10, 349–360. [Google Scholar]
  41. Alawadhi, E.M. Finite Element Simulations Using ANSYS; CRC Press: Boca Raton, FL, USA, 2009. [Google Scholar]
  42. Holman, J. Heat Transfer, 10th ed.; McGraw-Hill Education: New York, NY, USA, 2010. [Google Scholar]
  43. Askar, A.H.; Omle, I.; Kovács, E.; Majár, J. Testing Some Different Implementations of Heat Convection and Radiation in the Leapfrog-Hopscotch Algorithm. Algorithms 2022, 15, 400. [Google Scholar] [CrossRef]
  44. Munka, M.; Pápay, J. 4D Numerical Modeling of Petroleum Reservoir Recovery; Akadémiai Kiadó: Budapest, Hungary, 2001; ISBN 9630578433. [Google Scholar]
  45. Chen-Charpentier, B.M.; Kojouharov, H.V. An Unconditionally Positivity Preserving Scheme for Advection—Diffusion Reaction Equations. Math. Comput. Model. 2013, 57, 2177–2185. [Google Scholar] [CrossRef]
  46. Amando, R.S. Numerical Computation of Internal and External Flows, Volume 1: Fundamentals of Numerical Discretization; Hirsch, C., Ed.; John Wiley & Sons Ltd.: Hoboken, NJ, USA, 1989. [Google Scholar]
  47. Gourlay, A.R.; McGuire, G.R. General Hopscotch Algorithm for the Numerical Solution of Partial Differential Equations. IMA J. Appl. Math. 1971, 7, 216–227. [Google Scholar] [CrossRef]
  48. Kovács, E. New Stable, Explicit, First Order Method to Solve the Heat Conduction Equation. arXiv 2019, arXiv:1908.09500. [Google Scholar] [CrossRef]
  49. Weather-Online Website. Available online: https://www.worldweatheronline.com/miskolc-weatherhistory/miskolc/hu.aspx (accessed on 20 November 2022).
  50. Askar, A.H.; Kovács, E.; Bolló, B. ANN modeling for thermal load estimation in a cabin vehicle. In Vehicle and Automotive Engineering 4: Select Proceedings of the 4th VAE2022, Miskolc, Hungary; Springer: Cham, Switzerland, 2022; pp. 357–373. [Google Scholar]
  51. Duffie, J.A.; Beckman, W.A. Solar Engineering of Thermal Processes; Wiley: New York, NY, USA, 1980. [Google Scholar]
Figure 1. The leapfrog–hopscotch (LH) structure.
Figure 1. The leapfrog–hopscotch (LH) structure.
Energies 16 04604 g001
Figure 2. (A) One-layer wall. (B) Two-layer wall with brick-insulator and a thermal bridge.
Figure 2. (A) One-layer wall. (B) Two-layer wall with brick-insulator and a thermal bridge.
Energies 16 04604 g002
Figure 3. Visualization and Arrangement of the generalized variables.
Figure 3. Visualization and Arrangement of the generalized variables.
Energies 16 04604 g003
Figure 4. The maximum errors as a function of the time step size for the tested methods in the case of the coarse mesh.
Figure 4. The maximum errors as a function of the time step size for the tested methods in the case of the coarse mesh.
Energies 16 04604 g004
Figure 5. The maximum errors as a function of the running time for the tested methods in the case of the coarse mesh.
Figure 5. The maximum errors as a function of the running time for the tested methods in the case of the coarse mesh.
Energies 16 04604 g005
Figure 6. The temperature distribution contour in Kelvin units for the measurement of the running time in the case of the coarse mesh (left) for ANSYS and MATLAB on (right).
Figure 6. The temperature distribution contour in Kelvin units for the measurement of the running time in the case of the coarse mesh (left) for ANSYS and MATLAB on (right).
Energies 16 04604 g006
Figure 7. The maximum errors as a function of the running time for the tested methods in the case of the medium mesh.
Figure 7. The maximum errors as a function of the running time for the tested methods in the case of the medium mesh.
Energies 16 04604 g007
Figure 8. The maximum errors as a function of the running time for the tested methods in the case of the fine mesh.
Figure 8. The maximum errors as a function of the running time for the tested methods in the case of the fine mesh.
Energies 16 04604 g008
Figure 9. (A) One-layer wall, (B) Two-layers wall with brick-insulator, (C) Two-layers with straight thermal bridge, (D) Two-layers with bent thermal bridge.
Figure 9. (A) One-layer wall, (B) Two-layers wall with brick-insulator, (C) Two-layers with straight thermal bridge, (D) Two-layers with bent thermal bridge.
Energies 16 04604 g009
Figure 10. Illustration of the points and grid arrangements that present the results for the different cases.
Figure 10. Illustration of the points and grid arrangements that present the results for the different cases.
Energies 16 04604 g010
Figure 11. The temperature distribution in Celsius units as a function of time in days for the long-term simulation of the one-layer wall.
Figure 11. The temperature distribution in Celsius units as a function of time in days for the long-term simulation of the one-layer wall.
Energies 16 04604 g011
Figure 12. The temperature distribution in Celsius units as a function of time in days for the long-term simulation of the two-layer wall.
Figure 12. The temperature distribution in Celsius units as a function of time in days for the long-term simulation of the two-layer wall.
Energies 16 04604 g012
Figure 13. The temperature distribution in Celsius as a function of time in days for the long-term simulation of the wall with straight thermal bridging.
Figure 13. The temperature distribution in Celsius as a function of time in days for the long-term simulation of the wall with straight thermal bridging.
Energies 16 04604 g013
Figure 14. The temperature distribution contour in Kelvin units for the long-term simulation in the case of the medium mesh case (left) for straight thermal bridge and bent thermal bridge (right) at the end of the last day.
Figure 14. The temperature distribution contour in Kelvin units for the long-term simulation in the case of the medium mesh case (left) for straight thermal bridge and bent thermal bridge (right) at the end of the last day.
Energies 16 04604 g014
Figure 15. The temperature distribution in Celsius units with days for the long-term wall simulation in case of exterior points of thermal bridging (straight and bent) compared with one-layer and two-layer cases.
Figure 15. The temperature distribution in Celsius units with days for the long-term wall simulation in case of exterior points of thermal bridging (straight and bent) compared with one-layer and two-layer cases.
Energies 16 04604 g015
Figure 16. Total heat loss distribution in Watt units as a function of time in days for all cases of the long-term wall simulation for all cases.
Figure 16. Total heat loss distribution in Watt units as a function of time in days for all cases of the long-term wall simulation for all cases.
Energies 16 04604 g016
Table 1. Names and short descriptions of the methods used.
Table 1. Names and short descriptions of the methods used.
AbbreviationName/Description of the Method
DF-DDufort–Frankel, where the Sigma and K terms are only in the denominator
DF-MDufort–Frankel, where the Sigma and K terms are present in a mixed way
DF-KDDufort–Frankel, where the K term turns up only in the denominator, and the Sigma term is present in a mixed way
DF-SDDufort–Frankel, where the Sigma term turns up only in the denominator, and the K term is present in a mixed way
OOEHOriginal odd–even hopscotch
NS-OEHOOEH with the non-standard treatment of convection and radiation
LH Pseudo-ImpLeapfrog–hopscotch scheme with pseudo-implicit treatment of the Sigma term
LH InsideLH with inside treatment of the Sigma term (it appears in the numerator)
LH MixedLH with a combination of the pseudo-implicit and inside treatment
ANSYS-MMixed solver (controlled by ANSYS)
ANSYS-DDirect solver
ANSYS-IIterative solver
ode15sA 1 to 5-order numerical differentiation formula with variable-step and variable order (VSVO), that was designed for stiff problems
ode23tApplies the trapezoidal rule with a free interpolant
ode23tbUses trapezoidal rule in the first stage and a backward differentiation formula in the second one
ode23Second (third) order Runge–Kutta–Bogacki–Shampine method
ode45A fourth (fifth) order Runge–Kutta–Dormand–Prince solver
ode1131-13 order VSVO Adams–Bashforth–Moulton solver
Table 2. The properties of the materials used.
Table 2. The properties of the materials used.
ρ k g m 3 c J k g 1 K 1 k W m 1 K 1
Brick19008400.73
Rigid Polyurethane Foam32014000.023
Steel Beam780084016.2
Table 3. The convection, radiation, and heat source parameters on both sides of the wall components in the case of the two-layer wall.
Table 3. The convection, radiation, and heat source parameters on both sides of the wall components in the case of the two-layer wall.
h c W m 2 K ε σ * W m 2 K 4 × 1 0 8 q * W / m 2
Right Elements220.84.5435.39
Left Elements90.95.1360.95
Table 4. The comparison between MATLAB and ANSYS at different mesh sizes and boundary conditions.
Table 4. The comparison between MATLAB and ANSYS at different mesh sizes and boundary conditions.
The Type of Boundary
Condition
Grid TypeThe Maximum Error between Analytical
Solution
MATLABANSYS
Zero Dirichletcoarse 9 . 114 × 10 6 3.577 × 10 5
medium 2 . 279 × 10 6 2.097 × 10 5
fine 1.013 × 10 6 1.863 × 10 5
Only convection on the surfacemedium 1 . 24 × 10 10 8.49 × 10 5
Only radiation on the surfacemedium 4.04 × 10 9 6.344 × 10 5
Table 5. The convection, radiation, and heat source parameters are on both sides of the wall for all the previously mentioned wall types [42].
Table 5. The convection, radiation, and heat source parameters are on both sides of the wall for all the previously mentioned wall types [42].
h c · W m 2 K 4 ε σ * · W m 2 K 4 × 1 0 8
Right Elements0.6–22.450.95.1
Left Elements90.73.97
Table 6. Heat loss through the 1 m2 piece of wall, and the energy cost in HUF and USD.
Table 6. Heat loss through the 1 m2 piece of wall, and the energy cost in HUF and USD.
One LayerTwo
Layers
Two Layers with a Straight BridgeTwo Layers with Bent Bridge
Heat loss (full month, kWh)19.141.995.295.01
The cost in HUF717.1974.63198.24187.8
The cost in USD1.90.20.530.5
Table 7. Heat loss through the wall and energy cost in HUF and USD for the coldest day.
Table 7. Heat loss through the wall and energy cost in HUF and USD for the coldest day.
One LayerTwo
Layers
Two Layers with a Straight BridgeTwo Layers with Bent Bridge
Heat loss (one day, kWh)0.7130.0760.2130.2
The cost in HUF26.72.857.967.56
The cost in USD0.0710.00760.0220.02
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Omle, I.; Askar, A.H.; Kovács, E.; Bolló, B. Comparison of the Performance of New and Traditional Numerical Methods for Long-Term Simulations of Heat Transfer in Walls with Thermal Bridges. Energies 2023, 16, 4604. https://0-doi-org.brum.beds.ac.uk/10.3390/en16124604

AMA Style

Omle I, Askar AH, Kovács E, Bolló B. Comparison of the Performance of New and Traditional Numerical Methods for Long-Term Simulations of Heat Transfer in Walls with Thermal Bridges. Energies. 2023; 16(12):4604. https://0-doi-org.brum.beds.ac.uk/10.3390/en16124604

Chicago/Turabian Style

Omle, Issa, Ali Habeeb Askar, Endre Kovács, and Betti Bolló. 2023. "Comparison of the Performance of New and Traditional Numerical Methods for Long-Term Simulations of Heat Transfer in Walls with Thermal Bridges" Energies 16, no. 12: 4604. https://0-doi-org.brum.beds.ac.uk/10.3390/en16124604

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop