− Anisotropic Models

Anisotropic models are applied with a partially cloudy sky, which can vary between clear and cloudy sky. Below are three models for this case.

### (a) Temps and Coulson Model

The Temps and Coulson model of anisotropic distribution for clear skies was developed in [14] and can be calculated with Equation (8).

$$\mathbf{I}\_{\mathbf{s}} = \frac{1}{2} \mathbf{I}\_{\mathbf{d}} (1 + \cos \theta) \left[ 1 + \operatorname{sen}^3 \left( \frac{\mathfrak{d}}{2} \right) \right] \left( 1 + \cos^2 \theta \operatorname{sen}^3 \theta\_{\mathbf{z}} \right) \tag{8}$$

This model includes correction factors for the isotropic diffuse irradiation model which take into account the zones of anisotropy in diffuse irradiation. Factor [1 + sen<sup>3</sup> (β/2)] is included to explain the increase in skylight observed near the horizon on clear days, the factor (1 + cos2θ sen3θz) approximates the brightness of the sky near the Sun, and the third factor represents the reflection of the earth in a better way.

### (b) Klucher Model

The Klucher model of anisotropic distribution for all of sky types [37] modifies the formulation of the Temps and Coulson model by including a factor F = 1 <sup>−</sup> (Id/I)<sup>2</sup> , as indicated in Equation (9).

$$\mathbf{I}\_{\mathbf{s}} = \frac{1}{2} \mathbf{I}\_{\mathbf{d}} (1 + \cos \theta) \left[ 1 + \mathbf{F} \sin^{3} \left( \frac{\mathfrak{d}}{2} \right) \right] \left( 1 + \mathbf{F} \cos^{2} \theta \,\mathrm{sen}^{3} \theta\_{\mathbf{z}} \right) \tag{9}$$

When the sky is completely cloudy, F = 0 (i.e., the equation returns to the isotropic model), and when the sky is clear, F = 1 (i.e., the Temps and Coulson model is used), thus improving the estimates for all types of sky.

### (c) Hay Model

The Hay model utilizes a circumsolar component that comes directly from the direction of the Sun and another component of diffuse irradiation that is distributed isotropically from the rest of the celestial vault, as calculated by Hay [38,39] with Equation (10).

$$\mathbf{I}\_{\mathbf{s}} = \mathbf{I}\_{\mathbf{d}} \left\{ \frac{\mathbf{I} - \mathbf{I}\_{\mathbf{d}}}{\mathbf{I}\_{0}} \mathbf{r}\_{\mathbf{b}} + \frac{1}{2} (\mathbf{1} + \cos \beta) \left[ \mathbf{1} - \frac{\mathbf{I} - \mathbf{I}\_{\mathbf{d}}}{\mathbf{I}\_{0}} \right] \right\} \tag{10}$$

These two components are weighted according to an isotropy index (i.e., a ratio of the direct horizontal solar irradiation on the earth and the extraterrestrial solar irradiation).

#### 2.2.4. Global Solar Irradiation Incident on an Inclined Plane The total amount of the solar radiation incident on an inclined plane is made up of

The total amount of the solar radiation incident on an inclined plane is made up of three components (direct, reflected from the ground, and diffused from the sky) which are then combined once the individual values are known. In places where the hourly solar irradiation (global and diffuse) on horizontal surfaces is known or the latter can be estimated, the global irradiation on an inclined plane can be written, as in Iqbal [35], with Equations (11) or (12). three components (direct, reflected from the ground, and diffused from the sky) which are then combined once the individual values are known. In places where the hourly solar irradiation (global and diffuse) on horizontal surfaces is known or the latter can be estimated, the global irradiation on an inclined plane can be written, as in Iqbal [35], with Equations (11) or (12). Iஒஓ = ሺI−Iୢሻrୠ + I୰ + Iୱ (11)

$$\mathbf{I}\_{\mathbf{\hat{B}}\,\mathbf{\hat{Y}}} = (\mathbf{I} - \mathbf{I}\_{\mathbf{d}})\mathbf{r}\_{\mathbf{b}} + \mathbf{I}\_{\mathbf{r}} + \mathbf{I}\_{\mathbf{s}} \tag{11}$$

$$\mathbf{I}\_{\mathsf{\B}\,\mathbf{\varDelta}} = \mathbf{I}\_{\mathsf{\B}\,\mathbf{\varDelta}\,\mathbf{\varDelta}} + \mathbf{I}\_{\mathsf{\mathsf{d}}\,\mathbf{\varDelta}} \tag{12}$$

where: where:

I: global hourly irradiation incident on a horizontal surface; Iβγ: global hourly irradiation incident on an inclined and oriented plane; Ibβγ: direct hourly irradiation incident on an inclined and oriented plane; Idβ: diffuse hourly irradiation incident on an inclined and oriented plane (I<sup>r</sup> + Is). I: global hourly irradiation incident on a horizontal surface; Iβγ: global hourly irradiation incident on an inclined and oriented plane; Ibβγ: direct hourly irradiation incident on an inclined and oriented plane; Idβ: diffuse hourly irradiation incident on an inclined and oriented plane (Ir + Is).

### *2.3. Simulink-MATLAB Methodology for Estimating the Solar Irradiation Incident on the Inclined Plane 2.3. Simulink-MATLAB Methodology for Estimating the Solar Irradiation Incident on the Inclined Plane*

For the calculation of solar irradiation on an inclined and/or oriented surface, existing models described in the literature need the horizontal global solar irradiation, horizontal diffuse solar irradiation, solar height, and angle of incidence as input variables, with their values either measured or estimated. When other models are used for the calculation (e.g., horizontal diffuse solar irradiation), the value of extraterrestrial solar irradiation is provided through parameters such as the solar constant, a correction factor for Earth's eccentricity, the day of the year, the solar declination, and the solar time angle, as well as others pertaining to the location of the greenhouse, such as its geographical latitude and longitude, and the inclination and orientation of the solar panels. For the calculation of solar irradiation on an inclined and/or oriented surface, existing models described in the literature need the horizontal global solar irradiation, horizontal diffuse solar irradiation, solar height, and angle of incidence as input variables, with their values either measured or estimated. When other models are used for the calculation (e.g., horizontal diffuse solar irradiation), the value of extraterrestrial solar irradiation is provided through parameters such as the solar constant, a correction factor for Earth's eccentricity, the day of the year, the solar declination, and the solar time angle, as well as others pertaining to the location of the greenhouse, such as its geographical latitude and longitude, and the inclination and orientation of the solar panels.

