Design of an Algorithm for Modeling Multiple Thermal Zones Using a Lumped-Parameter Model

Abstract

The generation of mathematical models for the analysis of buildings with multiple thermal zones is a large and complex task. Furthermore, the order and complexity of the dynamical model are increased by the number of included thermal zones. To overcome this problem, this paper presents an algorithm to define the mathematical model automatically, using the geometric and physics parameters as inputs. Additionally, the spatial position of each thermal zone must be recorded in an arrangement called a contact matrix. The algorithm for modeling systems with multiple thermal zones is the main contribution of this work. This algorithm is presented in pseudocode format and as an annex, an implementation in MATLAB software. One of the advantages of this methodology is that it allows us to work with parallelepipeds and not necessarily cubic thermal zones. The algorithm allows us to generate mathematical models with symbolic variables, starting from the knowledge of how many thermal zones compose the system and its geometric organization. This information must be organized in a matrix arrangement called a contact matrix. Different arrays of thermal zones were constructed with wooden boxes to verify the functionality of the models generated with the algorithm. Each case provided information that allowed us to adjust the mathematical models and their simulations, obtaining a range of errors between experimental and simulated temperatures from $2.08$ to $5.6$, depending on the number of thermal zones studied.

1. Introduction

Building modeling is an increasingly important field of research. In the last decade, several papers have been published in this field [1,2,3,4]. This increase has been driven in part by international agreements for environmental protection, such as the Kyoto Protocol [5] and the Paris Agreement [6]. Countries participating in this agreement have committed to formalizing various policies to reduce energy consumption and CO₂ emissions, especially in cities and buildings, which, according to additional research, can consume up to $40\%$ of the annual energy production [7,8,9].

The energy consumption in buildings depends on different factors, such as occupant activity [9] and heating, ventilation and air conditioning (HVAC) systems [10], among others. However, to understand the causes of energy consumption, buildings must be studied individually, taking into account the minimum spaces and energy requirements to ensure occupant comfort [11,12]. In this sense, mathematical models and simulators play an essential role, allowing researchers to analyze different configurations and situations in buildings, including huge interior spaces in buildings [13] or green roofs [14].

Regarding tools and strategies to generate and study the mathematical models of buildings, several works have been carried out [15,16,17,18]. Many works use specific simulators such as Energy Plus [19], Contam [20] and Trnsys [21] to understand and model problems such as choosing the best cooling system, airflow and energy savings [4,10,11]. These issues are considerably more complex if the building has multiple thermal zones. Papers such as in [22,23] mention the use of machine learning strategies to tackle these problems.

In different investigations, the lumped-parameter technique (LPM) is used to analyze thermal zones [7,13]. This technique allows a deeper numerical understanding of the thermal behavior in a closed enclosure, as it allows the user to generate a clear mathematical model, where each term can be explained and studied individually. These models are known as gray-box models and are the equivalent of the black-box models used by academic simulators that do not clearly show the mathematical model used for the simulation [24,25,26].

The LPM makes it possible to generate an equivalent electrical circuit to analyze a thermal zone and, using Kirchhoff’s laws, to define a mathematical model [27,28,29,30]. The order of the model depends on the number of resistors and capacitors used for each wall. However, the size and model requirements increase with the number of thermal zones studied. The need to consider multiple interacting thermal zones arises when studying more realistic buildings, such as residential houses or apartments [31,32,33,34].

Figure 1 shows different combinations of thermal zones, providing a complex heat- and air-transfer modeling problem. In this paper, an algorithm for modeling multiple thermal zones is presented. The algorithm can be extended to consider as many zones as necessary, requiring only the geometrical and physical parameters and information about the spatial position of each zone. The proposed methodology is not limited to thermal zones having a square geometry. This is a problem presented in similar algorithms [35]. In addition, the algorithm presented in [32] showed a statistical model for multiple thermal zones, but that type of algorithm uses a complex nomenclature and may be more difficult to use than the proposed methodology. Other works presented algorithms for modeling but were limited by the geometry, excessive information requirements and mathematical complexity [15,36,37].

The paper is organized as follows: Section 2 is devoted to describing the different organizations of the thermal zones, classified into cases. The description of the algorithm is given in Section 3. An experimental verification of the different mathematical models is presented in Section 4. Finally, conclusions and future work are presented in Section 6.

2. Mathematical Model

