Modification of building energy simulation tool TRNSYS for modelling nonlinear heat and moisture transfer phenomena by TRNSYS/MATLAB integration

Software for numerical simulation of various types of energy used in buildings, i.e. building energy simulation (BES), have become an essential tool for recent research pertaining to building physics. TRNSYS is a well-known BES used in both academia and the construction industry for a wide range of simulations, such as the design and performance evaluation of buildings and related facilities for heating, cooling, and ventilation. TRNSYS has a modular structure comprising various components, and each component is interconnected and compiled through a common interface using a FORTRAN compiler. Its modular structure enables interactions with various external numerical simulation tools, such as MATLAB, Python, and ESP-r. For ordinary simulations of building energy load using TRNSYS, the generic module Type 56 is usually recommended, which provides detailed physics modelling of building thermal behaviours based on unsteady energy conservation equations and Fourier’s law for each building material. However, Type 56 explicitly depends on the transfer function method to discretise the original differential equations; therefore, it cannot model nonlinear phenomena, such as latent heat and moisture transfer between a building surface and ambient air. In other words, the current TRNSYS cannot be used to estimate the effectiveness of evaporation during cooling, which is a typical passive design method. Hence, the authors developed a MATLAB/TRNSYS integration scheme, in which TRNSYS was modified to model simultaneous heat and moisture transfer from the wet roof surface of a building. This scheme enabled TRNSYS to calculate the rate of evaporative heat and moisture transfer dynamically from the roof surface, assuming a control volume approximation of the roof surface. Finally, the effect of evaporative cooling on the thermal performance of an Indian building was estimated using the modified model.


Introduction
In recent years, several building energy simulation (BES) tools have become popular for evaluating the energy performance of new and existing buildings along with the increasing importance of energy saving and emission reduction in the building sector. In general, BES describes building thermal behaviours, enables the estimation of heating and cooling loads, indoor thermal climate, solar gains, and air infiltration rate; hence, it can be used for the design of heating, ventilation, and airconditioning (HVAC) systems [1][2][3]. TRNSYS is one of these BES tools, which is well known for its numerous validation studies and ability to model buildings incorporated with other systems, such as HVAC and renewable energy sources [4,5]. It has a modular structure that enables it to interact with various external numerical simulation tools, such as MATLAB, Python, and ESP-r [6].
TRNSYS employs the conduction transfer function (CTF) method developed by Mitalas and Stephenson to calculate transient heat conduction through slabs [7]. It involves the determination of transfer functions for one-dimensional heat transfer through building components by solving conduction equations using Laplace and Z transform theories [8]. However, the current implementation of the CTF in TRNSYS multizone building models, known as Type 56, cannot simultaneously analyse heat and moisture transport through building components. Type 56 has no specified module to handle the moisture content inside building walls or roofs.
Nevertheless, moisture transport through building envelopes is important for building hygrothermal performance, which affects mould-growth-related problems. Furthermore, the appropriate modelling of moisture transfer and the dynamic variation of moisture content of building envelopes are essential for the analysis of passive evaporative cooling.
Evaporative cooling by spraying water over roof surface is a well-known passive cooling technique widely used to reduce cooling loads in hot and dry climates in tropical and subtropical regions [9]. Various experimental and numerical studies have been performed to evaluate the effectiveness of roof spray evaporative cooling systems [10][11]. Furthermore, combined heat and moisture transfer models have been proposed to quantify the effect of evaporation from wet roof surfaces [12][13]. However, the integration of these models with building energy simulation tools is rare owing to various difficulties, such as complexity of existing software packages, difference in time step to model heat and moisture transport, nonlinearity in modelling, and numerical convergence.
Purohit et al. [14] were the first to model roof surface evaporative cooling using TRNSYS. However, simplification was adopted in their modelling, in which the saturation vapour pressure of water was assumed as a linear function of temperature within a narrow temperature range. In addition, the effective heat and transfer coefficients and effective temperatures (combining both convective and evaporative heat transfer) were calculated separately and fed to TRNSYS as the underside heat transfer coefficient and roof temperature, respectively. Although this approach provides a method to analyse evaporative cooling using TRNSYS, it does not consider the effect of moisture dynamics of the roof materials. Furthermore, validation based on experimental or numerical results has not been reported. Spanaky et al. [15] investigated the passive cooling effect of a ventilated roof pond protected with a reflecting layer and presented its numerical modelling using TRNSYS. They utilised module Type 344b in TRNSYS, which was originally designed for outdoor swimming pools, and modified it to simulate a roof pond using Type 67. Their modelling presented good agreement with experimental data However, this approach cannot be directly applied for the modelling of evaporation from a wet roof or a roof with intermittent water spray.
Hence, the authors developed a component that enables the hygrothermal behaviour of a wet roof to be simulated by communicating with module Type 56 and the interaction with the entire simulation of a building model in TRNSYS. Herein, the roof surface evaporative cooling model is first presented briefly. In the subsequent section, the methodology that enables this model to interact with Type 56 is discussed with the associated complexities and assumptions. Next, the model was tested for a building with dry roof conditions, and the results were compared with the simulation results obtained for the same building using TRNSYS. Finally, the evaporative cooling effect reproduced by this model under the climatic condition of New Delhi, India are presented.