In this study, the methodology was developed with models from different authors that have been accepted in the literature; however, to facilitate its use, the Simulink-MATLAB platform for programming with visual objects was deployed. The modular nature of its design provides the possibility of using different existing models for different case studies, depending on which is the most suitable, or creating new models to obtain the various variables required. In this study, the methodology was developed with models from different authors that have been accepted in the literature; however, to facilitate its use, the Simulink-MATLAB platform for programming with visual objects was deployed. The modular nature of its design provides the possibility of using different existing models for different case studies, depending on which is the most suitable, or creating new models to obtain the various variables required.

In this study, the value of the hourly global solar irradiation on an inclined plane, according to the scheme proposed in Figure 1, was calculated with six functional blocks: In this study, the value of the hourly global solar irradiation on an inclined plane, according to the scheme proposed in Figure 1, was calculated with six functional blocks:


This methodology received the following as input variables: This methodology received the following as input variables:


The rest of the parameters were set inside each block and identified the position of the solar panels of the photovoltaic system and the location of the greenhouse: The rest of the parameters were set inside each block and identified the position of the solar panels of the photovoltaic system and the location of the greenhouse:


*Agronomy* **2021**, *11*, x FOR PEER REVIEW 8 of 23

**Figure 1.** Methodology for estimating the measured values of global irradiance on a horizontal surface with the estimates received by the inclined plane of the solar panel. **Figure 1.** Methodology for estimating the measured values of global irradiance on a horizontal surface with the estimates received by the inclined plane of the solar panel.

The following sections detail the implementation of each block of the methodology developed in Simulink-MATLAB. The following sections detail the implementation of each block of the methodology developed in Simulink-MATLAB.

2.3.1. Hourly Extraterrestrial Solar Irradiation Block 2.3.1. Hourly Extraterrestrial Solar Irradiation Block

The Sun emits a flow of energy which is considered to be almost constant except for small variations due to sun spots. This solar energy, when it reaches the top of the atmosphere, receives the name of extraterrestrial solar irradiation and it is this solar irradiation that would reach a certain point on Earth if the atmosphere that protects it did not exist. Extraterrestrial solar irradiation for various latitudes can be estimated from the fol-The Sun emits a flow of energy which is considered to be almost constant except for small variations due to sun spots. This solar energy, when it reaches the top of the atmosphere, receives the name of extraterrestrial solar irradiation and it is this solar irradiation that would reach a certain point on Earth if the atmosphere that protects it did not exist.

lowing parameters: the solar constant, the solar declination, and the time of year. For hourly or shorter periods, the solar angle at the beginning and the end of the period has to be considered (see the work of Allen [40], recommended by the FAO for the development of calculations for crop evapotranspiration, Duffie et al. [41], and Kalogirou [42]). This is done with Equation (13). Extraterrestrial solar irradiation for various latitudes can be estimated from the following parameters: the solar constant, the solar declination, and the time of year. For hourly or shorter periods, the solar angle at the beginning and the end of the period has to be considered (see the work of Allen [40], recommended by the FAO for the development of calculations for crop evapotranspiration, Duffie et al. [41], and Kalogirou [42]). This is done with Equation (13).

$$\mathbf{I}\_{0} = \left(\frac{12 \cdot 60}{\pi}\right) \mathbf{I}\_{\rm sc} \to \left[ (\omega\_{2} - \omega\_{1}) \text{ sen } \phi \text{ sen } \delta + \cos \phi \text{ cos } \delta (\text{sen } \omega\_{2} - \text{sen } \omega\_{1}) \right] \tag{13}$$

I0: hourly extraterrestrial solar irradiation incident on a horizontal surface, MJ/(m2·h); where:

Isc: solar constant, 0.082 MJ/(m2·min); E0: correction factor for the eccentricity of the Earth (r0/r)2, with E0 = 1 + I0: hourly extraterrestrial solar irradiation incident on a horizontal surface, MJ/(m<sup>2</sup> ·h); Isc: solar constant, 0.082 MJ/(m<sup>2</sup> ·min);

(0.033·cos(2·π·J/365)); r0: average Sun–Earth distance, 1 ua; E0: correction factor for the eccentricity of the Earth (r0/r)<sup>2</sup> , with E<sup>0</sup> = 1 + (0.033·cos(2·π·J/365)); r0: average Sun–Earth distance, 1 ua;

r: current Sun–Earth distance, ua; r: current Sun–Earth distance, ua;

ua: astronomical unit, 1496 × 108 km; ua: astronomical unit, 1496 <sup>×</sup> <sup>10</sup><sup>8</sup> km;

J: day number of the year, 1 for January 1, …, 365 for December 31; J: day number of the year, 1 for January 1, . . . , 365 for December 31;

ϕ: geographic latitude (rad), north (+) and south (−): −90° ≤ ϕ ≤ 90° where (rad) = π/180·(°decimal places); δ: solar declination (rad), with δ = 0.409·sen((2·π·J/365) − 1.39). Defined as the angular poφ: geographic latitude (rad), north (+) and south (−): −90◦ ≤ φ ≤ 90◦ where (rad) = π/180·( ◦decimal places);

sition of the Sun at solar noon—that is, when the Sun is in the local meridian—in relation

to the plane of the Earth's equator, north (+) and south (−): −23.45° ≤ δ ≤ 23.45°;