The lumped-parameter technique (LPM) represents the thermal zones with an RC circuit and defines a mathematical model using circuit analysis strategies. The size of the mathematical model can change depending on the structure used. The structure presented in Figure 2 was used in previous work, showing a good behavior and adaptability, and thus was used as equivalent for each thermal zone studied in this work [37,38,39].

The resistance and capacitors are calculated as a function of the geometric and physical parameters of the thermal zones, such as wall dimensions (L for thickness and A for surface area), density ($\rho$), thermal conductivity ($kt$) and specific heat ($Ce$). Equations (1)–(6) show the relationships for each surface. The subscripts i and j have a key function: i represents each wall in the thermal zone. Because of that, $i=1,2,\dots ,6$ in all the cases. Furthermore, the subscript j represents the number of analyzed thermal zones (m), so $j=1,2,\dots ,m$.

$$R_{i,j}=\frac{L_{i,j}}{2k{t}_{i,j}{A}_{i,j}}$$

(1)

$$R_{i{n}_j}=\frac{1}{h{i}_j{A}_{i,j}}$$

(2)

$$R_{e{x}_j}=\frac{1}{h{e}_j{A}_{i,j}}$$

(3)

$$R{w}_j=\frac{{\prod }_{k=1}^6(R_{k,j}+R_{i{n}_j})}{{\sum }_{i=1}^6{\prod }_{k=1}^6\frac{R_{k,j}+R_{i{n}_j}}{R_{i,j}+R_{i{n}_j}}}$$

(4)

$$C{w}_j=\sum _{i=1}^6{\rho }_{i}C{e}_i{A}_i{L}_i$$

(5)

$$C{r}_j={\rho }_{a}C{e}_a{V}_a$$

(6)

The resistance $R_{i,j}$ represents the thermal resistance of the exterior half-wall and can be extended to include heat transfer to the surroundings by adding Equation (2) for cases where the surface is exposed to outside air. The resistance $R{w}_j$ corresponds to the parallel reduction of the resistances of the inner half-walls. The constants $hi$ and $he$ are the heat transfer coefficients between the walls and the indoor and outdoor air, respectively. These values must be adjusted to each specific environmental situation [40,41]. Finally, $C{w}_j$ and $C{r}_j$ are used to represent the thermal capacity of the walls and internal air, respectively.

Traditionally, the process of developing a mathematical model using LPM consists of three stages:

• The identification of geometrical and physical parameters of the thermal zones, including internal loads.
• RC circuit construction and the calculation of resistors and capacitors
• Define the equations of the dynamic system, using the theory of energy flow transfer.

These three steps can be arduous if the system includes several thermal zones. In addition, this section represents the final two steps of the analysis of a thermal zone with different spaces.

2.1. Case 1: A Single Thermal Zone (m = 1)

To facilitate understanding of the notation, the same wall numbering is always used in the remainder of this paper. Figure 3 shows the chosen assignment, taking the front and back faces as walls one and two, the right and left walls are denoted by walls three and four. Meanwhile, surfaces five and six are the top and bottom surfaces, respectively.

The circuit for the single thermal zone case is presented in Figure 2. Using classical electrical circuit analysis strategies, Equations (7) and (8) were defined, where $Tw$ and T represent the wall and indoor-air temperature, respectively.

$$\begin{array}{c}C{w}_1\dot{T{w}_1}=\frac{T_1}{R{w}_1}-T{w}_1\left(\frac{1}{R_{1,1}}+\frac{1}{R_{2,1}}+\frac{1}{R_{3,1}}+\frac{1}{R_{4,1}}+\frac{1}{R_{5,1}}+\frac{1}{R_{6,1}}+\frac{1}{R{w}_1}\right)\hfill \\ \hfill +\left(\frac{1}{R_{1,1}}+\frac{1}{R_{2,1}}+\frac{1}{R_{3,1}}+\frac{1}{R_{4,1}}+\frac{1}{R_{5,1}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,1}}\end{array}$$

(7)

$$C{r}_1\dot{T_1}=-\frac{T_1}{R{w}_1}+\frac{T{w}_1}{R{w}_1}$$

(8)

