In this study, an algorithm, based on which it is possible to estimate the total shading caused by photovoltaic modules placed on the roof of an even-span greenhouse is presented. The algorithm is essentially a parent function, consisting of five nested functions, as well as three sections of individual commands for the number and position of photovoltaics on the roof. The final extracted results of the algorithm concern the calculation of the shaded surface by the photovoltaics within the desired construction unit of the greenhouse, the percentage of shading on the total surface of the construction unit, as well as two graphs, the first with the 3D visualization of the greenhouse, the photovoltaics and of the formed shadow in a Cartesian coordinate system and the second with the top view of all of the above in a 2D Cartesian coordinate system. The estimation of the shadow may be done up to a time step, t, equal to 10 minutes at any time of year, with the algorithm incorporating criteria that prohibit it from operating for intermediate time intervals or circumstances where it does not exist on a calendar basis, but also geographically. The above criteria are presented in Table 1. All the procedures used for the functions are outlined below, along with the functions as they were implemented using the MATLAB programming language.
Table 1
Calendar basis and geographical criteria in the algorithm.
| Variable | Criterion | MATLAB Code |
Date/Time | Days | Day = [1:1:last day of month] | if d < 1 || d > eomday(yr,mth) error(); end |
Hours | Hour = [0:1:23] | if hr < 0 || hr > 23 error(); end |
Minutes | Minutes = [0:10:50] | if mnt~=[0:10:50] error(); end |
Seconds | Seconds = 0 | if sc ~ = 0 error(); end |
Location | Longitude | Longitude = [0,360] | if lon < 0 || lon > 360 error(); end |
Latitude | Latitude = [-90,90] | if lat<-90 || lat > 90 error(); end |
The nested functions implemented, concern the Sun’s position, as it can be calculated in a spherical coordinate system (Function 1), the distance between the Earth and the Sun (Function 2), the results of which are used to transform the Sun’s position from spherical to Cartesian coordinates (Function 3). Then, the coordinates of the points that form the greenhouse are found in a Cartesian coordinate system (Function 4) based on greenhouse basic characteristics, such as the number of construction units, which are repeated either widthwise or lengthwise, the length of the construction unit, its width, the gutter, and ridge height. The last nested function created and used is finding the shadow formed by the photovoltaics (Function 5).
The individual sections of the algorithm implemented concern the coordinates of the points that correspond to the 4 corners of each photovoltaic unit, depending on their position on the roof of the greenhouse, both for the first and the second unit in which they have placed the photovoltaics (Section 1), the binary representation of the surface covered by the greenhouse, and the greenhouse’s ground surface shaded by the photovoltaics (Section 2), and finally the calculation of the shaded area within the arable area of the greenhouse (Section 3). In Fig. 1, the order in which Functions and Sections appear in the algorithm is presented.
Although Functions 4 and 5 can be used for any case of a greenhouse and the location of photovoltaics on its roof, their use in this research concerns a specific case of a greenhouse and photovoltaics. The greenhouse used in the research is located on the campus of the University of Patras (38°17′27.9′′ N, 21°47′23.9′′ E). It is a real-scale, even-span greenhouse with an East-West ridge. The East-West orientation of the greenhouse is intended to be physically identical to that of productive greenhouses. The greenhouse consists of four different units, of which only two are considered in this work, the unit north of the greenhouse (northern construction unit) and the unit next to it (southern construction unit). The northern construction unit is also one available for cultivation. Into the south-sloped level of the roof of each unit, a total of 12 photovoltaic modules have been integrated. The southern placement has been made for greater energy production, especially in winter, when the sun is at a low altitude. More specifically, 8 of the 12 photovoltaic modules have been integrated into the roof of the northern construction unit, while 4 of them are on the roof of the unit next to it. Assuming that 8 positions are available for the photovoltaics on the roof of the greenhouse, the positions covered in the northern construction unit are 1 to 8 (Fig. 2a), while for the southern greenhouse unit, the positions covered are 3,4 (near the west side of the greenhouse) and 7,8 (near the east side of the greenhouse) (Fig. 2b). Table 2 presents the characteristics of the greenhouse construction units, as well as the dimensions of the photovoltaics. It should be mentioned that although the photovoltaics that have been integrated into the greenhouse roof are semi-transparent, in the context of this study they have been considered as completely opaque.
Table 2
Greenhouse construction unit and PV modules characteristics.
Greenhouse construction unit characteristics | Values |
Orientation | East – West |
Number of construction units – widthwise | 2 |
Number of construction units – lengthwise | 0 |
Width | 3.20 m |
Length | 16.39 m |
Gutter height | 3.22 m |
Ridge height | 3.90 m |
Photovoltaic modules dimensions | |
Width | 1.033 m |
Length | 2.048 m |
2.1 Calculations for the position of the sun
Starting with one of the most fundamental factors, the location of the sun in the sky can be defined by two separate angles in a spherical coordinate system, the solar zenith angle – SZA, and the solar azimuth angle – SAA (Marco and Giuseppe 2018). The SZA is defined as the angle formed by the line passing between the sun and the study area on Earth, and the perpendicular line to the same area (the z-axis in a Cartesian coordinate system). It can also be defined as the angle at which the beam radiation strikes the earth. The SAA is defined as the angle formed by the location of the meridian and the line projecting the sun to the observer on the horizontal plane and can range from − 180° to 180°. The angle can range between 90° and − 90° in mid-latitude locations and on days less than 12 hours. However, on days longer than 12 hours, the angle may be more than 90° or less than − 90°, depending on whether it is approaching sunset or sunrise. The SAA can be calculated based on two different systems, depending on the position that takes its zero value. In the first system, where its zero value is taken towards the South, the value of the azimuth angle becomes increasingly negative as the observer approaches the East, while its value becomes increasingly positive approaching the West. A more practical method is to base the angle's sign on the sign of the solar hour angle, with the azimuth angle being positive when the hour angle is positive and negative when the hour angle is negative. The calculation of the solar zenith angle is presented in Eq. 1, and according to the above solar azimuth angle is calculated by Eq. 2 (Duffie et al. 2020):
$$SZA={\text{cos}}^{-1}\left[\text{cos}\left(\phi \right)\bullet \text{cos}\left(\delta \right)\bullet \text{cos}\left(\omega \right)+\text{sin}\left(\phi \right)\bullet \text{sin}\left(\delta \right)\right]$$
1
$$SAA=sign\left(\omega \right)\bullet \left|{\text{cos}}^{-1}\left[\frac{\text{cos}\left(SZA\right)\bullet \text{sin}\left(\phi \right)-\text{sin}\left(\delta \right)}{\text{sin}\left(SZA\right)\bullet \text{cos}\left(\phi \right)}\right]\right|$$
2
where SZA is the solar zenith angle (°), SAA is the solar azimuth angle (°), φ is the latitude of the area under study (°, positive towards the Northern Hemisphere), δ is the solar declination (rad), and ω is the solar hour angle (°). In this study, the SAA values used are based on the aforementioned second system, where its zero value is taken to-ward the North and increases as we move to the East. In such a system, the SAA takes only positive values in a range between 0° and 360°. The conversion of the SAA from system 1 to system 2 was achieved by adding to the first, value of 180°.
The calculation of the aforementioned angles (SZA & SAA) requires the calculation of additional parameters, such as the solar declination – δ, and the solar hour angle – ω which are characterized by increased complexity due to the desired short timestep in the present study, which is equal to 10 minutes. Starting with the calculation of the fractional year – γ – the expression of a given period into decimal format is achieved. The fractional year is calculated by Eq. 3:
$$\gamma =\frac{2\bullet \pi }{365}\bullet \left(N-1+\frac{h-12}{24}\right)$$
3
where γ is the fractional year (rad), N is the day of the year, which ranges between 1 and 365 or 366 (for leap years), and h is the needed hour of the day. Through the fractional year, it is possible to calculate both the equation of time and the solar declination with a 10-minute timestep. The equation of time – Et considers the perturbations in the earth’s rotation rate to determine when the sun crosses the observer’s meridian (Duffie et al. 2020). For the calculation of the equation of time, Eq. 4 (Duffie et al. 2020) was used:
$${E}_{t}=229.2\bullet \left[0.000075+0.001868\bullet \text{cos}\left(\gamma \right)-0.032077\bullet \text{sin}\left(\gamma \right)-0.014615\bullet \text{cos}\left(2\bullet \gamma \right)-0.04089\bullet \text{sin}\left(2\bullet \gamma \right)\right]$$
4
where Et is the equation of time (min), and γ is the fractional year (rad). In Eq. 4, the equation of time is represented by the term in square brackets on the right-hand side of the equation, while the multiplier value 229.2 converts it into minutes. The results obtained using Eq. 4 are accurate to within 0.0025 rad, which is approximately equivalent to 35 seconds (Iqbal 1983). The solar declination – δ refers to the sun’s angle with respect to the equator during solar noon, when it is in the north positive direction and is on the local meridian. δ takes values in a range of -23.45° to 23.45° and is calculated by Eq. 5 (Duffie et al. 2020):
$$\delta =0.006918-0.399912\bullet \text{cos}\left(\gamma \right)+0.070257\bullet \text{sin}\left(\gamma \right)-0.006758\bullet \text{cos}\left(2\bullet \gamma \right)+0.000907\bullet \text{sin}\left(2\bullet \gamma \right)-0.002697\bullet \text{cos}\left(3\bullet \gamma \right)+0.00148\bullet \text{sin}\left(3\bullet \gamma \right)$$
5
where δ is the solar declination (rad), and γ is the fractional year (rad). The results obtained by Eq. 5 are accurate to within 0.0006 rad (Spencer 1971).
Another important parameter that needs to be calculated is the true solar time – TST. The time employed in all sun-angle connections is known as solar time, which is different from local clock time. Standard time must be converted into true solar time, which can be accomplished by applying two modifications. The first modification concerns the difference that is presented in the longitude between the observer’s meridian and the meridian that serves as the foundation for local standard time – LST. One degree of longitude is traversed by the sun in 4 minutes. The second adjustment comes from the equation of time. The TST is given by Eq. 6:
where TST is the true solar time (min), LST is the local solar time (min) and time offset is the Coordinated Universal Time (UTC) offset, which represents the difference in time, measured in minutes, between a location’s LST and UTC. The time offset is calculated by Eq. 7 (Duffie et al. 2020):
$$time offset=4\bullet \left({L}_{st}-{L}_{loc}\right)+{E}_{t}$$
7
where Lst is the standard meridian for the local time zone, Lloc is the longitude of the lo-cation under study (both in degrees and positive west of the Prime Meridian, so that 0⁰ < L < 360⁰), and Et is the equation of time given by Eq. 4. Lst can be calculated by multiplying the time zone of the study area by 15 degrees. This is because the solar hour angle varies at a rate of 15 degrees every hour (Wang 2019).
The final parameter that must be calculated is the solar hour angle – ω. The solar hour angle is the sun's angular displacement from the local meridian east or west, as a result of the earth's 15° per hour rotation on its axis. Hour angle can have negative and positive values in the morning and the afternoon, respectively, and can be calculated by Eq. 8 (Duffie et al. 2020):
$$\omega =\frac{TST}{4}-180^\circ$$
8
where TST is the true solar time (min).
Function 1 (SunPos) provides the main points of the solar position computation, according to Equations 1–8.
Function 1: Sun Position in Spherical Coordinates |
function [SZA,SAA] =… SunPos(yr,mth,d,hr,mnt,sc,l_loc,lat,timezone) dt = datetime(yr,mth,d,hr,mnt,sc,l_loc,lat,timezone) γ← Eq. 3 Et ← Eq. 4 δ← Eq. 5 ω← Equations 6–8 SZA ← Eq. 1 SAA ← Eq. 2 for k = 1:length(SAA) if ω(k) > 0, SAA(k) = SAA(k); elseif ω(k) < 0, SAA(k) = -SAA(k); end end end |
2.2 Transformation from Spherical to Cartesian Coordinate System
To identify the lines on which both the boundary points of the photovoltaics and the position of the sun will correspond, it is necessary that the second one, which is calculated in a spherical coordinate system S(r,θ,φ), must be transformed into a Cartesian coordinate system C(x,y,z), where r is the radial distance, or in this case the distance between the point (0,0,0) on Earth and the Sun’s position, θ is the polar angle or the solar zenith angle – SZA, and φ is the azimuthal angle or the solar azimuth angle – SAA. The variables x, y, and z are defined as the coordinates of any point on the Cartesian coordinate system, with respect to the x-, y-, and z-axis, respectively.
The system of the equations used to implement the above transformation of the sun's position into the Cartesian coordinate system is presented by Equations 9, 10, and 11.
$${x}_{sun}=r\bullet \text{sin}\left(SZA\right)\bullet \text{cos}\left(SAA\right)$$
9
$${y}_{sun}=r\bullet {sin}\left(SZA\right)\bullet \text{sin}\left(SAA\right)$$
10
$${z}_{sun}=r\bullet {sin}\left(SZA\right)$$
11
where xsun, ysun, and zsun are the coordinates of the Sun’s position with respect to the x-, y-, and z-axis, respectively.
One parameter from the above, which also shows a dependence on time, albeit with relatively minor variation in such a small timestep, is the Earth-Sun’s distance, r. For the calculation of the Earth-Sun distance, the Kepler orbit theory was used. A Kepler orbit can be described as the motion of a body in relation to another. Such an orbit forms a two-dimensional curve (ellipse, parabola, or hyperbola), in a three-dimensional space, and is characterized by elements such as the semimajor axis, α, eccentricity, e, and the inclination, i. In the Earth-Sun system, the Earth moves in an elliptical orbit (0 < e < 1) around the sun, while there are distances, such as aphelion and perihelion, which correspond to the shortest and the greatest distance of the Earth from the Sun, respectively, and their values have been identified. Taking into consideration the above information, the Earth-Sun distance calculation is feasible by using specific equations, even for very short periods (Krüger and Grün 2014; Pälike 2005). In Table 3, the characteristics of Earth’s elliptical orbit that were used for the calculation of r, are presented.
Table 3
Characteristics of Earth’s elliptical orbit around the Sun.
Element | Value |
Eccentricity | 0.01671 |
Perihelion distance | 0.983 AU or 147,054,706,898 m |
Based on the relation, which connects the Earth-Sun’s distance at the perihelion to the semimajor axis, α, and the eccentricity, and basic general equations for each curve in the form of an ellipse, the parameters α (Eq. 12), c (Eq. 13) and b (Eq. 14) are calculated (Krüger and Grün 2014).
$$b=\sqrt{{a}^{2}-{c}^{2}}$$
14
where a is the semimajor axis (m), p is the Earth-Sun’s distance at the perihelion (m), e is the eccentricity (-), c is the linear eccentricity or the distance between the center of the ellipse and the focal (m), and b is the semi-minor axis of the ellipse (m). Assuming that the above ellipse is centered at the point C(0,0) (Fig. 3) in a two-dimensional Cartesian coordinate system, and the Sun is on the foci 1 of the ellipse (point F1 in Fig. 3), the parametric equation that gives the earth's orbit around the sun is given by Eq. 15, with the orbital coordinates given by Eq. 16 for the x-component and from Eq. 17 for y-component (Krüger and Grün 2014; Xu and Xu 2013). Eq. 18 presents the parametric equation of the distance between foci 1 and the point (x,y), or in the Earth-Sun system, the distance between Earth and Sun. At the same time, in Fig. 3, the position of the perihelion and the aphelion are presented by points P and A, respectively.
$$\frac{{x}^{2}}{{a}^{2}}+\frac{{y}^{2}}{{b}^{2}}=1$$
15
$$x=a\bullet \text{cos}\left(t\right)$$
16
$$y=b\bullet \text{sin}\left(t\right)$$
17
$$r=\sqrt{{\left(x-c\right)}^{2}+{y}^{2}}$$
18
where a is the semi-major axis of the ellipse, b is the semi-minor axis of the ellipse, and t is an angle that varies in a range from 0 to 360°.
For the estimation of the angle t, it is considered that the Earth requires 365 days for a complete rotation of 360°, therefore for each day of the year, the Earth traverses an angle equal to 0.986°, while for each 10-minute period, which is the desired time step in this study, this angle takes a value equal to 0.986°/144, as one day corresponds to 144 10-minute time intervals. Finally, the parametric equation resulting from the above and used to calculate the Earth-Sun distance (distance F1E in Fig. 3) for every 10 minutes of the year is given by Eq. 18, which by substituting Equations 16 and 17 turns into Eq. 19.
$$r\left({t}^{{\prime }}\right)=\sqrt{{\left[a\bullet \text{cos}\left(\frac{0.986}{144}\bullet {t}^{{\prime }}\right)-c\right]}^{2}+{\left[b\bullet \text{sin}\left(\frac{0.986}{144}\bullet {t}^{{\prime }}\right)\right]}^{2}}$$
19
where r is the Earth-Sun distance for 10 minutes around the year (m), and t’ is the index of the 10-minute period (-) throughout the year.
In Function 2 (EarthSunDist) the main points of the Earth-Sun distance computation are presented, according to Equations 12–14 and 19.
Function 2: Earth-Sun distance |
function r = EarthSunDist(yr,mth,d,hr,mnt,sc) dt_year = datetime(yr,1,1,0,0,0):minutes(10):… datetime(yr,12,31,23,50,0); ten_min_idx =… (1:1:find(dt_year = = datetime(yr,12,31,23,50,0)))’; dt_inp = datetime(yr,mth,d,hr,mnt,sc); ten_min = ten_min_idx(dt_inp = = dt_year); p = 0.983*149597870700; % Perihelion in meters a ← Eq. 12 c ← Eq. 13 b ← Eq. 14 r ← Eq. 19 end |
The calculated Earth-Sun distance by Function 2 is used in conjunction with the SZA and SAA angles calculated from Function 1 to transform the spherical coordinates of the Sun's position to Cartesian, which is presented by Function 3 (SunPosCart). For the times when the solar altitude angle, which is referred to as the complementary of the SZA, is less than zero, thus describing the nighttime hours, the coordinates for the position of the sun do not exist.
Function 3: Spherical to Cartesian Coordinates |
function [x_sun,y_sun,z_sun] = SunPosCart(r,SZA,SAA) x_sun ← Eq. 9 y_sun ← Eq. 10 z_sun ← Eq. 11 x_sun((90 - SZA) < 0) = NaN; y_sun((90 - SZA) < 0) = NaN; z_sun((90 - SZA) < 0) = NaN; end |
2.3 Greenhouse coordinates in the Cartesian coordinate system
The fourth function (Function 4) concerns the basis of the whole study, the greenhouse. This specific function gives the possibility to map an even-span greenhouse and find the points that form it in a three-dimensional Cartesian coordinate system. The orientation of the system (East-West-North-South) was determined based on the results of Function 1, where the position of the Sun is estimated. As a result, the North-South and the East-West direction, are represented by the x-axis and the y-axis, respectively, while the altitude from the surface is represented by the z-axis. Function 4 receives as inputs, basic parameters concerning a basic greenhouse construction unit, such as the orientation (ridge direction), the width (short-side length), the length (long-side length), and the gutter and ridge heights from the ground surface. Simultaneously, the repeatability of the basic construction unit to create the overall greenhouse, either lengthwise or widthwise, is considered. All the inputs of Function 4 are presented in Table 4.
Table 4
Inputs of greenhouse characteristics in Function 4.
Greenhouse Characteristics | Function 4 inputs |
Orientation | dir |
Number of construction units – widthwise | numCU_wdwise |
Number of construction units – lengthwise | numCU_lnwise |
Width | grh_width |
Length | grh_length |
Gutter height | gut_hght |
Ridge height | ridge_hght |
To begin with, it is necessary to define the reference point, based on which the greenhouse will be formed. This is regarded the origin of the axes (x,y,z) = (0,0,0) in this function, which corresponds to the bottom southwest corner of the greenhouse. Based on this point, all the other points that represent the edges of the frame of the first construction unit of the greenhouse are defined. Once these points have been calculated, they are stored in a two-line cell-type variable in MATLAB, the first one corresponds to the greenhouse units in the case of widthwise repentance, while the second one is to the greenhouse units in the case of lengthwise repentance. The size of the cell varies according to the number of construction units. Finally, each element of the cell consists of an array of three columns, one for each direction (x,y,z), and ten rows, one for each point of the greenhouse.
Function 4 (GreenhouseCoord) contains two limitations, the first concerns the ridge orientation, for which there are only two options, that of the East-West direction (Case 1) and that of the North-South direction (Case 2). The second limitation concerns the fact that the number of construction units should be a positive integer number, but also the definition of at least one construction unit. For the detection of the greenhouse frame points, except for the two previously stated cases, it was essential for additional sub-cases which involve scenarios for the number of construction units with widthwise or lengthwise repetition. More specifically, these cases describe:
1. One or more construction units widthwise and none or one unit lengthwise (Sub-case 1). The situation of one unit in length is identical to the first unit in width, but it must be defined due to a different variable in the Function’s input.
2. One or more construction units lengthwise and none or one unit widthwise (Sub-case 2)
3. More than one construction unit both widthwise and lengthwise (Sub-case 3).
Finally, three different arrays are created for the points, one for each coordinate based on the characteristic parameters of the greenhouse. Each of them consists of ten lines, one for each point and they are the output of the nested function “GrhXYZ”. The use of this nested function shortens and speeds up the code. The repeatability of units for the widthwise and lengthwise cases is achieved through a “for loop” for the x- and y-coordinate, respectively for each sub-case. The “for loops” for each sub-case are presented in Table 5. In contrast, in the case of the North-South direction of the ridge, the widthwise repetition concerns the y-coordinate, while the lengthwise repetition concerns the x-coordinate. Function 4 is presented only as the case for an East-West oriented greenhouse.
Table 5
Sub-cases for the number of construction units with widthwise or lengthwise repetition.
Sub-case | Definition code | “for loops” code |
1 | numCU_wdwise > 0 && (numCU_lnwise = = 0 || numCU_lnwise = = 1) | for k = 0:numCU_wdwise-1 for l = 0 |
2 | numCU_lnwise > 0 && (numCU_wdwise = = 0 || numCU_wdwise = = 1) | for k = 0 for l=… 0:numCU_lnwise-1 |
3 | numCU_wdwise > 1 && numCU_lnwise > 1 | for k = 0:numCU_wdwise-1 for l=… 90:numCU_lnwise-1 |
Function 4: Greenhouse coordinates in the Cartesian coordinate system |
function grh_points =… GreenhousCoord(dir,numCU_wdwise,numCU_lnwise,… grh_width,grh_length,gut_hght,rigde_hght) % limitation 1 assert(any(strcmp(dir,{'NS','EW'})),… 'Greenhouse direction must be either North-… South or East-West') % limitation 2 function validate_positive_integer_of_CU(direction,… n9umCU) if numCU < 0||abs(numCU-round(numCU)) > 0.0001,… error('%s must be a positive integer',… direction); end end validate_positive_integer_of_CU('Number of… costrucion units',numCU_wdwise); validate_positive_integer_of_CU('Number of… construction units', numCU_lnwise); if numCU_wdwise = = 0 && numCU_lnwise = = 0,… error('Must define at least one construction… unit.'); end function grh_points = GrhXYZ(dir,k,l,grh_width,… grh_length,gut_hght,rigde_hght) grh_z_coord=[0,0,0,0,gut_hght,gut_hght,… gut_hght,gut_hght,rigde_hght,rigde_hght]'; % Case 1 if dir=='EW' grh_x_coord=[grh_width*k,grh_width*(k + 1),… grh_width*k,grh_width*(k + 1),grh_width*k,… grh_width*(k + 1),grh_width*k,… grh_width*(k + 1),(grh_width/2)*(2*k + 1),… (grh_width/2)*(2*k + 1)]'; grh_y_coord=[grh_length*l,grh_length*l,… grh_length*(l + 1),grh_length*(l + 1),… grh_length*l,grh_length*l,grh_length*l,… grh_length*(l + 1),grh_length*(l + 1),… grh_length*l,grh_length*(l + 1)]'; grh_points=… table(grh_x_coord,grh_y_coord,grh_z_coord); end end grh_points{l + 1,k + 1} ← Sub-Case 1 grh_points{l + 1,k + 1} ← Sub-Case 2 grh_points{l + 1,k + 1} ← Sub-Case 3 end |
2.4 Calculation of shadow coordinates
The points that form the shaded surface by the photovoltaics are located by Function 5. To find this surface on the ground, it is necessary to find the points that form the photovoltaics in the three-dimensional Cartesian coordinate system (xPV,yPV,zPV), as well as the corresponding point of the sun's position (xsun,ysun,zsun), as computed by Functions 1–3. In Section 1 of the algorithm (Calculation of PV Points) the points of the photovoltaic location are calculated according to the Cartesian coordinate system.
Firstly, in Section 1, the positions of the photovoltaics installed in each construction unit are defined. Based on the length of the photovoltaics and the length of the greenhouse, the total available positions that the photovoltaics can take on the greenhouse roof are eight. However, some limitations have been set, based on which the algorithm stops, giving the corresponding error messages. These limitations correspond to cases where the total length of the photovoltaic array exceeds the total length of the greenhouse roof (Limitation 1), as well as cases where the width of the photovoltaics exceeds the length of the inclined plane of the roof (Limitation 2). The definition codes of these limitations are presented in Table 6.
Table 6
Limitations for the possible photovoltaic positions.
Limitation | Definition code |
1 | if PV_width>(sqrt((rigde_hght-ut_hght).^2… + (grh_width).^2)) || PV_width < = 0, error("The width of the PV Unit must be a… positive integer and not exceed the… length of the inclined plane of the… roof"); end |
2 | if grh_length < numel(pos_PV(n))*PV_length, error("The PV units exceed the limits of… the greenhouse roof"); end |
At the same time, in Section 1 two more cases are distinguished, one for the placement of the photovoltaics in a greenhouse with an East-West direction and the other for the placement of the photovoltaics in a greenhouse with a North-South direction. These cases are separated due to the repetition required to find the coordinates of the PVs. In the first case, a repetition is required for the x-coordinate, as it differs between the construction units, since they are repeated parallel to the x-axis. Simultaneously, with the placement of the photovoltaics done with their long side parallel to the ridge, the y-coordinate changes, with its value being calculated through repetition based on the length of the photovoltaic and the position where each element is placed. The y-coordinate for the case of the different units in the case of the East-West direction remains constant. For the second case, where the greenhouse has a North-South direction, the y-coordinate performs like the x-coordinate of the previous case, with it differing be-tween construction units, since now the greenhouse units are repeated parallel to the y-axis. Accordingly, the value of the x-coordinate is calculated through a repetition based on the length and position of the photovoltaic.
With the photovoltaics on the roof of the greenhouse placed so that their long side is the same as the ridge, some of their points are identical to the already calculated points of the greenhouse. The coordinates of the remaining points have been found based on simple geometric relationships and the characteristics of the greenhouse. Finally, the z-coordinate appears the same for both the first and the second case, as both the gutter and ridge height remains constant.
In Section 1, the code is presented only for the case of a greenhouse with an East-West orientation. For the 3D visualization of the photovoltaics on the roof of the greenhouse, it is necessary to create data grids based on the corner points for each photovoltaic. The main commands used are "linspace" and "meshgrid", while the "CU_tables" is a cell-array that contains the greenhouse points that have been extracted by Function 4.
Section 1: Calculation of PV Points |
pos_PV={pos_PV_CU_1,pos_PV_CU_2}; if grh_dir=='EW' for i = 1:size(CU_tables,2) for m = 1:size(CU_tables,1) for j = pos_PV{m,i} PV_x(:,j,i,m)=… [(cosd(atand(2*(rigde_hght-gut_hght)/grh_width))*… (sqrt(((rigde_hght-gut_hght).^2)+… ((grh_width.^2)/4))-PV_width))+(i-1)*grh_width,… CU_tables{1,i}.x("Point 9"),… (cosd(atand(2*(rigde_hght-gut_hght)/grh_width))*… (sqrt(((rigde_hght-gut_hght).^2)+… ((grh_width.^2)/4))-PV_width))+(i-1)*grh_width,… CU_tables{1,i}.x("Point 9")]’; PV_y(:,j,i,m)=… [PV_length*(j-1)+(m-1)*grh_length,… PV_length*(j-1) + (m-1)*grh_length,… PV_length*j + (m-1)*grh_length,… PV_length*j + (m-1)*grh_length]'; PV_z(:,j,i,m)=… [(sind(atand(2*(rigde_hght-gut_hght)/grh_width))*… (sqrt(((rigde_hght-gut_hght).^2)+… ((grh_width.^2)/4))-PV_width)) + gut_hght,… CU_tables{1,i}.z("Point 9"),… (sind(atand(2*(rigde_hght-gut_hght)/grh_width))*… (sqrt(((rigde_hght-gut_hght).^2)+… ((grh_width.^2)/4))-PV_width)) + gut_hght,… CU_tables{1,i}.z("Point 9")]’; PV_Points{:,j,i,m}=… [PV_x(:,j,i,m),PV_y(:,j,i,m),PV_z(:,j,i,m)]; PV_Points {m,i}=PV_Points(:,:,i,m); % PV visualization X1{m,i}=PV_x_coord(:,:,i,m); Y1{m,i}=PV_y_coord(:,:,i,m); Z1{m,i}=PV_z_coord(:,:,i,m); Xi{m,i}{j}=… linspace(min(X1{m,i}{:,j}),max(X1{m,i}{:,j}),2); Yi{m,i}{j}=… linspace(min(Y1{m,i}{:,j}),max(Y1{m,i}{:,j}),2); ZiX{m,i}{j}=… linspace(min(Z1{m,i}{j}),max(Z1{m,i}{j}),2); ZiY{m,i}{j}=… linspace(min(Z1{m,i}{j}),max(Z1{m,i}{j}),2)'; [Xii{m,i}{j},Yii{m,i}{j}]=… meshgrid(Xi{m,i}{j},Yi{m,i}{j}); Zii{m,i}{j}=meshgrid(ZiX{m,i}{j},ZiY{m,i}{j}); end end end end |
Based on the points related to the four corners of the photovoltaic module and the point related to the sun’s position, four lines are formed for each photovoltaic module. Each of these lines that are formed for each desired time t corresponds to one of the four photovoltaic points, while a common point is the position of the sun. For each of these four lines, there is a three-dimensional vector of the form (Eq. 20):
$$\overrightarrow{{r}_{sh.point,n}}=\left(\overrightarrow{{x}_{sh,n}},\overrightarrow{{y}_{sh,n}},\overrightarrow{{z}_{sh,n}}\right)=\left[{x}_{n}+k\bullet \left({x}_{s}-{x}_{n}\right),{y}_{n}+k\bullet \left({y}_{s}-{y}_{n}\right),{z}_{n}+k\bullet \left({z}_{s}-{z}_{n}\right)\right]$$
20
where \(\overrightarrow{{r}_{sh.point,n}}\) is the position vector of the shadow point produced by one of the four points of the PV, \(\overrightarrow{{x}_{sh,n}},\overrightarrow{{y}_{sh,n}},\overrightarrow{{z}_{sh,n}}\) are the position vectors of the shadow point with respect to each axis of the three-dimensional coordinate system, xs, ys, and zs are the coordinates that define the position of the sun, xn, yn, and zn, are the coordinates defining each point of the photovoltaic, with n = 1, 2, 3, 4, and k is a constant.
From Eq. 20, a system of three partial equations that represent the coordinates of the points that form the shadow, is obtained (Equations 21, 22, and 23).
$${x}_{sh,n}\left(k\right)={x}_{n}+k\bullet \left({x}_{s}-{x}_{n}\right)$$
21
$${y}_{sh,n}\left(k\right)={y}_{n}+k\bullet \left({y}_{s}-{y}_{n}\right)$$
22
$${z}_{sh,n}\left(k\right)={z}_{n}+k\bullet \left({z}_{s}-{z}_{n}\right)$$
23
Solving the above system of equations and setting the z-coordinate (zsh) equal to 0, since the shadow is on the ground (xy-plane on the three-dimensional coordinate system), two equations are obtained, one for the values of the x-coordinate of the shadow and one for the values of the y-coordinates (Equations 24 and 25, respectively).
$${x}_{sh}={y}_{n}-{c}_{n}\bullet {z}_{n}$$
24
$${y}_{sh}={x}_{n}-\frac{{z}_{n}}{{b}_{n}}$$
25
where bn and cn are given by Equations 26 and 27, respectively.
$${b}_{n}=\left(\frac{{z}_{s}-{z}_{n}}{{y}_{s}-{y}_{n}}\right)$$
26
$${c}_{n}=\left(\frac{{x}_{s}-{x}_{n}}{{z}_{s}-{z}_{n}}\right)$$
27
Function 5 (PVShadowPoint) provides the computations of the shadow points, according to Equations 24–27.
Function 5: Calculation of the points that form the PV shade. |
function [x_shd,y_shd]=… PVShadowPoint(PV_x,PV_y,PV_z,x_sun,y_sun,z_sun) b_n ← Eq. 26 c_n ← Eq. 27 x_shd ← Eq. 24 y_shd ← Eq. 25 end |
2.5 Binary analysis
To represent the area covered by the greenhouse and determine the portion shaded by the photovoltaics, arrays with logical values (0 and 1) have been created. In terms of the greenhouse, the value 0 (false) corresponds to an area that is not covered by a greenhouse, while the value 1 (true) refers to an area that is enclosed by the greenhouse's limits. Because of the nature of the issue, this array will always be the same for each time of year. Regarding the shade formed by the photovoltaics, the corresponding array will have the same dimensions as the corresponding one of the greenhouse. This is because the algorithm requests the shadow on the ground inside the greenhouse to the total ground covered by the greenhouse, rather than the entire ground. As before, values of 0 and 1 represent the unshaded and shaded areas, respectively.
An array of zero values (false) is initially created based on the dimensions entered for the greenhouse. More specifically, it is assumed that the rows of the array created are the values corresponding to the y-components of the greenhouse coordinates, while the columns of the array correspond to the x-components. Due to the need for integer values of the dimensions of the greenhouse (since each row-column corresponds to a specific point of the x-y plane of the Cartesian coordinate system), and since the accuracy of the entered values (length of the greenhouse, width of the greenhouse, etc.) is of the order of centimeters (up to two decimal places), all dimensions, both for the greenhouse and for the shadows, are converted from meters to centimeters. Another transformation of the actual points of the greenhouse, but also the shade, is the additional increase of each value by 1 cm since there is no zero columns or row in the created array. Finally, after the greenhouse limits have been calculated by Function 4, the array cells corresponding to values within the greenhouse limits are replaced with the value 1 (true). At the same time, it should be noted that value 1 concerns only the cultivable unit of the greenhouse (Construction unit 2 – CU 2) and not the entire greenhouse. The code implemented for the above is presented in Section 2 – Part A.
Section 2 – Part A: Binary analysis – Greenhouse representation |
binary_array_CU_2 = false(int16(grh_length*100) + 1),… int16 ((grh_width*numCU_wdwise*100) + 1)); CU_2_limits = [CU_arrays_points {1,2}(CU_bounds_2d,1),… CU_arrays_points{1,2}(CU_bounds_2d,2)]*100; CU_2_limits(:) = CU_2_limits(:) + 1; binary_array_CU_2(CU_2_limits(1,2):CU_2_limits(3,2),… CU_2_limits(1,1):CU_2_limits(3,1)) = 1; |
As for the shade formed by the photovoltaics, the same strategy is followed with some additional necessary points. The first point concerns the conversion of all the coordinates of the shadow boundaries into positive quantities. Considered a coordinate system in which the coordinates of the sun’s position take either positive or negative values, the coordinates of each shadow point can be either positive or negative, respectively. However, an array where each row-column corresponds to a point on the x-y plane requires that each value be positive and non-zero. Thus, every non-positive coordinate of the shadow is converted into 1. The shadow boundaries may become narrower, without affecting the final desired result as the boundaries of the greenhouse consist of positive values and the purpose of the algorithm is to calculate the shading within it. Finally, for the best performance of the algorithm, the limits of the shadow that present positive values and are greater than the limits of the greenhouse are equated with those of the greenhouse limits without again affecting the final result which, as mentioned above, is the estimation of the shadow inside the greenhouse and in relation to the total surface thereof. The additional necessary points of the code implemented for the above are presented in Section 2 – Part B.
Section 2 – Part B: Binary analysis – Shadow representation |
if shading_limits(t,1) < = 0,shading_limits(t,1) = 1; end if shading_limits(t,2) < = 0,shading_limits(t,2) = 1; end if shading_limits(t,2) > = int16((grh_length*100) + 1),… shading_limits(t,2) = int16((grh_length*100) + 1); end if shading_corners(t,1)>=… int16((grh_width*numCU_wdwise*100) + 1),… shading_corners(t,1)=… int16((grh_width*numCU_wdwise*100) + 1); end |
2.6 Calculation of the greenhouse shaded surface
As mentioned above, the final desired result extracted by the algorithm concerns the surface of the greenhouse's cultivable unit that is covered by the shadow of photovoltaics. By creating arrays of logical values from the previous steps, there is the possibility of knowing if every square centimeter of the greenhouse's surface is also covered by the shadow. The basic idea is that since the number of cells of the "binary_array_CU_2" and "binary_array_shade" is the same and each cell represents the same square centimeter of land, then the number of cells of the "binary_array_CU_2" that are in the same position as the cells of the "binary_array_shade" and both contain the logical value 1 (true), represents the size of the shaded area of the greenhouse in cm2. Finally, dividing by the total number of cells of the "binary_array_CU_2" containing the value 1, the percentage of the shaded area is calculated. The above is described in Section 3 of the code.
Section 3: Greenhouse-shaded surface |
CU_2_grh_surface=… length(find(binary_array_CU_2(:,:) = = 1)); shaded_surface = length(find(binary_array_CU_2(:,:) = = 1 &… binary_array_shade(:,:) = = 1))/10.^4; shaded_surface_perc=… (shaded_surface.*10.^4/CU_2_grh_surface)*100; |
*The multiplication by 104 is for the conversion of cm2 to m2.
2.7 Complete algorithm
As mentioned above, the algorithm consists of five distinct functions and three individual sections. To derive the result from the algorithm, as input values for each sub-sequent function or section, the output values of the previous function are used. By entering the desired values for the time (year, month, day, hour, minutes, seconds, time zone), for the greenhouse (orientation, number of construction units in width and length, width and length of construction unit, gutter, and ridge height), for the geo-graphical coordinates (latitude and longitude of the location of the greenhouse), for the dimensions of the photovoltaics (length and width), and finally for the position of the photovoltaics on the roof there is the possibility, based on the way the algorithm is configured, to extract the surface (m2) that is covered by the shade of the photovoltaics within the greenhouse cultivation unit with a precision of two decimal points, the per-centage of the shaded surface in relation to the total area of the unit (%) and finally two representative two- and three-dimensional graphics. The complete algorithm (PVShading) is presented below.
Algorithm – PVShading |
Data: yr, mth, d, hr, mnt, sc, timezone, dir, numCU_wdwise, numCU_lnwise, grh_width, grh_length, gut_hght, rigde_hght, l_loc, lat, PV_width, PV_length, pos_PV_CU_1, pos_PV_CU_2 |
function [shaded_surface,shaded_surface_perc]=… PVShading(data) Function 1 Function 2 Function 3 Function 4 Section 1 Function 5 Section 2 Section 3 end |