2.
Modelling of Roof surface evaporative cooling Previously, the authors have developed a simplified numerical model to estimate the cooling effect produced by transient evaporation from asphalt sheet roofing after a water spray considering the time variation of moisture content within a thin surface layer of roofs [17]. The mechanism can be explained as follows. The evaporation rate can be expressed as a function of water content of a roof surface layer. The water content is dependent on the amount of sprayed water and the subsequent evaporation. Hence, a water balance equation can be established for the surface layer and coupled dynamically with the heat conductive equation. The water balance equation and the evaporation rate are expressed by Equations 1 and 2, respectively. (1) The heat transfer through the roof after a water spray is modelled using a one-dimensional unsteady heat transfer equation [18], as shown in Equation 3. The conduction, convection, and evaporative heat loss from the top surface of the roof are described by Equations 4-9. (3) The convective heat transfer between the indoor surface of the roof and the room air can be expressed as show in Equation 10. The convective heat transfer coefficient of the indoor surface of a roof was maintained constant at 9.6 W/m 2 K. The outside convective heat transfer coefficient was calculated using a regression equation that was developed previously by one of the authors [19] and shown as Equation 11.
More detailed information regarding the numerical model and the associated parameters are available in a previous study by the authors [17].
TRNSYS is the acronym for transient system simulation, a simulation programme primarily used in renewable energy and building simulation for passive and active solar designs [20]. It comprises a graphical user interface called Simulation Studio, which is a simulation kernel, and various component modules called 'Type'. A simulation project was formulated by interconnecting a set of component modules, such as solar collector, pump, and multizone building in a definite manner to represent the actual physical phenomena and then saved as a project definition file or 'deck' file. Subsequently, the simulation kernel accessed this deck file, performed an unsteady simulation, and created the desired outputs. The simulation was performed iteratively by the built-in solvers embedded in the kernel using a fixed time step until convergence was attained.
The modularity of TRNSYS provides flexibility for adding external components modelled with different tools, such as MS Excel, Python, R, and MATLAB. The Type 155 module establishes a link between MATLAB and TRNSYS. It enables TRNSYS to communicate with the MATLAB Engine through a component object model interface. Type 155 has two different calling modes: (i) standard iterative component, and (ii) realtime controller. In the first mode, MATLAB is invoked at each call of each time step, whereas in the second mode, it is referred at the first call of each time step with the converged inputs from the previous time step. In the scheme we developed, this feature was used to model heat and moisture transfer from a wet roof using TRNSYS. The evaporation from a wet roof surface was modelled in MATLAB using the previously described numerical model and integrated with multizone building component Type 56 using Type 155.

Implementation in TRNSYS
Although Type 155 enables to external MATLAB components to be called from TRNSYS, many problems are involved in the design principle of the integrated environment, such as the compatibility of the MATLAB component, time step, and parameter settings.
In this study, we used a procedure similar to those of Type 399 of Nostand [21] and 'slab on grade' type [22] in the TESS library. Figs. 1-3 illustrate the principle of integration between the wet roof component in MATLAB and the multizone building in TRNSYS. It can be summarised as follows: i. The wet roof was modelled in MATLAB with all the required thermophysical parameters.
ii. The inputs to the wet roof component were as follows: solar heat flux (short wave radiation), wind velocity, atmospheric pressure, relative humidity, dry bulb temperature, and amount of water spray.
iii. The desired building to assess the effect of evaporative cooling was designed with Type 56.
iv. The roof of the building in Type 56 was designed as a dummy roof with high thermal conductivity and low back side heat transfer coefficient (Hback < 0.001) [22] to ensure a direct contact with the wet roof designed in MATLAB .
v. The temperature of the underside surface of a wet roof calculated using MATLAB was used as the boundary temperature for the dummy roof.
vi. The indoor air temperature and the heat gain of shortwave radiation absorbed through the wall of Type 56 building was used as an input to the MATLAB wet roof.

