TRNSYS Simulation and Experimental Validation of Internal Temperature and Heating Demand in a Glass Greenhouse

The energy demand in greenhouses is enormous, and high-performance covering materials and thermal screens with varying radiometric properties are used to optimise the energy demand in building energy simulations (BES). Transient System Simulation (TRNSYS) software is a common BES tool used to model the thermal performance of buildings. The calculation of the greenhouse internal temperature and heating demand in TRNSYS involves the solution of the transient heat transfer processes. This study modelled the temperature and heating demand of two multi-span glass greenhouses with concave (farm A) and convex (farm B) shapes. This study aims to investigate the influence of the different BES longwave radiation modes on greenhouse internal temperature in different zones and the heating demand of a conditioned zone. The standard hourly simulation results were compared with the experimental data. The results showed that the standard and detailed modes accurately predicted greenhouse internal temperature (the Nash–Sutcliffe efficiency coefficient (NSE) > 0.7 for all three zones separated by thermal screens) and heating demand (NSE > 0.8) for farms A and B. The monthly heating demand predicted by the simple and standard radiation modes for farm A matched the experimental measurements with deviations within 27.7% and 7.6%, respectively. The monthly heating demand predicted by the simple, standard, and detailed radiation modes for farm B were similar to the experimental measurements with deviations within 10.5%, 6.7%, and 2.9%, respectively. In the order of decreasing accuracy, the results showed that the preferred radiation modes for the heating demand were standard and simple for farm A, and detailed, standard, and simple for farm B.


Introduction
Greenhouse farming aims to grow warm-season crops in winter and has been utilised to solve food shortages year-round. The control of indoor air, especially under adverse weather conditions, is crucial for crop production in greenhouses [1][2][3]. The energy consumption for the indoor air control contributes to 50% of the total cost of greenhouse production in many countries [4], making the greenhouse sector one of the most energy-intensive agricultural sectors, with significant economic and environmental impacts [5,6].
Passive energy-saving strategies reduce greenhouse energy consumption. This involves selecting high-performance greenhouse-covering materials and thermal screens to limit the amount of solar radiation the crop receives in summer and the rate of heat loss in winter [7]. Using thermal screens reduces the heat loss rate by 23-24% [8], and they are placed either internally or externally to ensure a favourable microclimate for the crop [9]. Most growers prefer internal thermal screens use owing to the increased heat retention [10] and reduced surface-to-volume ratio of the heating zone, leading to a reduction in the heating cost [11].
Several studies have used numerical modelling to assess the impact of different covering materials and thermal screens on greenhouse microclimate and energy consumption [12][13][14]. Zhang et al. [15] solved unsteady state equations using MATLAB to study the thermal performance of a greenhouse by altering the cover transmittance and absorbance for beam, diffuse, and ground-reflected radiations in a glass greenhouse. The surface radiative properties change with time and impact energy demand. Another study [16] investigated the effect of different plastic materials with spectral radiative properties ranging from 0.2 to 28 μm. The annual heating demand of the greenhouse was reduced by 6.3% using appropriate greenhouse materials in the developed model.
EnergyPlus and TRNSYS are building energy simulation (BES) tools for modelling the thermal performance of different types of buildings [17]. Within the simulation environment of these tools, there exist pre-existing components to construct a model [18]. The users have the flexibility of varying different parameters on each component to suit the application the model is intended for. The major setback of these tools for greenhouse modelling is their inability to account for the influence of crop and their environmental control systems [19,20]. Unlike EnergyPlus, TRNSYS permits users to develop new components and link with other simulation tools [7]. This makes TRNSYS, which was primarily developed for residential and commercial buildings, able to be continuously improved for greenhouse modelling. The Gembloux Dynamic Greenhouse Climate Model based on heat and mass balance of greenhouse layers was linked with TRNSYS as a sub-model for greenhouse air renewal rate and canopy resistance in a naturally ventilated greenhouse [21]. The model accurately predicted the greenhouse microclimate parameters, but there was an overestimation of the ventilation rate, mostly during night time. Ahamed et al. [22] also developed a TRNSYS model to predict the heating demand in a Chinese-style solar greenhouse (CSG). The results were compared with a heating simulation model (CSGHEAT) from MATLAB. The heating load from TRNSYS resulted in significant errors from assumptions such as constant infiltration rate, fixed scheduled for thermal screens and moisture gain.
A solar-radiation-based control for thermal screens as operated in the CSGHEAT model was implemented in TRNSYS using the equation editor in the simulation studio [23]. Based on the control, the multi-layer thermal screens reduced greenhouse energy consumption compared with double and single thermal screens investigated. The thermal screens also ensured a temperature difference of approximately 5 °C between the three zones in the multi-span greenhouse. However, with a unidirectional heat transfer model, there was no significant temperature difference within four zones of the greenhouse with moveable thermal screens on the north and south walls [24]. Asa'd et al. [25] performed a parametric analysis using TRNSYS to study the effect of the cooling and heating setpoint temperatures, cover materials, rock-bed air flow rate, and mechanical ventilation to improve the overall design of a greenhouse. The study concluded that the U-value of the materials and the rock-bed airflow rate have the most significant effect on the greenhouse internal temperature. Rabiu et al. [26] used TRNSYS and hotbox to evaluate the greenhouse energy-saving screen properties. The investigation also concluded that the screens' infiltration has the most significant impact on their U-values.