The state variable for this case can be shortened to $X={\left[\begin{array}{cc}T{w}_1& T_1\end{array}\right]}^T$, the inputs to the system are the ambient temperature denoted by $T_{exterior}$ and the ground temperature denoted by $T_{ground}$; for the specific situation where the thermal zone is separated from the ground, the ground temperature must be substituted for the environmental temperature. In addition, Equation (8) can be extended to include thermal loads by adding the term $\pm Q$, taking Q as the power of the load, positive when the load is a heating load, and negative when the load is a cooling system. This case was studied in previous works and was not analyzed in depth in this paper [40,42].

2.2. Case 2: Two Thermal Zones (m = 2)

Figure 4 shows the spatial location of the thermal zones. In this configuration, surface 3 of zone 1 is in contact with surface 4 of zone 2, according to the numbering presented in Figure 3.

For the case with two thermal zones, the circuit presented in Figure 2 must be extended to include the second thermal zone; Figure 5 shows the new circuit, taking $T{w}_2$ and $T_2$ as the temperatures for the second zone. Using the circuit for two thermal zones it is possible to derive Equations (9) to (12).

$$\begin{array}{c}C{w}_1\dot{T{w}_1}=\frac{T_1}{R{w}_1}-\left(\frac{1}{R{w}_1}+\frac{1}{R_{1,1}}+\frac{1}{R_{2,1}}+\frac{1}{R_{4,1}}+\frac{1}{R_{5,1}}+\frac{1}{R_{6,1}}+\frac{1}{R_{3,1}+R_{4,2}}\right)T{w}_1\hfill \\ \hfill +\frac{T{w}_2}{R_{3,1}+R_{4,2}}+\left(\frac{1}{R_{1,1}}+\frac{1}{R_{2,1}}+\frac{1}{R_{4,1}}+\frac{1}{R_{5,1}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,1}}\end{array}$$

(9)

$$C{r}_1\dot{T_1}=-\frac{T_1}{R{w}_1}+\frac{T{w}_1}{R{w}_1}$$

(10)

$$\begin{array}{c}C{w}_2\dot{T{w}_2}=\frac{T_2}{R{w}_2}-\left(\frac{1}{R{w}_2}+\frac{1}{R_{1,2}}+\frac{1}{R_{2,2}}+\frac{1}{R_{3,2}}+\frac{1}{R_{5,2}}+\frac{1}{R_{6,2}}+\frac{1}{R_{4,2}+R_{3,1}}\right)T{w}_2\hfill \\ \hfill +\frac{T{w}_1}{R_{4,2}+R_{3,1}}+\left(\frac{1}{R_{1,2}}+\frac{1}{R_{2,2}}+\frac{1}{R_{3,2}}+\frac{1}{R_{5,2}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,2}}\end{array}$$

(11)

$$C{r}_2\dot{T_2}=-\frac{T_2}{R{w}_2}+\frac{T{w}_2}{R{w}_2}$$

(12)

For this case, the state variables are $X={\left[T{w}_1{T}_{1}T{w}_2{T}_2\right]}^T$ and the system inputs are the same as in case 1. Analyzing this mathematical model, it is easy to deduce the relationship between the size of the mathematical model (N) and the number of thermal zones. Using the relation $N=2m$, the size of the dynamic system can be known in future cases.

2.3. Case 3: Three Thermal Zones ($m=$ 3)

The organization of three thermal zones is presented in Figure 6. In this case, thermal zone 1 is next to two independent zones, zone 2 on the left and zone 3 on the right.

Following the same process described in the two previous cases, the circuit designed for this configuration is presented in Figure 7. Resistors three and four in zone 1 are dedicated to linking the first thermal zone to the adjacent zones. Six differential equations ($N=6$) are required to model this circuit, which are derived and presented in the following set of equations:

$$\begin{array}{c}C{w}_3\dot{T{w}_3}=\frac{T_3}{R{w}_3}-\left(\frac{1}{R{w}_3}+\frac{1}{R_{1,3}}+\frac{1}{R_{2,3}}+\frac{1}{R_{3,3}}+\frac{1}{R_{4,3}+R_{3,1}}+\frac{1}{R_{5,3}}+\frac{1}{R_{6,3}}\right)T{w}_3\hfill \\ \hfill +\frac{T{w}_1}{R_{4,3}+R_{3,1}}+\left(\frac{1}{R_{1,3}}+\frac{1}{R_{2,3}}+\frac{1}{R_{3,3}}+\frac{1}{R_{5,3}}\right)+\frac{T_{ground}}{R_{6,3}}\end{array}$$

(13)

$$C{r}_3\dot{T_3}=-\frac{T_3}{R{w}_3}+\frac{T{w}_3}{R{w}_3}$$