δ: solar declination (rad), with δ = 0.409·sen((2·π·J/365) − 1.39). Defined as the angular position of the Sun at solar noon—that is, when the Sun is in the local meridian—in relation to the plane of the Earth's equator, north (+) and south (−): −23.45◦ ≤ δ ≤ 23.45◦ ;

ω1: solar hour angle at the beginning of the period (rad), with ω<sup>1</sup> = ω − ((π·t1)/24);

ω: solar hour angle at the moment when the midpoint of the considered period occurs (rad), with ω = (π/12)·[(t + 0.06667·(L<sup>z</sup> − Lm) + Sc) − 12];

t1: duration of the period considered (hours), 1 for hourly periods, 0.5 for periods of 30 min; t: standard time at the midpoint of the period considered (hours) (e.g., for a period between 2:00 p.m. and 3:00 p.m., t = 14.5);

Lz: geographic longitude of the center of the local time zone, degrees west of Greenwich: 75◦ East, 90◦ Central, 105◦ Rocky Mountain, 120◦ Pacific USA, 0◦ Greenwich, 330◦ Cairo, 255◦ Bangkok, 345◦ Spain (Iberian Peninsula);

Lm: geographic longitude of the measurement area, degrees west of Greenwich;

Sc: seasonal correction for solar time (hours), with S<sup>c</sup> = 0.1645·sen(2b) − 0.1255·cos(b) − 0.025·sen(b) and b = 2·π·(J − 81)/364;

ω2: solar hour angle at the end of the period (rad), with ω<sup>2</sup> = ω + ((π·t1)/24). If, in the morning, ω < −ωs, or, in the afternoon, ω > ωs, this indicates that the Sun is below the horizon such that I<sup>0</sup> = 0.

It is necessary to take into account the advance of the clock time in the official time (i.e., wintertime UTC+1 from the last Sunday in October and summertime UTC+2 from the last Sunday in March).

### 2.3.2. Horizontal Diffuse Solar Irradiation Block

The horizontal diffuse solar irradiation modeling from the hourly clearness index (k<sup>t</sup> = I/I0) and an hourly diffuse fraction (k<sup>d</sup> = Id/I) was carried out with a third-order model (that of Miguel et al. [43]) and was performed using data from several countries in the northern Mediterranean, resulting in Equation (14).

$$\mathbf{k}\_{\rm d} = \begin{cases} 0.995 - 0.081 \,\mathbf{k}\_{\rm t} & \mathbf{k}\_{\rm t} \le 0.21\\ 0.724 + 2.738 \,\mathbf{k}\_{\rm t} - 8.32 \,\mathbf{k}\_{\rm t}^2 + 4.967 \,\mathbf{k}\_{\rm t}^3 & \text{for} \quad 0.21 < \mathbf{k}\_{\rm t} \le 0.76\\ 0.18 & 0.76 < \mathbf{k}\_{\rm t} \end{cases} \tag{14}$$

### 2.3.3. Solar Height and Zenith Angle Block

Solar height, also called solar elevation, is the angular height of the Sun above the observer's celestial horizon, which is an angle between 0◦ and 90◦ . The zenith angle, also called the zenith distance, is the angle between the local zenith and the line that joins the observer and the Sun; the value of the angle is between 0◦ and 90◦ . The solar height is the complement of the zenith angle. For a given geographic position, in the absence of atmospheric refraction of the earth, the trigonometric relationship between the Sun and a horizontal surface provided by Iqbal [35] is defined by Equation (15).

$$\cos\theta\_{\mathbf{z}} = \text{sen } \delta \text{ sen } \phi + \cos\delta \text{ } \cos\phi \text{ } \cos\omega = \text{sen } \mathfrak{a} \tag{15}$$

where:

α: solar height—the angle of elevation of the Sun above the true horizon; θz: zenith angle—the angular position of the Sun in relation to the local vertical, θ<sup>z</sup> = 90◦ − α.

### 2.3.4. Angle of the Incidence of Solar Irradiation on the Solar Panel Block

The angle of incidence is the angle formed between the normal to the inclined plane and the Sun–Earth vector. There are two cases: the inclined plane may be oriented towards the equator or with an arbitrary orientation towards east or west.

− Modeling the angle of incidence for an inclined solar panel oriented towards the equator

The angle of incidence for an inclined surface oriented towards the equator can be described, as detailed by Liu and Jordan [12], with Equation (16).

$$\cos\theta\_0 = \sin\delta \,\mathrm{sen}\,(\phi - \beta) + \cos\delta \,\cos(\phi - \beta) \,\cos\omega \tag{16}$$

where:

θ0: angle of incidence for an inclined surface oriented towards the equator; β: inclination of the surface to the horizontal position.

− Modeling the angle of incidence for an arbitrarily inclined and oriented solar panel

The relationship of the angle of incidence for a surface inclined and oriented in any direction with the local meridian is trigonometric (see Berod and Bock [44], Kondratyev [45], and Coffari [46]) and can be described with Equations (17) or (18).

$$\begin{aligned} \cos \theta &= (\text{sen } \phi \cos \beta - \cos \phi \text{sen } \beta \cos \gamma) \text{sen } \delta \\ &+ (\cos \phi \cos \beta + \text{sen} \phi \text{sen } \beta \cos \gamma) \cos \delta \cos \omega \\ &+ \cos \delta \text{sen} \beta \text{sen } \gamma \text{ sen } \omega \end{aligned} \tag{17}$$

where:

θ: the angle of incidence for an arbitrarily inclined and oriented surface.

$$\cos\theta = \cos\beta\cos\theta\_{\mathbf{z}} + \text{sen }\beta\cos\theta\_{\mathbf{z}}\cos(\psi-\gamma) \tag{18}$$

where:

γ: azimuth angle of the surface, orientation. Defined as the deviation of the normal to the surface of the solar collector from the local meridian in the directions west (−), south (0), and east (+);