Scope of the Study
BES tools are used to investigate the effect of the different parameters on greenhouse microclimate and predict the most favourable operating conditions that minimise the energy consumption. With TRNSYS, several studies have reported the effect of the regime type, environmental control systems, covering materials, greenhouse shape, and orientation on the thermal performance of the greenhouse, with no consideration to the radiation modes [27][28][29]. The accurate modelling of the solar heat gains in the greenhouse with the different radiation modes in TRNSYS is also an important aspect that should be considered. To this end, a multi-span glass greenhouse area of 3942 m 2 was divided in two (farm A, 2160 m 2 , and farm B, 1782 m 2 ) based on their geometry to study the effect of the radiation modes on the internal temperature and heating demand of the greenhouse. The resulting TRNSYS simulation results were validated with the experimental measurements of farm A (concave-shaped) and farm B (convex-shaped) greenhouses. The major findings from this research will assist greenhouse growers, researchers, and engineers to choose the best radiation mode to model the thermal performance of a greenhouse depending on the greenhouse shape.

Description of the Experimental Site
The experimental greenhouse (Purme social farm, Figure 1) is located in Ohakdong, Yeoju-si, a small town in Gyeonggi-do, western South Korea. The experimental site was in a temperate climate with lengthy, hot, rainy, and partially cloudy summers, while winters are short, cold, snowy, and largely clear [30]. The farm covers a gross area of 4374 m 2 in three parts, farm A (2160 m 2 ), farm B (1782 m 2 ), and a packaging room (432 m 2 ), represented by A, B, and C, respectively, in    The greenhouse was fitted with passive, eco-friendly materials (transparent and opaque) for internal climate management. The sides and roof of the greenhouse were covered with horticultural glass (HG) and fluorine film. Retractable horizontal (Tempa and Luxous) and vertical (Obscura) thermal screens were installed ( Figure 4). The white surface of Obscura helps to reflect light onto the crop to enhance crop growth. Tempa is a shading screen used during the summer to minimise excessive solar radiation during A B C the day. It reflects and diffuses sunlight to moderate the solar heat load and cools the greenhouse. However, because most shading screens are used as thermal screens during winter, Tempa is also deployed at night to enhance heat retention. Luxous is a thermal screen that aids retaining heat at night to maintain a stable greenhouse temperature. The horizontal thermal screens divided the greenhouse into conditioned (zone 1) and unconditioned (zones 2 and 3) zones to reduce the greenhouse heating area. Sensors were fixed in each zone to monitor the effect of the screens on the macro-environment.  Outside weather data were recorded during winter between December 2021 and March 2022. The recorded data included the solar radiation, temperature, and relative humidity. Other required outside weather data for the BES model, such as wind speed, wind direction, and ambient pressure, were downloaded from the Korean Meteorological Administration (KMA), Icheon (Station 203, 37.26° N 127.48° E) [31]. The details of the recorded data are listed in Table 1. The windspeed data were measured at a height that was different from the eave height of the greenhouse. The data were modified according to the power law in Equation (1): where ℎ is the calculated wind speed (ms −1 ) at height h (m) of the greenhouse, is the downloaded KMA station wind speed (ms −1 ) at height ℎ (m), and ∝ is the power law exponent. The power law exponent is an empirically derived coefficient that increases in value with increase in terrain roughness. In an experiment conducted by Jung et al. [32] Obscura Tempa Luxous in a coastal region (Buan-gu) similar to the experimental site, the power law exponent was 0.28. The greenhouse design combines different strategies to supply heating energy via a fan coil unit, growing tubes, and tube rails ( Figure 5) to ensure adequate thermal comfort of crops under low ambient temperature conditions. The water temperature and flow rate were monitored to calculate the energy consumed using Equation (2). The locations and specifications of the water temperature and flow rate sensors are detailed in Figure 3 and Table 1, respectively.
where Q is the amount of energy consumed (kcalh −1 ), ̇ is the mass flow rate (kgh −1 ), is the specific heat capacity of water (kcalkg −1 −1 ), and ∆ is the difference between inlet and outlet water temperatures ( ).