(14)

$$\begin{array}{c}C{w}_2\dot{T{w}_2}=\frac{T_2}{R{w}_2}-\left(\frac{1}{R{w}_2}+\frac{1}{R_{1,2}}+\frac{1}{R_{2,2}}+\frac{1}{R_{3,2}+R_{4,1}}+\frac{1}{R_{4,2}}+\frac{1}{R_{5,2}}+\frac{1}{R_{6,2}}\right)T{w}_2\hfill \\ \hfill +\frac{T{w}_1}{R_{3,2}+R_{4,1}}+\left(\frac{1}{R_{1,2}}+\frac{1}{R_{2,2}}+\frac{1}{R_{4,2}}+\frac{1}{R_{5,2}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,2}}\end{array}$$

(15)

$$C{r}_2\dot{T_2}=-\frac{T_2}{R{w}_2}+\frac{T{w}_2}{R{w}_2}$$

(16)

$$\begin{array}{c}C{w}_1\dot{T{w}_1}=\frac{T_1}{R{w}_1}+\frac{T{w}_2}{R_{4,1}+R_{3,2}}+\frac{T{w}_3}{R_{3,1}+R_{4,3}}+\left(\frac{1}{R_{1,1}}+\frac{1}{R_{2,1}}+\frac{1}{R_{5,1}}\right)T_{exterior}\hfill \\ \hfill -\left(\frac{1}{R{w}_1}+\frac{1}{R_{1,1}}+\frac{1}{R_{2,1}}+\frac{1}{R_{3,1}+R_{4,3}}+\frac{1}{R_{4,1}+R_{3,2}}+\frac{1}{R_{5,1}}+\frac{1}{R_{6,1}}\right)T{w}_1+\frac{T_{ground}}{R_{6,1}}\end{array}$$

(17)

$$C{r}_1\dot{T_1}=-\frac{T_1}{R{w}_1}+\frac{T{w}_1}{R{w}_1}$$

(18)

For this case, $T{w}_3$ and $T_3$ correspond to the temperatures of the third thermal zone. In this model, the central thermal zone is influenced by the adjacent zones. This means that each zone can be affected by a maximum number of 6 thermal zones in the current configuration.

2.4. Case 4: Four Thermal Zones ($m=$ 4)

The last case study is present in Figure 8; in this configuration, thermal zone 1 is influenced by the three adjacent thermal zones, using surfaces 1, 3 and 4 as contact areas and reducing the ambient temperature over the central zone.

The equivalent circuit for this configuration is shown in Figure 9. In this structure, it is evident that zones 2, 3 and 4 have no effect on each other. The only relationship is an indirect one, using zone 1 as a conductance.

The mathematical model established for this case is present in Equations (19)–(26); it should be noted that only Equation (25) includes the wall temperature of all zones, allowing this zone to be a driver of the indirect relationship between adjacent zones.

$$\begin{array}{c}C{w}_4\dot{T{w}_4}=\frac{T_4}{R{w}_4}-\left(\frac{1}{R{w}_4}+\frac{1}{R_{1,4}}+\frac{1}{R_{2,4}+R_{1,1}}+\frac{1}{R_{3,4}}+\frac{1}{R_{4,4}}+\frac{1}{R_{5,4}}+\frac{1}{R_{6,4}}\right)T{w}_4\hfill \\ \hfill +\frac{T{w}_1}{R_{2,4}+R_{1,1}}+\left(\frac{1}{R_{1,4}}+\frac{1}{R_{3,4}}+\frac{1}{R_{4,4}}+\frac{1}{R_{5,4}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,4}}\end{array}$$

(19)

$$C{r}_4\dot{T_4}=-\frac{T_4}{R{w}_4}+\frac{T{w}_4}{R{w}_4}$$

(20)

$$\begin{array}{c}C{w}_3\dot{T{w}_3}=\frac{T_3}{R{w}_3}-\left(\frac{1}{R{w}_3}+\frac{1}{R_{1,3}}+\frac{1}{R_{2,3}}+\frac{1}{R_{3,3}}+\frac{1}{R_{4,3}+R_{3,1}}+\frac{1}{R_{5,3}}+\frac{1}{R_{6,3}}\right)T{w}_3\hfill \\ \hfill +\frac{T{w}_1}{R_{4,3}+R_{3,1}}+\left(\frac{1}{R_{1,3}}+\frac{1}{R_{2,3}}+\frac{1}{R_{3,3}}+\frac{1}{R_{5,3}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,3}}\end{array}$$