ψ: solar azimuth with cosψ = ((senα·senφ-senδ)/(cosα·cosφ)), with 0◦ ≤ ψ ≤ 90◦ for cosψ ≥ 0, and with 90◦ ≤ ψ ≤ 180◦ for cosψ ≤ 0. This is the angle at the local zenith between the plane of the observer's meridian and the plane of a great circle passing through the zenith and the Sun in the directions west (−), south (0), and east (+): −180◦ ≤ 0 ◦ ≤ +180◦ .

Measurement of the angle of incidence of the direct solar irradiation can be done by constructing a simple device. This consists of mounting a normal pointer to a flat surface, on which graduated concentric circles are marked. The length of the shadow cast by the pointer can be measured using concentric circumferences and can be used to determine the angle of incidence according to the international standard ISO 9806:2017 [47]. The device should be located on the plane and to one side of the solar panel.

2.3.5. Angle of the Incidence of the Solar Irradiation on the Solar Panel Block

This block was used to unify the three components of solar radiation that affect the solar panel in order to obtain the amount of the global solar radiation incident on the inclined surface.

### 2.3.6. Conversion from Hourly Global Irradiance to Hourly Mean Solar Irradiance Block

Finally, the conversion of the global hourly solar irradiance (i.e., energy value) obtained for the inclined and oriented plane to the corresponding hourly mean solar irradiance (i.e., power value) was undertaken. Through interpolation, these hourly values can be used to create a database with a time interval of one minute in order to carry out a more detailed simulation of the operation of the photovoltaic system under study (e.g., as input values for the simulation of a model of a greenhouse photovoltaic system).

### **3. Results**

This section presents the results of the selected models for the estimation of the global solar irradiation value on the solar panel inclined plane, based on the global solar irradiation on the horizontal surface (which is usually recorded in meteorological stations) and using the following Simulink-MATLAB blocks:


### *3.1. Result of the Hourly Extraterrestrial Solar Irradiation Block 3.1. Result of the Hourly Extraterrestrial Solar Irradiation Block*  The extraterrestrial solar irradiation modeling was applied to determine the temporal

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 11 of 23

The extraterrestrial solar irradiation modeling was applied to determine the temporal evolution of the extraterrestrial solar irradiation at the top of the atmosphere (which would be found to reach a certain point on Earth if the attenuation and scattering effects which are produced by the atmosphere when the sun's rays pass through it are not considered). evolution of the extraterrestrial solar irradiation at the top of the atmosphere (which would be found to reach a certain point on Earth if the attenuation and scattering effects which are produced by the atmosphere when the sun's rays pass through it are not considered).

The results of the methodology for the calculation of hourly extraterrestrial solar irradiation, obtained for the 15th day of each month of the year at latitude 42◦ N and longitude 5.6◦ W, are shown in (Figure 2). The results of the methodology for the calculation of hourly extraterrestrial solar irradiation, obtained for the 15th day of each month of the year at latitude 42° N and longitude 5.6° W, are shown in (Figure 2).

**Figure 2.** Hourly extraterrestrial solar irradiation calculated at latitude 42° N and longitude 5.6° W for the 15th day of each month. (**Top**) January to June; (**bottom**) July to December. **Figure 2.** Hourly extraterrestrial solar irradiation calculated at latitude 42◦ N and longitude 5.6◦ W for the 15th day of each month. (**Top**) January to June; (**bottom**) July to December.

### *3.2. Result of the Hourly Horizontal Diffuse Solar Irradiation Block*

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 12 of 23

The different correlations proposed for the hourly diffuse fraction in the introduction are shown in Figure 3. Generally, each model includes three relationships for each sky type, as detailed by Iqbal [35], that determine the daily clearness index in order to define sky conditions such that: The different correlations proposed for the hourly diffuse fraction in the introduction are shown in Figure 3. Generally, each model includes three relationships for each sky type, as detailed by Iqbal [35], that determine the daily clearness index in order to define sky conditions such that: *Agronomy* **2021**, *11*, x FOR PEER REVIEW 12 of 23 The different correlations proposed for the hourly diffuse fraction in the introduction are shown in Figure 3. Generally, each model includes three relationships for each sky type, as detailed by Iqbal [35], that determine the daily clearness index in order to define sky conditions such that:


**Figure 3.** Different correlations of the diffuse solar irradiation on the horizontal surface proposed in the literature. **Figure 3.** Different correlations of the diffuse solar irradiation on the horizontal surface proposed in the literature.

Figure 4 shows the results of the methodology used for the conversion of the hourly diffuse solar irradiation to the hourly global solar irradiation, using data recorded by the State Meteorological Agency (AEMet) at La Virgen del Camino station (León, Castile and León region, Spain), with regard to the relation between the hourly diffuse fraction and the hourly clearness index, with values recorded for the central eight hours of the day during 2011, along with the model selected (Miguel et al. [43]). Figure <sup>4</sup> shows the results of the methodology used for the conversion of the hourlydiffuse solar irradiation to the hourly global solar irradiation, using data recorded by the State Meteorological Agency (AEMet) at La Virgen del Camino station (León, Castile and León region, Spain), with regard to the relation between the hourly diffuse fraction and the hourly clearness index, with values recorded for the central eight hours of the day during 2011, along with the model selected (Miguel et al. [43]). diffuse solar irradiation to the hourly global solar irradiation, using data recorded by the State Meteorological Agency (AEMet) at La Virgen del Camino station (León, Castile and León region, Spain), with regard to the relation between the hourly diffuse fraction and 

**Figure 4.** AEMet data from La Virgen del Camino (León, Castile and León region, Spain) for the hourly diffuse fraction vs. the hourly clearness index for the eight central hours of the day during the year 2011, along with the horizontal diffuse solar irradiation estimation model (Miguel et al. [43]). **Figure 4.** AEMet data from La Virgen del Camino (León, Castile and León region, Spain) for the hourly diffuse fraction solar irradiation estimation model (Miguel et al. [43]). **Figure 4.** AEMet data from La Virgen del Camino (León, Castile and León region, Spain) for the hourly diffuse fraction vs. the hourly clearness index for the eight central hours of the day during the year 2011, along with the horizontal diffuse solar irradiation estimation model (Miguel et al. [43]).