Numerical Method
The governing equations of heat and moisture transfer through a wet roof expressed by Equations 1 to 11 cannot be solved algebraically because of their nonlinearity. Hence, the control volume method was used for space discretisation and the finite difference method time discretisation. It is noteworthy that, to integrate the wet roof modelled in MATLAB with Type 56 of TRNSYS, the time step of both the simulations must be same. However, suitable time steps must be selected by considering the following: v. The lower the time base, the more coefficients are required to generate the transfer function equations. Therefore, a practical limitation is imposed on the selection of the lowest possible value for the time base for the problem under consideration.
vi. As a requirement, the Type 56 time base must be an integral multiple of the simulation time step.
For the given problem, the lowest possible value of the Type 56 wall time base is 1 h, below which TRNSYS cannot generate transfer function coefficients. Hence, the time base was maintained as 1 h, and the simulation time step was selected as 36 s or 0.01 h (ratio of time base to time step was 100) after several trial and error. Exterior

Performance test of developed TRNSYS-MATLAB integration for dry roof condition
The performance of the developed integrated module was tested by comparing the original TRNSYS model with the default Type 56 for a dry roof condition.
The simulation was performed for a single-room building. The building structure and thermophysical properties of the building materials are described in Table 1. Typical Meteorological Year (TMY) weather data of New Delhi, India (28.6139° N, 77.2090° E) for 6 d from 1 st to 6 th May were used for the boundary conditions, and the simulation results of the period from 3 rd to 6 th May were used for analysis to avoid the effects of the initial conditions. Fig. 4 depicts the hourly variations of solar radiation, relative humidity, ambient temperature, and wind speed. As the weather data were available at 1 h intervals, linear interpolation was employed to generate weather changes at each 36 s interval and used as input for the numerical simulations. The combined convective and radiative heat transfer coefficient for the inside surface of the building was assumed to be 9.6 W/m 2 K. The target building was assumed as ventilated with an air infiltration rate of 0.6 /h and no internal heat gain/sink.  Finally, the indoor air temperature of both models are compared in the top plot. The difference between the TRNSYS and TRNSYS-MATLAB integration models varied in a nearly cyclic manner from -1.2 °C to 2.1 °C. The mean absolute deviations in temperature of the room air, roof top, and bottom surface were calculated to be 0.3 °C, 0.7 °C, and 0.5 °C, respectively. This analysis proves that the new integrated model equipped with a roof modelled in MATLAB can successfully predict the thermal behaviour of a dry roof with minimal error. The discrepancy may be attributed to the lower time step selected for TRNSYS simulation, assumed values of material properties, and convective heat transfer coefficient. Therefore, this model can be used to analyse the heat transfer through a wet roof, which is otherwise impossible in TRNSYS.

