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

: 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.


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 onedimensional 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 diffusionor 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 AN-SYS' 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. Sections 3.1-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.

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: 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 are the specific heat, heat conductivity, and the density of the material, respectively.
Newton's law of cooling suggests that a term   a K u u  can describe convective heat transfer [42], where the ambient temperature a u (measured in Kelvin) can be considered as independent from the unknown variable u. Consequently, it makes sense to include the a Ku term in the expression q, which represents the heat source. Stefan-Boltzmann's law [43] says that a term 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 a Ku 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: Note that all terms in Equation (2) which is second order in x  , where 1 i ,..., 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: 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 In this case, the heat conduction equation can be discretized as follows 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 , 1 i i k  is the heat conductivity between node i and its right neighbor. The discretized equation attains the following form The dimensions of a cell, measured along its length and across its (typical) crosssection, are represented as x  and S . Since the volume of the cell can be given as 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): 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.

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 , and it is discretized as usual: 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.

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 where If we adapt it for Equation (7), we obtain where , , 1 , and , 1, ..., , 0, ..., 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 One should note that the sum in Equation (7) is separated into two terms: the first contains the actual cell variable, whereas the other, n i A , 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 which can be rearranged as 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 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 We exemplify these by presenting the "zeroth" stage equations as follows.
1. Pseudo-implicit treatment: 0   for the radiation term, which yields: 2. "Inside" treatment: 1   for the radiation term, which means that it is taken into account explicitly, which yields 3. Mixed treatment with an equal share of the previous two treatments, that is, for the radiation term, which yields

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, 1 i u has to be derived from 0 i u using some other approach. Here, we use the UPFD calculation method [43] for this purpose: 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: 5. DF-M: the Sigma and K terms in a mixed way: 6. DF-KD: the K term appears only in the denominator, and the Sigma term in a mixed way: 7. DF-SD: the Sigma term appears only in the denominator, and the K term in a mixed way: 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: Second stage: 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: The second stage is the same as in Equation (22).

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 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 1 T o l = 1 0  until an extremely small minimum value, which is 1 2 T o l = 1 0  . 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. Applies the trapezoidal rule with a free interpolant ode23tb Uses trapezoidal rule in the first stage and a backward differentiation formula in the second one ode23 Second (third) order Runge-Kutta-Bogacki-Shampine method ode45 A fourth (fifth) order Runge-Kutta-Dormand-Prince solver ode113 1-13 order VSVO Adams-Bashforth-Moulton solver

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.

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 2 A, 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  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 i properties are different, we can write it as follows

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 fin 20, 000 s t  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: 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 a u , 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: We assumed that the left and right elements have the following heat sources as follows: In terms of the left-hand side: In terms of the right-hand side: .
The constant initial temperature is prescribed:   , , 0 2 0 K 9 u x z t   .

Verification Using Analytical Solutions
For verification, we used three different cases. In the first one, we considered a onelayer 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: Zero Dirichlet boundary conditions are used:

One can quickly verify that the problem's analytical solution is
where the wave numbers are fixed to 1, 1 x z k k   , and substituting the physical properties of the brick, we obtained α. We used the analytical solution Equation (23) where the initial condition here is a constant temperature that equals to 290 K.

 Third verification (for radiation)
We have the ODE: 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 xcoordinate 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.

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.

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).

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.

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 100 t   s, where the total time was 31 24 3600 2 678 400 T , ,     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 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 0.0057 x   , whereas in the other cases, the widths of the cells are larger, 0.0076 x   , 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.

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=2678400 s, which indicates the simulation of a full winter month. We obtain the values of the coefficients in our equations as follows: [42], 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 onelayer, it yields a linear function of the x variable for the initial condition: 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:   mid For the insulator part:

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. Figure 10. Illustration of the points and grid arrangements that present the results for the different cases.

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. Figure 11. The temperature distribution in Celsius units as a function of time in days for the longterm simulation of the one-layer wall.

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).  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. 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.

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 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 twolayer cases. 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.

Calculation of Heat Loss
The total energy requirement Q in kWh units was calculated from Equation (25) cell, 1 1 ( ) where * cell cell Q q S   and * q is the current density of the heat loss, which can be calcu- 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.

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 60 t   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.

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 longterm 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, twolayer 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 m 2 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://www.mdpi.com/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.

Conflicts of Interest:
The authors declare no conflict of interest.

Symbols
Greek