Greenhouse Material Properties
Some properties of the greenhouse-covering materials have been presented in the literature. Valera et al. [33], Rafiq et al. [34], and Rabiu et al. [26] detailed the properties of HG, Tempa, and Luxous, respectively ( Table 2). Ground and steel are treated as opaque in TRNBuild (a TRNSYS building interface), with properties listed in Table 3. Fluorine film and Obscura are novel greenhouse materials whose properties were determined experimentally. The thickness (mm) of the fluorine film and Obscura was measured using an electronic digital calliper (TED PELLA, Inc., Redding, CA, USA), and the thermal conductivity (Wm −1 K −1 ) was measured using a QTM-500 (Kyoto Electronics MFG. Co., Ltd., Kyoto, Japan) thermal conductivity meter. The infiltration through the materials was measured using a procedure detailed in our previous study [26].
There are two main approaches for determining radiative properties: spectrophotometry and the radiation balance method (RBM) [35]. The suitability of fluorine film and Obscura in the BES model was confirmed using these two methods. The shortwave and longwave properties were determined in the laboratory using a haze meter (Murakami, HM-150) and Fourier transform infrared spectrophotometer. The RBM experimental procedure for longwave and shortwave radiation is shown in Figures 6 and 7.   The longwave radiation is calculated using Equations (3)-(5): where Qb is the upward longwave radiation from the material to the sky (Wm −2 ), Qc is the upward longwave radiation from the material to the black fabric during the night (Wm −2 ), Qd is the inward longwave radiation from the black fabric toward the material (Wm −2 ), Qa is the inward sky radiation toward the material (Wm −2 ), ρL is the longwave reflectance of the material, τL is the longwave transmittance of the material, and Es is the emissive power of the material (Wm −2 ) from the Stefan-Boltzmann's law. Es is calculated using Equation (6): Thermal screen Black fabric where is the Stefan-Boltzmann constant, and are the surface temperature and emissivity of the black fabric, respectively.
The shortwave radiation is calculated using Equations (7) and (8): where Sb is the outward shortwave radiation from the material toward the sky, Sc is the outward shortwave radiation from the material toward the black material, ρs is the reflectance of the material, τS is the transmittance of the material, and Sa and Sd are the downward shortwave radiation from the sky and the radiation from the black fabric toward the material, respectively.

Greenhouse Simulation Modelling in TRNSYS 18
Dynamic simulation of the greenhouse environment in TRNSYS was conducted by linking several components from the TRNSYS library. After connecting the required components, the time steps were specified to run the simulation. Figure 9 shows the simulation studio interface, and a dynamic flowchart of the modelling and simulation procedure is shown in Figure 10. Prior to the simulation, a three-dimensional (3D) greenhouse model was designed using the add-in TRNSYS3D for Google SketchUp ( Figure 11). In Google SketchUp, the greenhouse was divided into three zones separated by thermal screens, as shown in Figure 3.