Case study of roof surface evaporative cooling using modified TRNSYS tool.
To assess the effect of evaporative cooling by spraying water over the roof, a single-room building with a roof structure, as shown in Fig. 6, was selected as a case study. It was assumed that water was sprayed twice daily at 7:00 am and 12:00 pm at a uniform rate of 1/6 kg/min for 30 min. Considering the physical structure of the roof, the following assumptions were made for the numerical modelling: (i) no internal heat is generated within any material layer, (ii) each layer is homogenous with constant thermophysical properties, (iii) an appropriate contact exists between layers with no interface resistance, (iv) the thermophysical properties of the roof material is independent of water content. The modified TRNSYS model was used to analyse the temperature distribution through the roof with and without water spray. The simulation was performed using TMY weather data of New Delhi, India with a time step of 36 s. New Delhi is categorised as a composite climate zone that is extremely hot in the summer and severely cold in the winter. Generally, June is the hottest month, with the daily maximum temperature exceeding 40 °C. Therefore, a hot and clear week in June (6th to 13th) with a daily maximum temperature exceeding 35° and a relative humidity between 20% and 60% was selected for the analysis. Fig.  7 describes the weather data of that week. For a simple analysis, the simulation result for 8th June is discussed. Fig. 6. Schematic representation of the wet roof selected for the case study The effects of water spray on the indoor air temperature, moisture content of roof, and evaporative heat flux on 8 th June are depicted in Fig. 8. As shown from the graph depicting evaporative flux, evaporation begins immediately after the water spray at 7:00 am, reaches the maximum and then continues up to 10:00 am at a decreasing rate before diminishing to zero. Similar pattern of evaporation is seen after water spray at 12pm. But, the amplitude of the maximum evaporative flux at 12 pm is nearly four times that at 7 am. This is owing to the larger solar radiation and ambient temperature at noon compared with those in the morning. Furthermore, the effect of evaporation starting at 7:00 am remained for more than 2 h, whereas that starting from 12:00 pm diminished in 1.5 h. Owing to the high evaporation rate during noon time, the moisture content of the roof material decreased rapidly, resulting in a faster decrease in the evaporation flux. The maximum amount of evaporative flux was 2156 W/m 2 at 12:01 pm immediately after water spray.
From the plot of roof moisture content, it can be deduced that when water was sprayed over the roof, the moisture content increased rapidly and attained the maximum threshold value, which was 0.45 kg/m 2 for the roof material calculated in a previous study by the author [17]. Then as evaporation begins, the moisture content decreases and reduces to zero after a few hours. Comparing the plot of evaporative flux and moisture content, it can be inferred that, as the evaporation rate at 7:00 am is lower than that at 12 pm, the moisture content requires more time to diminish completely at 7:00 am than at a later time. The top plot shows the variation in outside surface temperature of the roof with and without employing evaporation cooling. As shown, before water was sprayed, the outside surface temperature of the roof for both cases were similar. However, as water was sprayed at 7 am, the corresponding roof temperature decreased. The decrease in roof temperature continued for approximately 2 h until the effect of evaporation cooling remained prevalent. When evaporation ceased, it again increased and immediately before water was sprayed again at 12:00 pm, it was approximately the temperature of the roof without employing evaporation. After the second water spray, the temperature of the outside roof surface decreased rapidly, and the difference between two roofs attained a maximum value of approximately 15 °C at 12:32 pm. As the effect of evaporation decreased, the temperature of the wet roof increased gradually and approached that of the dry roof after approximately 6 h. Furthermore, the second water spray at 12:00 pm exerted a more pronounced effect on the roof top surface temperature in terms of maximum temperature drop and continuation of cooling effect than the first spray at 7:00 am. This was because the rate of evaporation was higher at noon time than that in the morning. Therefore, it can be concluded that the time of water spray significantly affects the effectiveness of evaporative cooling. The variation in indoor air temperature for both cases are shown in the second plot along with the ambient temperature. Unlike the roof top surface temperature, the difference in room air temperature of the building with and without evaporation remained nearly uniform. The maximum difference between the two was approximately 1.4 °C at 4 pm. This was because the room air temperature was strongly affected by the assumption of the roof-area-to-room-volume ratio and the thermal mass of the building roof and walls.

Discussion and Conclusions
The water content of building materials are important for the evaporative cooling of buildings. Hence, the dynamic variation of moisture content of roof materials must be considered for the numerical modelling of evaporation from the roof surface after water spray. However, the current version of TRNSYS cannot simultaneously model heat and moisture transport through building envelopes. Therefore, a modified tool was developed by integrating TRNSYS with MATLAB using the Type 155 module.
The integration of MATLAB with TRNSYS involves defining a dummy roof in TRNSYS with a low backside heat transfer coefficient such that a direct contact is established between the building designed with Type 56 and the wet roof modelled in MATLAB. The room temperature obtained from TRNSYS simulation was provided as an input to the MATLAB simulation, whereas the temperature of the inside surface of the wet roof modelled in MATLAB was considered as the temperature of the dummy roof of the building designed using TRNSYS. This strategy is consistent with the 'slab on grade' type [22] in the TESS library.
The numerical solution of simultaneous heat and moisture transfer equations by MATLAB requires a small-time step in the order of a few seconds. However, a small-time step requires a smaller time base, which presents some practical limitations in terms of the generation of wall transfer function coefficients. Hence, considering all the constraints, the time base was assumed to be 1 h and the time step was fixed at 36 s.
The performance of the modified model was tested by comparing it with the original TRNSYS model for a dry roof condition without evaporation. The results of the modified model agreed well with that of TRNSYS with mean absolute deviations in temperature of only 0.3 °C, 0.7 °C, and 0.5 °C for the room air, roof top, and roof bottom surface, respectively.
A single-room building located in the climatic condition of New Delhi, India was selected as a case study to analyse the effect of evaporative cooling by spraying water over a roof surface at a rate of 1/6 kg/min for 30 min twice daily. The maximum decrease in external roof surface temperature owing to water spray at 12:00 pm on 8 th June was approximately 15 °C. Even though this model can be used successfully to analyse the evaporation from the building roof after water spray, its accuracy still depends highly on the appropriate selection of time base and time step of the simulation. A lower time base value (less than one hour) will improve the accuracy, but an improved algorithm for CTF generation in Type 56 is required.

Cp
Specific heat [J kg -