### *3.3. Result of the Hourly Solar Height Block*

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 13 of 23

The results of the methodology used for the modeling of the hourly solar height, obtained for the 15th day of each month of the year at latitude 42◦ N and longitude 5.6◦ W, are shown in Figure 5. The results of the methodology used for the modeling of the hourly solar height, obtained for the 15th day of each month of the year at latitude 42° N and longitude 5.6° W, are shown in Figure 5.

**Figure 5.** Hourly solar height calculated at latitude 42° N and longitude 5.6° W for the 15th day of each month. (**Top**) January to June; (**bottom**) July to December. **Figure 5.** Hourly solar height calculated at latitude 42◦ N and longitude 5.6◦ W for the 15th day of each month. (**Top**) January to June; (**bottom**) July to December.

Figure 6 shows that the maximum values for the solar height were obtained at solar midday. Figure 6 shows that the maximum values for the solar height were obtained at solar midday.

### *3.4. Result of the Hourly Incidence Angle Block*

Figure 7 shows the results of the methodology used for the modeling of the hourly incidence angle obtained for the surface of a solar panel with an inclination of 45◦ and oriented towards the equator, for the 15th day of each month of the year at latitude 42◦ N and longitude 5.6◦ W.

**Figure 6.** Solar height at solar noon calculated at latitude 42° N and longitude 5.6° W for each day of the year. **Figure 6.** Solar height at solar noon calculated at latitude 42◦ N and longitude 5.6◦ W for each day of the year. incidence angle obtained for the surface of a solar panel with an inclination of 45° and oriented towards the equator, for the 15th day of each month of the year at latitude 42° N and longitude 5.6° W.

*3.4. Result of the Hourly Incidence Angle Block* 

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 14 of 23

**Figure 7.** Hourly incidence angle calculated at latitude 42° N and longitude 5.6° W for the 15th day of each month. (**Top**) January to June; (**bottom**) July to December. **Figure 7.** Hourly incidence angle calculated at latitude 42◦ N and longitude 5.6◦ W for the 15th day of each month. (**Top**) January to June; (**bottom**) July to December.

**Incidence Angle at Solar Noon: 42ºN 5.6ºW. Tilted 45º. Oriented 0º**

**Figure 8.** Incidence angle at solar noon calculated at latitude 42° N and longitude 5.6° W for all days of the year.

**01/01 01/02 01/03 01/04 01/05 01/06 01/07 01/08 01/09 01/10 01/11 01/12 31/12 Day of the Year**

> the following four models. − Inclined Model 1

**0**

**5**

**10**

**15**

**Incidence Angle, º**

**20**

**25**

**30**

*3.5. Results for Hourly Global Irradiance and Hourly Mean Solar Irradiance on the Solar Panel* 

of the different methods of obtaining diffuse solar irradiation are shown below.

The data recorded for the global solar irradiation on the horizontal surface during 2011 in the agrometeorological station located in the town of Mansilla Mayor (León, Castile and León region, Spain), as part of the SIAR network, were used to apply the methodology developed in Simulink-MATLAB for the estimation of the global solar irradiation value on a surface with an inclination of 45° and oriented towards the equator. The results

Figure 9 shows the values of the daily global solar irradiation obtained from the horizontal SIAR network, together with the calculations carried out for the estimation of the daily global solar irradiation on a surface inclined at 45° and oriented to the equator using January to June; (**bottom**) July to December.

**6h 7h 8h 9h 10h 11h 12h 13h 14h 15h 16h 17h 18h Solar Hour of the Day**

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 15 of 23

**<sup>100</sup> Hourly Incidence Angle: 42ºN 5.6ºW. Tilted 45º. Oriented 0º**

The minimum values of the angle of incidence, which were obtained at solar noon for each day of the year, are presented in Figure 8. The minimum values of the angle of incidence, which were obtained at solar noon for each day of the year, are presented in Figure 8.

**July August September October November December**

**Figure 8.** Incidence angle at solar noon calculated at latitude 42° N and longitude 5.6° W for all days of the year. **Figure 8.** Incidence angle at solar noon calculated at latitude 42◦ N and longitude 5.6◦ W for all days of the year. nent, therefore obtaining different results.

*3.5. Results for Hourly Global Irradiance and Hourly Mean Solar Irradiance on the Solar Panel 3.5. Results for Hourly Global Irradiance and Hourly Mean Solar Irradiance on the Solar Panel* Inclined Model 2 used the Liu and Jordan isotropic model, presented in Equation (7),

The data recorded for the global solar irradiation on the horizontal surface during 2011 in the agrometeorological station located in the town of Mansilla Mayor (León, Castile and León region, Spain), as part of the SIAR network, were used to apply the methodology developed in Simulink-MATLAB for the estimation of the global solar irradiation value on a surface with an inclination of 45° and oriented towards the equator. The results The data recorded for the global solar irradiation on the horizontal surface during 2011 in the agrometeorological station located in the town of Mansilla Mayor (León, Castile and León region, Spain), as part of the SIAR network, were used to apply the methodology developed in Simulink-MATLAB for the estimation of the global solar irradiation value on a surface with an inclination of 45◦ and oriented towards the equator. The results of the different methods of obtaining diffuse solar irradiation are shown below. for the estimation of the diffuse component. − Inclined Model 3 Inclined Model 3 used the anisotropic model of Temps and Coulson, presented in Equation (8), for the estimation of the diffuse component. − Inclined Model 4

of the different methods of obtaining diffuse solar irradiation are shown below. Figure 9 shows the values of the daily global solar irradiation obtained from the horizontal SIAR network, together with the calculations carried out for the estimation of the daily global solar irradiation on a surface inclined at 45° and oriented to the equator using Figure 9 shows the values of the daily global solar irradiation obtained from the horizontal SIAR network, together with the calculations carried out for the estimation of the daily global solar irradiation on a surface inclined at 45◦ and oriented to the equator using the following four models. Inclined Model 4 used the Klucher anisotropic model, presented in Equation (9), for the estimation of the diffuse component. − Inclined Model 5 Inclined Model 5 used Hay's anisotropic model, presented in Equation (10), to esti-