Solar radiation processor
The screens use the same shading device principle as that in a residential building to limit the conductive and radiative heat losses from the greenhouse. Similar to a residential building where the shading device is designed to open or close based on solar radiation, this study controls the thermal screens based on solar radiation using Equations (9) and (10). Equations (9) and (10) are for Tempa and Luxous control, respectively: The Boolean controller controls the closing (bool = 1) and opening (bool = 0) of the screens; SR is the solar radiation, and Ti is the greenhouse's internal temperature.
The 3D model was generated in the Intermediate Data Format using Google SketchUp and imported into TRNBuild. TRNBuild was used to define the building orientation and material properties. The TRNBuild database did not include the properties of the greenhouse covering materials and thermal screens used in this study. Therefore, the material properties detailed in Section 2.2 were defined using Berkeley Lab WIN-DOW 7.4 program, generating a DOE-2 file, and the file was imported into TRNBuild.
The building with airflow (Type 56) was selected in the TRNSYS 18 simulation studio interface to assess the internal greenhouse thermal comfort based on outside weather conditions. The recorded outside weather data were then input into the weather data component (Type 9). The solar radiation recorded was the total horizontal radiation on the horizontal surface. To calculate the direct, diffuse, and reflected solar radiation, Type 16c was linked with Type 9 to accurately account for the solar heat gain of the greenhouse. Type 16c uses different modes to calculate the tilted surface radiation. Rasheed et al. [36] calibrated a TRNSYS model by varying the Type 16c mode. They concluded that mode 2 was the most accurate for modelling the greenhouse internal temperature when the total horizontal radiation was input in Type 9.
Occupancy and equipment were implemented as internal heat gains in TRNBuild and were provided with a schedule using the Type 14 forcing function in the simulation studio. The ASHRAE standard 130 W gain type was chosen for the occupancy, and a radiative and convective power of 0.5 kJh −1 was used for the equipment gain [37]. However, internal gain due to plant evapotranspiration was not considered.
The sensible heat balance of the air within the thermal zones is described as follows: ̇, =̇, +̇, +̇, +̇, , +̇, +̇, +̇, , where ̇, is the sensible heat flux of the zone (kJh −1 ), ̇, is the convective gain from surfaces (kJh −1 ), ̇, is the infiltration gains (kJh −1 ), ̇, is the ventilation gains (kJh −1 ),̇, , is the internal convective gains (kJh −1 ), ̇, is the gains due to interconnected air nodes (kJh −1 ), ̇, is the solar radiation entering an air node through external windows, transformed immediately into a convective gain to the internal air (kJh −1 ), and ̇, is the solar radiation absorbed on all internal shading devices of the zone transformed immediately into a convective gain to the internal air (kJh −1 ). The latent heat flux of the air node is described as follows [37]: where ̇, is the latent energy flux of the zone (kJh −1 ), ℎ is the heat of vaporisation of water (kJkg −1 ), ̇ is the mass flow rate of infiltration air (kgm −3 ), is the ambient humidity ratio ( kg kg −1 ), is the air node humidity ratio ( kg kg −1 ), ̇ is the mass flow rate of ventilation air (kgm −3 ), is the humidity ratio of ventilation air (kg kg −1 ), is the internal humidity gain (kg h −1 ), ̇ is the mass flow rate due to couplings of two zones (kgm −3 ), is the adjacent air node humidity ra-tio (kg kg −1 ), is the effective moisture capacitance (kg), and ∆ is the change in time step.
The airflow due to natural ventilation was estimated using TRNFLOW. In TRN-FLOW, fans, straight ducts, cracks, and large opening are the different ways of modelling the airflow network. Since the ventilation in the greenhouse involves the opening and closing of the screen, the large opening window was selected. The equation that governs the large opening airflow network is given as follows [38]: where ̇ is the mass flow rate of air (kgh −1 ), 1 → 2 is the flow direction from one air node to another, is the discharge coefficient, H is the total height of the window (m), ( ) is the air density (kgm −3 ) at height z, z is the height of the opening (m), is the angle of opening (°), W is the width of rectangular opening (m), and f(z) is the pressure difference at height z (pa).
The thermal zone energy balance is described as follows [39]: where is the temperature of the masses (°C), is the thermal capacitance of the zone masses (Kj°C −1 ), ̇, is the convective heat flux between the zone and the inner surface (kJh −1 ) , ̇, is the longwave radiation exchange between two inner surfaces (kJ h −1 ) , ̇, and ̇, are the radiative solar and internal gains, respectively (kJh −1 ),̇, is the convective heat flux between the external surface and ambient (kJh −1 ), ̇ is the absorbed solar gain on the outside opaque surfaces (kJh −1 ), ̇, is the longwave radiation emitted by the outside surfaces to their surroundings (kJh −1 ), and ̇ is the heat conduction through the building envelope (kJh −1 ). Different radiation modes are available for computing the direct and diffuse shortwave and longwave radiations within a thermal zone in the TRNbuild. The radiation modes are explained in the following subsections.