(21)

$$C{r}_3\dot{T_3}=-\frac{T_3}{R{w}_3}+\frac{T{w}_3}{R{w}_3}$$

(22)

$$\begin{array}{c}C{w}_2\dot{T{w}_2}=\frac{T_2}{R{w}_2}-\left(\frac{1}{R{w}_2}+\frac{1}{R_{1,2}}+\frac{1}{R_{2,2}}+\frac{1}{R_{3,2}+R_{4,1}}+\frac{1}{R_{4,2}}+\frac{1}{R_{5,2}}+\frac{1}{R_{6,2}}\right)T{w}_2\hfill \\ \hfill +\frac{T{w}_1}{R_{3,2}+R_{4,1}}+\left(\frac{1}{R_{1,2}}+\frac{1}{R_{2,2}}+\frac{1}{R_{4,2}}+\frac{1}{R_{5,2}}\right)T_{exterior}+\frac{T_{ground}}{R_{6,2}}\end{array}$$

(23)

$$C{r}_2\dot{T_2}=-\frac{T_2}{R{w}_2}+\frac{T{w}_2}{R{w}_2}$$

(24)

$$\begin{array}{c}C{w}_1\dot{T{w}_1}=\frac{T_1}{R{w}_1}+\frac{T{w}_2}{R_{3,2}+R_{4,1}}+\frac{T{w}_3}{R_{3,1}+R_{4,3}}+\frac{T{w}_4}{R_{1,1}+R_{2,4}}+\left(\frac{1}{R_{2,1}}+\frac{1}{R_{5,1}}\right)T_{exterior}\hfill \\ -\left(\frac{1}{R{w}_1}+\frac{1}{R_{1,1}+R_{2,4}}+\frac{1}{R_{2,1}}+\frac{1}{R_{3,1}+R_{4,3}}+\frac{1}{R_{4,1}+R_{3,2}}+\frac{1}{R_{5,1}}+\frac{1}{R_{6,1}}\right)T{w}_1\\ \hfill +\frac{T_{ground}}{R_{6,2}}\end{array}$$

(25)

$$C{r}_1\dot{T_1}=-\frac{T_1}{R{w}_1}+\frac{T{w}_1}{R{w}_1}$$

(26)

The process described in this section is traditionally used to analyze single and multiple thermal zones [26,40]. In all cases, the model must be adjusted to reproduce the experimental measurements, which can be a complex task depending on the size of the model and the number of parameters considered. In this section, different numbers of thermal zones were presented, starting with a single $m=1$ zone up to $m=4$, moving from a $2\times 2$ to an $8\times 8$ system, increasing the tuning complexity and allowing for a more significant percentage error with respect to the experimental data.

3. Algorithm Design

The mathematical models presented in Section 2 and constructed manually require time and effort, but in the case of changes in the position of the thermal zones, steps 2 and 3 must be repeated.

A mathematical modeling algorithm was developed to avoid this problem, using the evidenced pattern. A pair of differential equations was added to the system for each thermal zone in all cases. The structure of these equations is presented in Equations (27) and (28). In this structure, the coefficients $\alpha$, $\beta$ and $\varphi$, which correspond to the impact of adjacent zones, the environment and the ground temperature in each zone, must be calculated.

$$C{w}_j\dot{T{w}_j}=\frac{T_j}{R{w}_j}+\sum _{z=1}^m{\alpha }_{z}T{w}_z+\beta T_{exterior}+\varphi T_{ground}$$

(27)

$$C{r}_j\dot{T_j}=-\frac{T_j}{R{w}_j}+\frac{T{w}_j}{R{w}_j}$$

(28)

To capture the spatial organization of the thermal zones, a matrix called a contact matrix was used. This matrix had the dimensions $m\times 6$; each row recorded the information of a different zone, thus the first row had information about zone 1, the second row about the second thermal zone, continuing up to the last thermal zone in row m.

The columns in the contact matrix maintained the surface numbering described in Figure 3, using column 1 to indicate whether there was a thermal zone adjacent to the front surface, column 2 to record the existence of a thermal zone on the rear surface and continuing this assignment for all six surfaces; in case there was no adjacent thermal zone, the elements in the column for these directions were replaced by zero.