ance values, given in energy units (MJ/m2), on the inclined plane into the average hourly

In Figure 10, the hourly variation of the solar irradiance is represented for one day in April (spring) and another day in October (autumn) in order to compare the solar irradiance on the inclined plane calculated with the five conversion methods with the evolution of the data from the global horizontal solar irradiance SIAR located in Mansilla Mayor (León, Castile and León region, Spain), specifically the solar irradiance recorded by the pyranometer of the solar panel. Table 1 lists the statistical results for four random sunny days, comparing the estimates made with the isotropic model and the three anisotropic models with the measured values for the inclined plane oriented towards the equator.

values of the solar irradiance, given in power units (W/m2), on the inclined plane.

the following four models. − Inclined Model 1 mate the diffuse component.

− Inclined Model 2

− Inclined Model 1

**Figure 9.** Global solar irradiation measured on the horizontal surface during the year 2011 and estimated for a the solar panel inclined at 45° and oriented to the equator with Inclined Models 1, 2, 3, and 4. **Figure 9.** Global solar irradiation measured on the horizontal surface during the year 2011 and estimated for a the solar panel inclined at 45◦ and oriented to the equator with Inclined Models 1, 2, 3, and 4.

Inclined Model 1 used the CENSOLAR [48] method, which provides a coefficient to obtain the value of the global solar irradiation on the solar panel, depending on the latitude, the inclination, and the month of the year.

The following models (Inclined Models 2, 3, 4, and 5) used the methodology described in the previous sections for the values of the direct component and the albedo of the solar irradiation, while different models were used for the value of the diffuse component, therefore obtaining different results.

− Inclined Model 2

Inclined Model 2 used the Liu and Jordan isotropic model, presented in Equation (7), for the estimation of the diffuse component.

− Inclined Model 3

Inclined Model 3 used the anisotropic model of Temps and Coulson, presented in Equation (8), for the estimation of the diffuse component.

− Inclined Model 4

Inclined Model 4 used the Klucher anisotropic model, presented in Equation (9), for the estimation of the diffuse component.

− Inclined Model 5

Inclined Model 5 used Hay's anisotropic model, presented in Equation (10), to estimate the diffuse component.

Finally, the conversion module was applied to convert the global hourly solar irradiance values, given in energy units (MJ/m<sup>2</sup> ), on the inclined plane into the average hourly values of the solar irradiance, given in power units (W/m<sup>2</sup> ), on the inclined plane.

In Figure 10, the hourly variation of the solar irradiance is represented for one day in April (spring) and another day in October (autumn) in order to compare the solar irradiance on the inclined plane calculated with the five conversion methods with the evolution of the data from the global horizontal solar irradiance SIAR located in Mansilla Mayor (León, Castile and León region, Spain), specifically the solar irradiance recorded by the pyranometer of the solar panel. Table 1 lists the statistical results for four random sunny days, comparing the estimates made with the isotropic model and the three anisotropic models with the measured values for the inclined plane oriented towards the equator.

**Table 1.** Observed adjustment of the hourly solar irradiance on the plane inclined at 45◦ and oriented towards the equator, as estimated with Inclined Models 2, 3, 4, and 5 from data for the horizontal surface from the agrometeorological station SIAR in Mansilla Mayor (León, Castile and León, Spain), along with the values measured in León, for four days.


Inclined Model 2 utilized the Liu and Jordan isotropic model; Inclined Model 3 utilized the Temps and Coulson anisotropic model; Inclined Model 4 utilized the Klucher anisotropic model; Inclined Model 5 utilized the Hay anisotropic model; RMSE, root mean square error (W/m<sup>2</sup> ); R<sup>2</sup> , determination coefficient. The best results for each day are underlined.

11 April

11 October

**Figure 10.** Hourly solar irradiance obtained for a single day from the horizontal SIAR in Mansilla Mayor (León, Castile and León region, Spain) compared with that estimated for a solar panel inclined at 45° and oriented to the equator with Inclined Models 1, 2, 3, 4, and 5. (**Top**) 11 April 2011 (spring); (**bottom**) 11 October 2011 (autumn). **Figure 10.** Hourly solar irradiance obtained for a single day from the horizontal SIAR in Mansilla Mayor (León, Castile and León region, Spain) compared with that estimated for a solar panel inclined at 45◦ and oriented to the equator with Inclined Models 1, 2, 3, 4, and 5. (**Top**) 11 April 2011 (spring); (**bottom**) 11 October 2011 (autumn).

### **Table 1.** Observed adjustment of the hourly solar irradiance on the plane inclined at 45° and oriented towards the equator, as estimated with Inclined Models 2, 3, 4, and 5 from data for the horizontal surface from the agrometeorological station **4. Discussion**

SIAR in Mansilla Mayor (León, Castile and León, Spain), along with the values measured in León, for four days. **Day Inclined Model 2 Inclined Model 3 Inclined Model 4 Inclined Model 5**  Once the estimation of the solar irradiation received by the inclined solar panel has been undertaken, the results obtained were analyzed using the Simulink-MATLAB blocks (Figure 1).

**RMSE R2 RMSE R2 RMSE R2 RMSE R2** <sup>2011</sup>21.90 0.9751 25.09 0.9674 22.14 0.9746 17.83 0.9835 1 June 2011 38.51 0.8658 32.14 0.9065 31.54 0.9100 55.12 0.7250 30 June 2011 21.40 0.9415 40.11 0.7947 39.06 0.8052 11.95 0.9817 For the modeling of the extraterrestrial solar irradiation, symmetry can be observed throughout the day with regard to solar noon in Figure 2, along with a progressive increase in solar irradiation for the seasonal component from December to June and a decrease in solar irradiation from June to December (41.91-12.28 MJ/m<sup>2</sup> daily). The results were validated by comparing the sum of the hourly values of each day with those described by Allen [40] and resulted in a good approximation without the need for more statistics.