Beam Radiation Distribution
Standard and detailed modes are available for the beam radiation distribution. The standard mode uses a user-defined factor called GEOSURF to distribute beam radiation within a zone. GEOSURF is the fraction of beam radiation that strikes the individual surfaces with a maximum value of 1 per zone. The detailed mode for the beam radiation distribution was used to distribute the solar radiation entering the zone through external windows. An external program, TRNSHD, is used by TRNBuild to generate the shading and insolation matrices required for the detailed mode. However, for beam radiation entering the zone through an adjacent window, GEOSURF was also applied.

Diffuse Radiation Distribution
Standard and detailed modes based on different principles are available for the diffuse radiation distribution. The standard mode is based on the absorption-transmission weighted area ratios for all surfaces, whereas the detailed mode uses the Gebhart factor for shortwave diffuse radiation. An auxiliary program, TRNVFM, based on the view factor matrix is used in the detailed mode.

Longwave Radiation Distribution
Simple, standard, and detailed modes are available for the longwave radiation within a zone. The simple mode is a one-node model with a combined transfer function, assuming a constant convective coefficients and longwave radiation resistance for the heat transfer process.
The standard mode calculates simultaneously the longwave radiation exchange between the surfaces in an air node and the convective heat flux from the inside surfaces to the air node using a star network. This approach uses an artificial temperature of the air node (Tstar) to consider simultaneously the energy flow from a wall surface by convection to the air node and radiation to other surfaces. A mathematical description of the internal surface is given in Equations (15)-(17) [40]: , , = , , + , , , , , where , and , are the resistance of each surface and the equivalent resistance of all the surfaces (Ohms), respectively, , , is the convective heat flux (kJh −1 ), , , is the radiative heat flux (kJh −1 ), , , is the combined convective and radiative heat flux (kJh −1 ), , is the inside surface area (m 2 ), and , and are the surface and the equivalent air node temperatures (°C), respectively.
In addition, the total heat transfer at the outside surface is the sum of the convective and radiative heat fluxes calculated from Equations (18) to (21) [37]. The convective heat transfer at the outside surface is complicated due to wind direction, building size, local wind velocity, building surrounding, and building size [40].
, , = , , , , , , where , , is the convective heat flux to the outside surface (kJh −1 ), , , is the radiative heat flux to the outside surface (kJh −1 ), , , is the combined convective and radiative heat flux to the outside surface (kJh −1 ), ℎ , , is the convective heat transfer coefficient at the outside surface, , is the outside surface temperature (°C), , is the ambient temperature (°C), , is the longwave emissivity of the outside surface from the WINDOW library, is the Stefan-Boltzmann constant, is the fictive temperature difference between the ground and sky (°C), , is the view factor of the sky, is the fictive ground temperature (°C), 1 − , is the view factor of each external surface, and is the fictive sky temperature (°C). Compared with the standard mode, the artificial star node is not considered in the detailed mode because the convective heat flux is calculated separately. The detailed mode used the Gebhart factor ( ,1−2 ) for longwave radiation exchange, considering multi-reflection from surfaces. The Gebhart factor is the fraction of emission from surface 1 that reaches surface 2 and is absorbed. ,1−2 includes all paths (direct and multiple reflections) to reach 2 . The Gebhart matrix for longwave radiation is given by Equation (22): where and are diagonal matrices describing reflectivity and emissivity, I is an identity matrix, and F is the view factor, which is the fraction of diffusely radiated energy leaving surface 1 that is incident on surface 2; 'ir' represents the longwave range of the radiation spectrum (infrared).
The auxiliary matrix ( * ), given by Equation (23), was introduced to determine the longwave radiation (̇) in the enclosure: where is the transpose of , A is the diagonal matrix describing the surface areas, is the Stefan-Boltzmann constant, and is the longwave emissivity. Then, where T is the temperature vector of the enclosure. Considering that thermal screens and roofs are opened during the day to harness solar radiation, the effects of the beam and diffuse radiation modes are assumed to be negligible in greenhouses compared with residential buildings. However, during the night, the thermal screens and roofs are fully closed when energy is supplied to the conditioned zone of the greenhouse, and the properties of the thermal screens and covering materials determine the conditions of other zones. Therefore, different longwave radiation exchange modes impact the thermal comfort of the greenhouse.
The simple and standard longwave radiation modes cannot consider radiation exchange over more than one air node. The greenhouse was modelled as a single air node per zone as shown in Figure 11, to investigate the effects of the three longwave radiation modes. The shape of a building must be convex and closed to generate its view factor matrix for the detailed longwave radiation mode. Farm A is a concave-shaped greenhouse, and farm B is a convex-shaped greenhouse ( Figure 12). Consequently, for farm A, the simple and standard longwave radiation models were utilised, whereas for farm B, the simple, standard, and detailed modes were adopted.