$$u=\left[\begin{array}{cccccc}{u}_{1,1}& u_{1,2}& u_{1,3}& \cdots & & u_{1,6}\\ u_{2,1}& u_{2,2}& u_{2,3}& \cdots & & u_{2,6}\\ ⋮& & ⋮& \ddots & & ⋮\\ u_{m,1}& u_{m,2}& u_{m,3}& \cdots & & u_{m,6}\end{array}\right]$$

(29)

Equation (29) illustrates the structure of the contact matrix. Using the thermal zone configurations presented in Section 2, the contact matrices were built:

$$u_{m=1}=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right]$$

(30)

$$u_{m=2}=\left[\begin{array}{cccccc}0& 0& 0& 2& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right]$$

(31)

$$u_{m=3}=\left[\begin{array}{cccccc}0& 0& 3& 2& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right]$$

(32)

$$u_{m=4}=\left[\begin{array}{cccccc}4& 0& 3& 2& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0& 0\end{array}\right]$$

(33)

We had zone 4 adjacent to wall 1, so element $u(1,1)$ was 4; zone 3 was adjacent to the right surface in zone 1 and for this reason, $u(1,3)$ was 3; and applying the same logic, zone 2 was adjacent to the left surface in zone 1, so $u(1,4)=2$. The zeros represented the surfaces not connected to the thermal zones in the contact matrix. Rows 2, 3 and 4 contained information about zones 2, 3 and 4. In these cases, the only adjacent zone was zone 1; for this reason, most of the spaces were filled with zeros, and the position of zone 1 indicated the position of the first thermal zone with respect to the current zone.

Using the number of thermal zones m, the contact matrix u, the resistances $R_{i,j}$ and the geometrical and physical parameters as inputs, Algorithm 1 was used to determine the differential equations for modeling the multiple-thermal-zone system.

The objective of the algorithm was to generate a pair of equations for each thermal zone, so the $\beta$, $\varphi$ and $\alpha$ coefficients had to be determined for each pair of equations. However, they were initially assumed to be zero, and using the information from the contact matrix, they were recalculated.

The coefficient $\beta$ represented all surfaces in contact with the outside air and contained the sum of the corresponding resistances. The coefficient $\varphi$ was exclusively dedicated to regulate the heat transfer with the ground. Finally, the coefficients $\alpha$ summarized the influence of adjacent thermal zones on the current zone. However, there was a special case, the condition for $i=j$, where this coefficient summarized the impact of the thermal walls on the interior temperature. In case of walls without connection with other thermal zones, only their own resistance was considered; for this purpose, the $\delta$ function was used, using the following behavior:

$$\delta \left(x\right)=\left\{\begin{array}{c}1\mapsto x=0\\ 0\mapsto \forall x\ne 0\end{array}\right.$$

(34)

The constant $k_1$ grouped resistors for the walls with no connection to other zones, while the constant $k_2$ grouped resistors with adjacent zones.

Algorithm 1 Algorithm to design differential equation systems for multiple thermal zones

4. Experimental Development

To evaluate the behavior of the algorithm, a set of experimental tests was executed. Each test was designed to verify the different cases modeled in section two. These models were corroborated by the algorithm presented. The algorithm was programmed in MATLAB software and its results are presented in Appendix A.

The case study selected to evaluate the mathematical models was that of multiple thermal zones made with balsa wood. The physical and geometrical parameters are shown in

Table 1. Physical and geometrical parameters.

Geometrical Parameters
x	y	z	$L_d$	$L_w$
0.111 m	0.111 m	0.2 m	0.005 m	0.008 m
Physical Parameters
$\rho$	Ce		kt
$150.7485\frac{\mathrm{kg}}{{\mathrm{m}}^3}$	$3.0512\frac{\mathrm{kJ}}{\mathrm{kg}·°\mathrm{C}}$		$0.1882\frac{\mathrm{kJ}}{\mathrm{m}·\mathrm{h}·°\mathrm{C}}$

.

We used x, y and z as the horizontal length, depth and height of the thermal zone. $L_d$ and $L_w$ represented the door and wall thickness, respectively. The physical parameters were the density ($\rho$), specific heat ($Ce$) and thermal conductivity ($kt$) [43].

The experiment was run in various stages, increasing the number of thermal zones according to the case studies. In all cases, thermal zone 1 was equipped with an internal heating load. Figure 10 shows a 100 W incandescent bulb as the heat source.

Figure 11 shows the construction of multiple thermal zones according to cases 2, 3 and 4. In addition, each thermal zone had a DHT11 temperature sensor to record the internal temperature, and an additional sensor was dedicated to recording the ambient temperature.

The tests were conducted between 22 and 25 March 2021. Four databases were constructed, one for each case, in which charging processes were carried out with the internal heat source activated for 5 min then turned off and the thermal zones allowed to discharge, releasing the heat to the environment for 30 min.

Tuning Process

In order for the mathematical model to represent the phenomenon under study, a tuning process was carried out. The tuning process had two parts, the first one consisted of calculating fixed values, such as conduction resistances and thermal capacities. These parameters depended on fixed characteristics, such as geometrical and physical characteristics. The second part was dedicated to tuning the variable parameters, i.e., the heat transfer coefficients $he$ and $hi$. Each surface could calculate these coefficients, but in this work, it was assumed that all surfaces in a single thermal zone had the same coefficients $he$ and $hi$.

The pattern search algorithm was used to adjust the heat transfer coefficients, taking as objective function the error between the experimental data and the simulated results. The fitting process was performed for the loading and unloading phases independently, allowing a different set of coefficients for each phase and thermal zone.

Using the experimental database, the model was fitted for case 1. The results are presented in Figure 12, using a red line to represent the experimental data. The blue line represents the simulated results and finally, the green line is used to represent the ambient temperature. In that case, the error rate between experimental data and simulations was $2.0896\%$ over 3 h of processing.

Figure 13 shows the internal temperatures for two thermal zones, case 2. The blue and red lines represent the experimental and simulated temperature for zone 1. In addition, the green and black lines denote the experimental and simulated internal temperature in zone 2. In this case, the objective function used for the tuning process was the sum of the experimental and simulated results for each zone, i.e., $F_{objective}=e_1+e_2$, using $e_1$ as the percentage error for zone 1 and $e_2$ the error for zone 2. The errors obtained were 2.94% and $2.97\%$ for zone 1 and zone 2, respectively.

The results of the study on case 3 are presented in Figure 14. In that graph, the first row represents the internal temperature for zone 1, using the blue line for the simulation and the red line for the experiment, and the green line illustrates the environmental temperature. Similarly, the second and third rows represent the simulated and experimental temperatures for zones 2 and 3, respectively. For the case $m=3$, the objective function used was $F_{objective}={\sum }_{i=1}^m{e}_i$, considering the error rate simultaneously for the three zones and obtaining $e_1=3.97\%$, $e_2=3.7\%$ and $e_3=4.88\%$ for each thermal zone.

The last case studied ($m=4$) is presented in Figure 15. In that figure, the first row shows the temperature for thermal zone 1, the second, third and fourth rows are dedicated to thermal zones 2, 3 and 4, respectively. The blue line corresponds to the simulation and the red line to the experimental results. The green lines represent the ambient temperature. The errors between the experimental and simulated temperature were $e_1=4.37\%$, $e_2=5.6\%$, $e_3=4.86\%$ and $e_4=4.9\%$.

The mathematical model was tuned for the charging and discharging process.

Table 2. Heat transmission coefficients for the charging process.

	$\mathit{hi}(\frac{\mathbf{kJ}}{\mathbf{h}·\mathbf{m}·\mathbf{K}})$				$\mathit{he}(\frac{\mathbf{kJ}}{\mathbf{h}·\mathbf{m}·\mathbf{K}})$
Zone	1	2	3	4	1	2	3	4
Case 1	86.9	0	0	0	12.36	0	0	0
Case 2	34.88	0.0022	0	0	12.29	58.77	0	0
Case 3	179.54	0.0047	0	0	2.75	199.9	179.99	0
Case 4	37.44	0.0037	0.0047	0.0047	0.0001	60.95	61.70	103.14

shows the coefficients for the charging phase, explaining the coefficients used in each case. Similarly,

Table 3. Heat transmission coefficients for the discharging process.

	$\mathit{hi}(\frac{\mathbf{kJ}}{\mathbf{h}·\mathbf{m}·\mathbf{K}})$				$\mathit{he}(\frac{\mathbf{kJ}}{\mathbf{h}·\mathbf{m}·\mathbf{K}})$
Zone	1	2	3	4	1	2	3	4
Case 1	0.6576	0	0	0	0.0145	0	0	0
Case 2	0.0025	0.001	0	0	3.39	21.45	0	0
Case 3	0.0027	0.0018	0.00027	0	18.65	0.0008	0.8387	0
Case 4	0.0027	0.0008	0.0008	0.0008	0.0008	4.08	1.89	4.67

shows the heat transfer coefficients for the discharge process.

Figure 16 shows the comparison between the different errors calculated for each case. In all cases, the chosen objective function ($F_{objective}={\sum }_{i=1}^m{e}_i$) allowed us to simultaneously tune the multiple heat transfer coefficients seeking to reduce the error between simulation and experimental results in each case.

Considering all error percentages equally valuable implied that all errors were minimized without emphasizing any one. On the other hand, the error obtained, for example in thermal zone 1, which was reused in all cases, showed a higher error percentage. This problem could be solved by implementing a function for the fitting process.

5. Discussion

In real applications, most offices and houses are composed of multiple thermal zones connected by walls, windows or doors. This situation drives research to formulate and use multiple-thermal-zone models. This paper described an algorithm to generate mathematical models for systems with equal geometry and materials.

The model used to evaluate the mathematical model was an arrangement of wooden boxes in a laboratory space with limited sources of error. However, in real applications, the system can contain multiple heat sources with fixed and variable power, such as windows, doors and occupants. In these cases, heat sources can function as new sources of error that reduce the accuracy of the model and hinder the fitting process. In addition, in residential or larger buildings, the internal air temperature can vary due to internal and external phenomena. For example, the stratification produced by the natural organization of indoor air according to its temperature is especially evident in areas with a high roof. External conditions can produce a high-convection process on a single wall, resulting in a temperature imbalance in the thermal zone. This situation can be corrected by using IR systems or specific ventilation systems that regulate the thermal energy between the different thermal zones.

6. Conclusions

The study of buildings with multiple thermal zones involves the generation of large mathematical models. In this work, an algorithm for the automatic generation of mathematical models was proposed. Initially, a manual process of model generation was presented, considering four cases.

The paper presented the modeling of four systems with multiple thermal zones used as case studies. In each case the organization of the thermal zones, the equivalent electrical circuit and the set of differential equations that represented them were described. The calculation process for all elements of the equivalent electrical circuit was described, including fixed physical parameters such as density, conductivity and specific heat. Additionally, variable physical parameters such as convection coefficients were included, which depended on the conditions under which the experiment was recorded and therefore had to be adjusted for each case. The objective of this process was to demonstrate the complexity of analyzing systems with multiple thermal zones and then compare it with the use of the proposed algorithm.

The proposed algorithm was designed to use the growth pattern in the previously calculated differential equations and allowed us to directly formulate the set of differential equations without drawing the equivalent circuit and circuit analysis. The main input of the algorithm was the contact matrix, an array containing the geometric position of all individual zones according to the nomenclature defined in this article.

A set of experimental tests was used to calibrate the mathematical model generated by the proposed algorithm. These experiments were based on thermal zones with balsa wood for the walls and an internal heat source to ensure energy transfer between the different thermal zones and the environment. The experiments were conducted indoors to reduce the environmental impact. The four systems with multiple thermal zones studied previously were experimentally represented. In each case, the individual zones were equipped with temperature sensors.

The tuning process for the heat transfer coefficient of all thermal zones was carried out using the error temperature, considering the difference between the simulation and experimental data. The best error results was about $2.08\%$ for a single thermal zone. Meanwhile, the maximal error was $5.6\%$ obtained in the case of four thermal zones. That error increase can be a problem with a bigger number of zones but can be solved using a better objective function in the tuning process. It is expected to be possible to improve the tuning process in future works and obtain small errors even if large numbers of thermal zones are used.

The internal and external convection coefficient were adjusted in each case. The experiment with a single thermal zone obtained the best error rate between the experimental data and the simulations, obtaining a $2.0896\%$ error. The second experiment with two thermal zones yielded an average error of $2.95\%$. The third experiment showed errors in each thermal zone of $3.97\%$, $3.7\%$ and $4.88\%$. Finally, the experiment with four zones showed larger errors in each individual zone, reaching an error of $5.6\%$.

With the appropriate fitting process, the mathematical models generated with the proposed algorithm can represent any combination of thermal zones, providing the information is in the format of contact matrices. For the case of systems consisting of parallelepipeds, the proposed algorithm can be a useful tool.