<sup>2011</sup>62.36 0.7853 31.44 0.9454 32.51 0.9416 37.32 0.9230 Inclined Model 2 utilized the Liu and Jordan isotropic model; Inclined Model 3 utilized the Temps and Coulson anisotropic model; Inclined Model 4 utilized the Klucher anisotropic model; Inclined Model 5 utilized the Hay anisotropic model; RMSE, root mean square error (W/m2); R2, determination coefficient. The best results for each day are underlined. **4. Discussion**  Three differentiated zones can be seen in the horizontal diffuse solar irradiation modeling (Figure 4). In one, a maximum hourly diffuse fraction (kd) and a minimum hourly clearness index (kt) were obtained (i.e., corresponding to the hours of the day with the sky covered); another zone was characterized by a minimum hourly diffuse fraction (kd) and a maximum hourly clearness index (kt) (i.e., corresponding to the hours of the day with clear skies); finally, there was a third intermediate zone with indices k<sup>d</sup> and kt inversely variable with each other, depending on the degree of cloudiness of the sky (i.e., corresponding to the hours of the partially cloudy day). The model provided in the study by Miguel et al. [43], which was developed in Mediterranean areas, shows a good approximation with the model obtained with the data recorded by the State Meteorological Agency (AEMet) at La Virgen del Camino station (León, Castile and León region, Spain) during the central eight hours of the day for the whole year of 2011. It can be observed that the data for clear sky days was different from the Miguel model (the mean diffuse fraction on clear days was lower than that of the Miguel model), which indicates a low dispersion of the solar irradiation of the sky on clear days (i.e., a very clear and clean sky in the place under study).

In the solar height modeling, in Figure 5 symmetry can be observed throughout the day with regard to solar noon, along with a progressive increase in solar height from December to June and a decrease in solar height from June to December. The results were validated by comparison with those described by Duffie et al. [41], showing a good approximation without the need for more statistics. The maximum values of solar height occurring at solar noon, which were observed throughout the year (Figure 6), reached their maximum during the last fortnight of June at the highest position of the Sun (72◦ ), and the minimum solar height occurred during the last fortnight of December (24.5◦ ).

For the results for the modeling of the angle of incidence, symmetry can be observed throughout the day with respect to solar noon, for which the minimum value was obtained, in Figure 7. In Figure 8, the progressive decrease from December to March (0.5◦ ), the increase from March to June (26.5◦ ), the decrease from June to September (4.5◦ ), and the increase from September to December (20.5◦ ) can be seen. The results were validated by comparison with those presented by Bérriz and Álvarez [49] and Duffie et al. [41], showing a good approximation without the need for more statistics. The most favorable months of the year, with the position of the solar panel inclined 45◦ and oriented towards the equator, for capturing solar energy were February, March, April, September, and October, as the inclined surface was in a more perpendicular position with regard to the path of the sun's rays on those dates.

The study of the variation of the angle of incidence is very important since it provides the time of year in which the use of solar energy is maximum, according to the inclination and orientation of the solar panel. The results of this model, if compared with the profile of annual energy needs of a photovoltaic system for a particular greenhouse, offer the possibility of finding the relationship (inclination–orientation) for each time of the year in which solar capture is optimal.

In M'Sila, Algeria, Ihaddadene et al. [50] theoretically searched for the best angle of inclination (monthly, seasonally, and annually), examining the Liu and Jordan model, the circumsolar model, the Hay model, and the Reindl model as the most appropriate models. They decided to change the angle of inclination of solar conversion systems monthly by 7272.27 MJ/m<sup>2</sup> , or seasonally by 7184.94 MJ/m<sup>2</sup> , instead of fixing them at 6836.83 MJ/m<sup>2</sup> in order to increase the amount of energy for a year; they recommended that this be studied using hourly data.

In Beijing, China, the annual optimal tilt angle shows a downward trend (Shen et al. [51]) compared to the optimum in the 1960s (i.e., 38◦ ); the optimal tilt angle decreased by 2◦ from 2011-2015 (i.e., to 36◦ ), caused mainly by the decrease in the direct irradiation ratio, which is highly related to atmospheric conditions (i.e., pollution that increases the proportion of diffuse irradiation and decreases the direct irradiation).

Kondratyev and Manolova [52] have stated that direct solar irradiation has been studied in detail, but diffuse and reflected solar irradiation on inclined planes are quite complicated, requiring study in relation to the distribution of the angular intensity. Therefore, special attention must be taken when selecting the most suitable diffuse irradiation model for a particular place; furthermore solar panels are more sensitive to electrical energy conversion in photovoltaic systems than solar collectors in thermal systems.

In light of these points, in this study, different models of diffuse irradiation were analyzed for a specific place where a greenhouse is located, taking into account that

the diffuse solar irradiation of the sky that falls upon the inclined plane of the solar panel is estimated in a different way as a function of the model. The circumsolar model overestimated it in Equation (6); the isotropic model underestimated it on the opposite slope of the inclined surface in Equation (7); and in clear and partly cloudy sky conditions, which occur in many cases, the light of the sky is anisotropic, meaning that the Temps and Coulson model provided a good prediction for the clear sky in Equation (8) but overestimated the solar irradiance when used for the cloudy sky. In this case, the Klucher model made corrections by setting factor F to 0 when the sky was completely cloudy in Equation (9), thus returning to the isotropic model, and setting F to 1 when the skies were, resulting in the Temps and Coulson model; this improved the estimates for all types of the sky.

In Nakhon Pathom and Ubon Ratchathani, Thailand, Wattan and Janjai [53] investigated the performances of 14 models in estimating hourly diffuse solar irradiation on inclined (30◦ , 60◦ , and 90◦ ) and oriented (north, south, east, and west) surfaces at two tropical sites, finding that the Muneer and Gueymard models performed better.