BES Model Validation
The proposed BES model was validated using the Nash-Sutcliffe efficiency coefficient (NSE) represented by Equation (25). This coefficient quantifies the fit of the experimentally measured data with the simulated data in a 1:1 plot. Its value ranges from −∞ to 1, with values closer to 1 indicating substantial predictive potential of the model: where is the experimentally measured data, is the simulated data, and is the mean of the experimentally measured data.

Sensitivity Analysis
A sensitivity analysis was performed to determine the impact of the longwave radiation modes on the total monthly heating demand. In mathematical modelling, sensitivity analysis quantifies the effects of several independent variables on the dependent variable under a specific condition, and the sensitivity coefficient (SC) is a dimensionless factor widely used to characterise the error assessment. The simplest form of the sensitivity X Y X Y a b coefficients for comparative energy studies, as given by Lam and Hui [41], is presented in Equation (26): where OP is the output, and IP is the input. In addition to the radiation modes, several other factors, including infiltration, moisture change, and capacitance, affect the heating requirement of a greenhouse in the BES model. Rather than changing these parameters, the output of the heating requirements from the radiation modes is considered as the OP, whereas the measured heating consumption from the greenhouse was considered as the IP.

Radiometric Properties of the Novel Greenhouse Materials
The experiment based on RBM was conducted for several nights, and the results were compared with those by spectrophotometry. Figures 13 and 14 show the longwave radiative properties of the fluorine film and Obscura measured from 18:00 to 06:00 based on RBM. The reflectance, transmittance, and emittance of the fluorine film were 0.43, 0.94, and 0.02, respectively. The reflectance, transmittance, and emittance of Obscura were 0.96, 0.001, and 0.045, respectively. Because the Obscura has approximately zero transmittance in the longwave spectrum, spectrophotometry was used only for the fluorine film within the waveband of 2.5 to 25 μm (Figure 15a). The longwave transmittance result by spectrophotometry was the same as that by RBM. The consistency of the results using the two methods validates the accuracy in the determination of the radiative properties of the materials used in the BES modelling (Table 4).    -- The diffuse shortwave transmittance of the fluorine film, measured using an ultraviolet-visible spectrophotometer (Jasco V-730), is shown in Figure 15b. The measured diffuse transmittance is 45.86%. However, owing to the limitations of spectrophotometry, the total transmittance was measured using a haze meter. The total transmittance was 92.49% with the beam and diffuse transmissions contributing to 46.74% and 45.75%, respectively. The solar transmittance of the fluorine film indicates that it transmits more than 90% of the shortwave radiation, which is vital for plant growth. This value is higher than the transmittance of the commonly used greenhouse covering materials, polyethylene (71%) [34] and HG (89%) [33]. The solar transmittance obtained is similar to the longwave transmittance; this does not agree with Choab et al. [42], who recommended high and low transmittances in the shortwave and longwave spectra, respectively, for covering materials. However, the installed thermal screens have very low transmittance in the longwave spectrum (Tempa, 5%; Luxous, 33%) to limit heat loss before transmission through the film.