In Algiers and Ghardaia, Algeria, and Málaga, Spain, Takilalte et al. [54] estimated the inclined global irradiation in 5-min intervals, using global irradiation on the horizontal plane, geographical parameters, site albedo, and two cloudiness factors, based on a combination of two models (Perrin Brichambaut and Liu and Jordan) for which the parameters of the state of the sky transformed the isotropic models into an anisotropic model. The results of the proposed model for all sky conditions with regard to the normalized root mean square error (nRMSE), the relative percentage error (RPE), the normalized mean absolute error (nMAE), and the correlation coefficient (R<sup>2</sup> ) varied between 4.70–6.41%, 5.50–5.90%, 3.07–4.73%, and 0.97–0.99, respectively, which are very accurate results, especially for such short time scales in which there is no compensation or average effects, as occur when using monthly data.

Putting the three solar irradiation components together using Equation (11), with the solar panel at an inclination of 45◦ and oriented towards the equator, differences with regard to the incident on the horizontal surface shown in Figure 9 were observed according to the time of year:


In Athens, Greece, Raptis et al. [55] studied the ideal inclination to maximize the capture of solar irradiation, which was determined by the latitude and the time of year, with the horizontal surface receiving more irradiance than the inclined surface during the summer months and on cloudy winter days, due in this case to an anisotropy of the diffuse light, with a greater diffuse contribution coming from angles closer to the zenith; however, the inclined surface reached higher values than the horizontal one in winter, with the optimum angle found to be around 30◦ during the year.

The result of the comparison of the recorded measurements of solar irradiance on the horizontal surface (i.e., the SIAR data) and the inclined plane of the solar panel (i.e., the pyranometer in León) with the results obtained with Inclined Models 1, 2, 3, 4, and 5 (Figure 10) show a better approximation with the anisotropic models for four random sunny days. Specifically: Inclined Model 5, which used the Hay corrections, obtained better results on 4 November 2011 and 30 June 2011 with regard to the RMSE (17.83 and 11.95 W/m<sup>2</sup> ) and R<sup>2</sup> (0.9835 and 0.9817), respectively. Inclined Model 4, which used the Klucher corrections, performed best on 6 January 2011, with RMSE = 31.54 W/m<sup>2</sup> and R <sup>2</sup> = 0.9100. Inclined Model 3, which used Temps and Coulson corrections, obtained better results on 10 November 2011, with RMSE = 31.44 W/m<sup>2</sup> and R<sup>2</sup> = 0.9454.

The solar panel global solar irradiation for 2011 was 1.84 MWh/m<sup>2</sup> with the CEN-SOLAR model, which used coefficients; 1.84 MWh/m<sup>2</sup> for the Liu and Jordan model; 1.99 MWh/m<sup>2</sup> for the Temps and Coulson model; 1.93 MWh/m<sup>2</sup> for the Klucher model; and 1.86 MWh/m<sup>2</sup> for the Hay model.

The global solar irradiation for the horizontal SIAR measured in Mansilla Mayor (León, Castile and León region, Spain) during 2011 was 1.67 MWh/m<sup>2</sup> , thus resulting in a higher value on the solar panel using corrections of 10.17% for the model by Liu and Jordan, 19.16% for the Temps and Coulson model, 15.56% for the Klucher model, and 11.37% for the Hay model. However, this was distributed throughout the year, as shown in Figure 9, and as mentioned previously, was higher in the months with moderate solar irradiation (i.e., February, March, April, September, October, and November), moderate in high solar irradiation months (i.e., May and August), and lower in the months with very high solar irradiation (i.e., June and July).

In Adrar, Algeria, an increment in the performance of horizontally placed solar panels was achieved by Bailek et al. [56] with a fixed inclination of 20.61% monthly, 19.58% seasonally, 19.24% semi-annually, and 13.78% for annual adjustments, with the optimal tilt angles in each period.

### **5. Conclusions**

Glasshouses are agricultural productive structures intended to increase the production and quality of early bloomer crops. They can be energetically characterized as follows:


Thus, distributed generation PV systems are ideal for connection to glasshouses, either on their own or together with power generators where the value of the solar irradiance which falls upon the solar panels is the main variable to determine the performance of the PV system.

The literature pertaining to the estimation of the incidental solar irradiance on the solar panel plane (at an angle and/or oriented) in relation to the irradiance received on the horizontal surface (data registered in meteorological stations) is highly diverse, especially with regard to the types of sky in particular places. However, the practical use of this diverse information is complex, which is what incentivized the present work, in which the following methodology was adopted.


which may have been due to some cloudiness or changes in the reflections of the surrounding light.

As a final conclusion to this paper, it should be noted that the solar estimation for the inclined plane can be used together with the daily prediction of solar irradiance, as detailed by Diez et al. [57], to obtain the value of available solar energy in the glasshouse PV system and thus to enable more efficient management of the glasshouse electric demand.

**Author Contributions:** Conceptualization, F.J.D., L.M.N.-G. and L.C.-S.; Methodology, F.J.D., L.M.N.-G. and A.M.-R.; Software, F.J.D., A.M.-R. and R.A.; Validation, F.J.D., L.C.-S. and A.C.-G.; Formal analysis, L.M.N.-G., A.M.-R. and R.A.; Investigation, F.J.D., L.M.N.-G. and L.C.-S.; Resources, L.M.N.-G. and A.C.-G.; Writing—original draft preparation, F.J.D., A.M.-R. and R.A.; Writing—review and editing, L.M.N.-G., L.C.-S., A.C.-G. and R.A.; Visualization, F.J.D. and A.M.-R.; Supervision, F.J.D., A.C.-G. and R.A.; Project administration, L.M.N.-G. and A.M.-R.; Funding acquisition, L.M.N.-G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors wish to acknowledge CYTED (the Ibero-American Program of Science and Technology for Development) for supporting this work through collaboration with the RITMUS network.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

### **References**