Comparison of Results of the TRNSYS Model with the Experimental Measurements
The daily average outside temperature and average solar radiation on the horizontal surface are shown in Figure 16  The simulated internal temperatures in the different zones for both farms were compared with the experimentally measured temperatures of the greenhouse. Because the daily greenhouse conditioning is approximately the same during the winter season, only a few days of the analysis are discussed. Figure 17 shows the simulated zone 1 temperatures compared with the experimentally measured temperatures from 26 to 31 January. On most days, the greenhouse's experimental daytime temperature was approximately 30 °C for both the farms. For farm A, the simulated daytime temperatures differ from the observed temperatures by a maximum of 3.0 °C (above the observed temperature) and 3.0 °C (below the observed temperature) for the simple and standard modes, respectively. For farm B, the simulated maximum daytime temperatures differ from the observed temperatures by a maximum of 3.0 °C (above the observed temperature), 2.3 °C (below the observed temperature), and 0.94 °C (below the observed temperature) for the simple, standard, and detailed modes, respectively. The simple and standard modes over-predicted and under-predicted the daytime temperatures, respectively, for both greenhouses. The controlled simulated night temperature for the three modes was 15 °C. So, it was not possible to determine the impact of the modes on the night temperature in zone 1. The calculated NSE values for farm A were 0.85 and 0.78 for the simple and standard modes, respectively. The calculated NSE values for farm B were 0.73, 0.78, and 0.83 for the simple, standard, and detailed modes, respectively.  Figure 18 shows the simulated zone 2 temperatures compared with the experimentally measured temperatures. The thermal screens placed between zones 1 and 2 caused an approximately 5 °C temperature difference between the two zones during the night when zone 1 was being heated. For farm A, the simple and standard modes overpredicted the daytime temperatures by a maximum of 8 °C and 3 °C, respectively. Similarly, for farm B, the simulated daytime temperatures differ from the observed temperatures by a maximum of 3.4 °C, 1.72 °C, and 0.8 °C above the observed temperature for the simple, standard, and detailed modes, respectively. However, the simple and standard modes under-predicted the night temperature by a maximum of 5.34 °C and 2.54 °C, respectively, in Farm A. By contrast, for farm B, the maximum difference of the simple, standard, and detailed modes was 3.5 °C, 3.7 °C, and 3.4 °C below the observed night temperatures, respectively. The calculated NSE values for farm A were 0.79 and 0.89 for the simple and standard modes, respectively. The calculated NSE values for farm B were 0.89, 0.88, and 0.9 for the simple, standard, and detailed modes, respectively.  Figure 19 shows the simulated zone 3 temperatures compared with the experimentally measured temperatures. The thermal screens placed between zones 2 and 3 caused an 8℃ temperature difference between the two zones during the night. This temperature difference is higher than between Zone 1 and 2 owing to the low transmissive properties of Tempa compared with Luxous. The average temperature difference between zone 3 and the ambient was 10 °C. The calculated NSE values for farm A were 0.91 and 0.93 for the simple and standard modes, respectively. The calculated NSE values for farm B were 0.92, 0.91, and 0.92 for the simple, standard, and detailed modes, respectively.
As the NSE values for all radiation modes were higher than 0.5, the BES model was valid, and the standard and detailed models were recommended for predicting the internal greenhouse temperature for concave and convex-shaped greenhouses. However, the user-defined emissivities of the surface for longwave radiation can only be handled Because the overall performance of the standard and detailed radiation modes for farms A and B was superior, both were used to analyse the heating demand of the two greenhouses. However, the maximum heating demand varied for the different radiation modes (Table 5). Figure 20 compares the simulated heating demand with the experimental energy consumption during the period of the lowest ambient temperature from 21 to 27 December 2021. There were significant hourly changes between the simulated and measured results owing to the transient state of the experimental greenhouse and the intermittent application of the energy control systems. The maximum heating loads for farms A and B were 113.5 kcal/hm 2 and 120.7 kcal/hm 2 , respectively, on 26 December when the ambient temperature (−17.8 ℃) was minimum and the greenhouse set point temperature was 15 ℃. For the same greenhouse with similar hybrid heat pump systems using geothermal sources and solar heat, Jeon et al. [43] designed a greenhouse (1015 m 2 ) with a maximum heating load of 148.8 kcal/hm 2 for the lowest ambient temperature of −19 ℃ and greenhouse set point temperature of 23 ℃. The simulated maximum heating demands for both the farms with the thermal screens were high compared with Rasheed et al. [23], who obtained 109.9 kcal/hm 2 , 98 kcal/hm 2 , and 81.3 kcal/hm 2 using single-, double-, and triple layered-thermal screens, respectively, for a greenhouse floor area of 7572.6 m 2 . However, the simulated heating demand was lower than that of a greenhouse with a floor area of 391.2 m 2 and with a maximum heating load of 250 kcal/hm 2 [44]. This trend indicates that greenhouse heating demand per square meter decreases with increasing floor area. Heat is lost through all greenhouse surfaces, including walls, roofs, and floors, and the amount of heat lost per square meter increases as the wall-to-total surface area ratio increases. The total daily heating demand was compared with the total daily energy consumed for BES model performance. The maximum heating demand was on 19 January (farm A: 3547.8 Mcal/day, and farm B: 3118.9 Mcal/day) when the daily average temperature was −7.3℃. In the study by Ahamed et al. [22] using the CSGHEAT models with MATLAB, the average daily heating requirement in a CSG integrated with thermal screens was 900 MJ/day (215.2 Mcal/day). The CSG heating area was 210 m 2 . This is approximately ten times lower than that of the experimental glass greenhouse. The calculated NSE values for farms A and B were 0.89 and 0.9, respectively, indicating the model's ability to predict the greenhouse energy load.  Figure 21 shows the monthly heating demands predicted by the radiation modes and the experimental heating energy consumption. The radiation modes affected the simulated heating requirement of the greenhouse. The monthly heating demand predicted by the simple and standard radiation modes for farm A matched the experimental measurements with deviations within 27.7% and 7.6%, respectively. The monthly heating loads predicted by the simple, standard, and detailed radiation modes for farm B were similar to the experimental measurements with deviations within 10.5%, 6.7%, and 2.9%, respectively. For both greenhouses, the monthly heating requirements predicted by the simple mode were less than the experimental measurements. However, the standard mode under-and over-predicted the monthly heating requirements in farms A and B, respectively.  Figure 22 shows the SC for both the farms. The simple and standard modes for farm A had an average SC of 0.7 and 0.9, respectively, whereas the simple, standard, and detailed modes for farm B had an average SC of 0.9, 1.1, and 1.0, respectively. In addition, for farm B, as the monthly ambient temperature increases, the heating demand predicted by the standard and detailed modes converged, indicating that the differences in the predictions between the radiation modes were more pronounced during cold months for a convex-shaped greenhouse.

Limitation of the Study and Future Work
The proposed BES model has some limitations. Firstly, three weather data (windspeed, wind direction, and ambient pressure) used in the simulation were downloaded at a neighbouring station, 18 km from the experimental greenhouse. The results obtained could lead to a slight variation if the model is to be adapted for year-round simulation, as the greenhouse solely depends on wind forces for natural ventilation in the summer. However, the division of the greenhouse based on the geometry shows the variation in the radiation mode on the internal temperature and heating demand of the greenhouse. More so, in the absence of on-site data, several researchers have used neighbouring stations weather data. Choi et al. [45] had studied the windspeed trend of 31 KMA stations from 1985-2019 and concluded that the average wind speed in Korea does not change that much, while Kim et al. [46] developed a TRNSYS model for Jincheon using Cheongju weather data, 23.9 km away from Jincheon. Secondly, the planted crops were not considered in the simulation. Although they are the most important factor in greenhouse simulation [47], the absence of a stand-alone component in TRNSYS limited their consideration in our study. The present TRNSYS model could only access the greenhouse thermal performance based on radiation modes without including the internal heat gains due to evapotranspiration from the crop. Modified crop models to calculate the latent mass flow and latent power due to the presence of crop will be implemented in future.

Conclusions
The greenhouse industry is faced with challenges to increase production and reduces the use of resources. Mathematical modelling has been seen as alternative to identify potential solutions. In this study, TRNSYS was used to investigate the thermal performance of a glass greenhouse integrated with movable thermal screens. For the novel greenhouse covering materials and thermal screens considered in the study, the thermophysical, aerodynamics, and radiative properties were experimentally determined for accurate BES modelling. The internal greenhouse temperature and heating demand were simulated using different longwave radiation modes, and the simulated results were validated using transient experimental data. The performances of the standard and detailed modes were superior in concave-and convex-shaped greenhouses as the radia- tion modes were sensitive to ambient temperature. The findings of this study will aid greenhouse growers, researchers, and engineers in selecting the appropriate longwave radiation mode for modelling the thermal performance of a greenhouse based on the geometry of the greenhouse. Data Availability Statement: Data available on request due to restrictions.

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

Nomenclature Symbols
, Inside surface area, m 2 Thermal capacitance of the zone masses, kJ −1 Specific heat capacity of water, kcalkg −