Using a shading matrix to estimate the shading factor and the irradiation in a three-dimensional model of a receiving surface in an urban environment

Available online at www.sciencedirect.com

Solar Energy 92 (2013) 15–25 www.elsevier.com/locate/solener

Using a shading matrix to estimate the shading factor and the irradiation in a three-dimensional model of a receiving surface in an urban environment Emerson G. Melo a,⇑, Marcelo P. Almeida b, Roberto Zilles b, Jose´ A.B. Grimoni a a

Universidade de Sa˜o Paulo – USP, Escola Polite´cnica, Departamento de Engenharia de Energia e Automacßa˜o Ele´trica, Av. Prof. Luciano Gualberto – Travessa 3, 158, Sa˜o Paulo – SP 05508-900, Brazil b Universidade de Sa˜o Paulo – USP, Instituto de Eletrote´cnica e Energia, Av. Prof. Luciano Gualberto, 1289, Sa˜o Paulo – SP 05508-900, Brazil Received 11 October 2012; received in revised form 19 February 2013; accepted 20 February 2013 Available online 30 March 2013 Communicated by: Associate Editor Jan Kleissl

Abstract In an urban environment, grid-connected building integrated photovoltaic (PV) systems can be subject to complex shading patterns. The study of the shadows projected by nearby buildings and other elements around a PV surface permits cutting down energy losses due to the module’s partial shading and improving the system’s performance ratio, so that the energy production costs can be lower. This paper presents a methodology that estimates the shading factor and irradiation on a three-dimensional model of a receiving surface in an urban environment. The main innovations introduced by this methodology are the building of a shading matrix composed by direct shading factor values around the whole sky dome and the analysis of the shading impacts on direct beam, isotropic diffuse, circumsolar diffuse and horizon brightening diffuse solar radiation components. The shading matrix improves the time spent on long simulation periods and permits an easy numerical integration over the sky to obtain the diffuse shading factors. Using this feature, a plug-in to the Google SketchUp three-dimensional modeling software was built to test this methodology. A series of similar results were obtained between actual measurements and estimates conducted by the plug-in. Ó 2013 Elsevier Ltd. All rights reserved. Keywords: Shading factor; Irradiation; Three-dimensional model; Shading matrix

1. Introduction By the end of 2009, the world’s cumulative installed photovoltaic (PV) capacity was approaching 23 GW. One year later it was 40 GW. In 2011, more than 69 GW were installed globally and could produce 85 TWh of electricity every year (EPIA, 2012). The world’s largest market is still ⇑ Corresponding author. Address: Rua Anhanduı´ Mirim, 492, Sa˜o Paulo, SP 04904-200, Brazil. Tel.: +55 (11) 30915312; fax: +55 (11) 30915719. E-mail addresses: [email protected] (E.G. Melo), [email protected] (M.P. Almeida), [email protected] (R. Zilles), [email protected] (J.A.B. Grimoni).

0038-092X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.solener.2013.02.015

Europe, but the greatest growth potential in the coming years is located in countries from other continents, like China, the USA, Japan and India (EPIA, 2012). Grid-connected PV systems have the advantage of producing electrical power close to the consumers or even at the point of use, so these systems can assist in peak shaving, minimizing transmission and distribution losses and increasing grid capacity, especially in great commercial centers, where the daytime peaking load profiles are fairly coincident with solar generation profiles (da Silva Jardim et al., 2008). In urban environments, the installation of PV systems is typically integrated to the building architecture through substitution of materials usually used in building facßades, on

16

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

Nomenclature aS b DGT c cS k x / qg h Ahij AS AT e0 fB fD fH F S Dt FSd FSh FS G GB GDc

solar elevation angle (°) surface tilt angle (°) irradiance losses due to shading (W/m2) surface azimuth angle (°) solar azimuth angle (°) longitude, East positive and West negative (°) hour angle (°) latitude, North positive and South negative (°) ground reflectivity (–) solar incidence angle (°) small sky dome slice area (m2) shaded area (m2) total analyzed surface area (m2) sun-earth correction factor (–) direct shading factor (–) isotropic diffuse shading factor (–) horizon diffuse shading factor (–) shading factor on a time interval, Dt (–) daily shading factor (–) hourly shading factor (–) shading factor (–) irradiance on horizontal plane (W/m2) direct beam irradiance (W/m2) circumsolar diffuse irradiance (W/m2)

rooftops, in shading devices and other elements of the external architecture by PV modules. In countries such as Germany and Japan, studies about PV system performance in urban environments have included the module’s partial shading in the main causes for the poor performance rates of some facilities (Decker and Jahn, 1997; Kurokawa and Ikki, 2001; Ueda et al., 2009). The shading caused by obstacles around the facility as well as that caused by the facility elements, could yield an annual energy loss near to 10% (Decker and Jahn, 1997; DGS, 2008; Kurokawa and Ikki, 2001). The non-uniform irradiance distribution at the PV module surface, besides decreasing the energy available to the PV conversion, is also responsible for electrical incompatibilities between the PV cells, and, in some cases, for hot spots that can permanently damage the PV module (Bishop, 1988). A shading analysis including the elements in the surroundings of PV surfaces can eliminate or minimize the partial shading condition. Some methodologies were presented to calculate the shading factor and irradiance in PV systems subjected to partial shading. The Quaschning and Hanitsch (1998) and Drif et al. (2008) methodologies make a comparison between the solar positioning curve, determined by the solar azimuth cS and solar elevation aS angles, and a curve generated by the spherical coordinates of the existing obstacles around a specific point of the PV

GDh GDi GD GR GTs GT H TsDt HTs HT I ID ITs IT M fB rM

horizon brightening diffuse irradiance (W/m2) isotropic diffuse irradiance (W/m2) diffuse irradiance (W/m2) ground-reflected irradiance (W/m2) effective irradiance on tilted plane (W/m2) global irradiance on tilted plane (W/m2) time interval effective irradiation on tilted plane (kW h/m2) daily effective irradiation on tilted plane (kW h/ m 2) daily global irradiation on tilted plane (kW h/ m 2) hourly global irradiation on horizontal plane (kW h/m2) hourly diffuse irradiation on horizontal plane (kW h/m2) hourly effective irradiation on tilted plane (kW h/m2) hourly global irradiation on tilted plane (kW/ m 2) shading matrix (–) angular resolution scanning (°)

surface. If the curves overlay at some points, the direct irradiance will be totally or partially reduced at those points, and so the irradiance losses can be estimated. To get a precise prediction through this technique this analysis must be performed on many PV surface points. The main differences between these two methodologies are the means used to calculate the spherical coordinates of the obstacles that surround the PV surface and the diffuse irradiance losses. Alternatively, photographic methods to describe the surrounding objects have been presented (Cellura et al., 2012; Orioli and Gangi, 2012). Quaschning and Hanitsch (1995) and Cascone et al. (2011) had presented methodologies to calculate the direct shading factor through the polygon projection method and the diffuse shading factor through the radiance integration along the sky dome surfaces seen through a PV surface point. The polygon projection method decreases the complexity that exists in the analysis of large areas. The Cascone et al. (2011) methodology is more complete and also more complex than the Quaschning and Hanitsch (1995) methodology. It was developed using the MATLAB software tool and can import the external environment information from DXF files. This paper proposes a methodology to calculate the shading factor and irradiation on a three-dimensional model of receiving surfaces in an urban environment under clear, average or generic sky condition. The main

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

innovations presented by this methodology in relation to the previous ones, are the building of a shading matrix generated by the calculated direct shading factor values over the whole sky dome using the polygon projection method and the analysis of the shading effect on the direct beam, isotropic diffuse, circumsolar diffuse and horizon brightening diffuse solar radiation components. The main purpose of the shading matrix construction is the establishment of a performance versus time option, that can bring more flexibility to simulation systems with different complexities, making it possible to reduce the required time to estimate the shading factor and irradiation for large time periods, such as months and years. During the construction of the shading matrix, it is also possible to calculate, quickly and easily, the diffuse shading factors that will be used to compute the diffuse irradiance losses through the Perez et al. (1987) model, in a similar manner used by Drif et al. (2008). The main purpose of this paper is to present the proposed methodology and the results of estimates obtained through simulations realized by a plug-in built with the proposed methodology and deployed in the Google SketchUp 8 three-dimensional modeling software. Thus, measured values of a shaded PV reference cell subjected to generic sky conditions are compared with the estimated values for its three-dimensional model counterpart. 2. Solar and surface positioning in a three-dimensional environment The solar position in relation to a point on the Earth’s surface can be determined by the solar azimuth angle cS and by the solar elevation angle aS. There are various methods to calculate the value of these parameters with geo-location data, date and time (Blanco-Muriel et al., 2001; Michalsky, 1988; Reda and Andreas, 2004). Although, in a three-dimensional modeling environment using a cartesian coordinate system, the solar position, as well as any other element, is better described by a vector. In Fig. 1, a cartesian coordinate system is formed by the

17

orthogonal basis B = {x, y, z} whose origin (0, 0, 0), given by point O, is related to the latitude / and the longitude k from any location. The solar position, within that system, is given by the vector s, whose components haS, bS, cSi are given by (1). In this equation the Sun-Earth distance (1 AU) is adjusted by the sun-earth correction factor e0. haS ; bS ; cS i ¼ 1AUe0 hcos aS sin cS ; cos aS cos cS ; sin aS i

ð1Þ

Objects inside a three-dimensional modeling environment software may be created by sets of faces bounded by vertices and edges. In Fig. 1, the vertices P1, P2, P3 and P4 delimit the face M that represents a PV surface. The reading of the surface vertices allows obtaining vectors, line and plane equations as well as calculating distances, areas, volumes and angles. The vector o determines the position of the surface M in relation to the point O. The solar position relative to the surface is established by vector m = s  o. The direction of the plane containing the surface M is given by the normal vector n = (P1  P4)  (P2  P4). So, the vector n projected in the plane xy gives the vector p. The surface azimuth angle c is calculated using the dot product as in (2). In the same way, the surface tilt angle b is obtained with (3) and the solar incidence angle h calculated through (4). Finally, with (5) it is possible to calculate the area of any polygon, as the surface M, with k number of vertices.    yp 180 c ¼ arccos ð2Þ jyjjpj p      p np 180  arccos ð3Þ b¼ 2 jnjjpj p    mn 180 ð4Þ h ¼ arccos jmjjnj p    k2 1 X  ð5Þ A ¼  ðP iþ1  P 1 Þ  ðP iþ2  P 1 Þ  2  i¼1

3. Global irradiance on tilted plane Part of the solar radiation is scattered upon entering the Earth’s atmosphere. Thus, the irradiance on horizontal plane G is made up by the summation of the direct beam irradiance GB and diffuse irradiance GD. Considering the sky as anisotropic, the diffuse irradiance can also be divided in isotropic diffuse GDi, circumsolar diffuse GDc and horizon brightening diffuse irradiance GDh. In a tilted plane, there is also the ground-reflected component GR, referring to the solar radiation that is reflected by the ground into the tilted plane. Thereby, the global irradiance on tilted plane GT is composed by the direct beam, diffuse and ground-reflected component’s summation presented in (6). Fig. 1. Solar and surface’s positioning description in a three-dimensional environment.

GT ¼ GB þ GDi þ GDc þ GDh þ GR

ð6Þ

18

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

The irradiance integration over a time interval Dt results in values of irradiation measured in units of energy per square meter. The difference between the models used for estimating the irradiation on tilted plane (Gueymard, 1987; Liu and Jordan, 1960; Perez et al., 1987; Reindl et al., 1990) is basically the way in which the diffuse radiation is handled. In the proposed methodology, the Perez et al. (1987) model is used to estimate the hourly global irradiation on tilted plane IT. Despite the relative complexity, it can be easily implemented, produces good estimates on surfaces of any orientation (Duffie and Beckman, 2006; Utrillas and Martinez-Lozano, 1994) and enables the individual assessment of the shading influence on different components of the diffuse irradiance. The value of IT, in kW h/ m2, is numerically equivalent to the average irradiance for that hour in kW/m2. Thus, by adjusting the hour angle x accordingly it can be possible to estimate the GT values taking them from the IT values, that, in turn, could be estimated by the Perez et al. (1987) model. The hourly global irradiation on horizontal plane I as well as the hourly diffuse irradiation on horizontal plane ID data used as inputs in Perez et al. (1987) model, can be directly inserted in it through a comma separated text file. Thereby, information retrieved through clear sky models (Gueymard, 2008; Ineichen and Perez, 2002; Mueller et al., 2004), such as generic sky conditions data can be used. When just monthly global average irradiation data are available, synthetic hourly irradiation values are generated using the Aguiar et al. (1988) and Aguiar and Collares-Pereira (1992) models, and then, the diffuse irradiation components are derived from the Erbs et al. (1982) correlation. Thus, in surfaces where the solar radiation impinges homogeneously, the energy E available in a given area A during a time interval Dt can be calculated by means of (7). Z E ¼ GT ðtÞAdt ð7Þ

Fig. 2. GB, GDc, part of GDi and GDh in the surface M obstructed by the obstacle O.

To calculate the GTs values directly using the Perez et al. (1987) irradiance model and some geometrical parameters, some simplifications were assumed in the way that the different solar radiation components are attenuated by obstacles. Fig. 2 presents a surface M facing North. In this instance, the obstacle O will block GB, but the GDc, GDi and GDh blocked ratio will vary in accordance with sky conditions. So, considering O fully opaque, the occurrence of total blocking of GB and GDc was assumed as in Drif et al. (2008), and a geometrical parameter named direct shading factor fB will be equal to 100% in this situation. Another geometrical parameter is the isotropic diffuse shading factor fD, that is the ratio between the area C and the hemisphere H portion in front of the surface M plane, that in turn, will be used to approximate the GDi blocked amount. Finally, the GDh attenuation is related to a horizon diffuse shading factor fH, that is similar to the fD concept, but relates to the obstructed and unobstructed surfaces within a thin solid angle in the direction of the horizon.

4. Shading factor and irradiation Considering the shadow’s projection effect, the energy E available in a given area A during a time interval Dt will suffer a reduction and shall be expressed as in the Eq. (8). The shading factor FS represents the ratio between the irradiance losses due to shading projected by elements around the PV surface DGT and GT, as shown in (9). Thus, the effective irradiance on tilted plane GTs that impinges in a shaded surface is found with the Eq. (10), in other words, by subtracting DGT from the GT estimated value. Z E ¼ ½1  F S ðtÞGT ðtÞAdt ð8Þ DGT GT GTs ¼ GT  DGT ¼ GT  F S GT ¼ ð1  F S ÞGT FS ¼

ð9Þ ð10Þ

Fig. 3. Part of GDi and GDh in the surface M obstructed by the obstacle O.

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

19

In the observed condition in Fig. 3, the surface M is totally free of projected shadows, however, a portion of GDi and GDh will remain blocked. Thus, fB can assume different values throughout the day and year, but fD and fH are considered constant at any time. In (11), GT is expanded into its components GB, GDi, GDc, GDh and GR. All these components are calculated through the Perez et al. (1987) model and the shading factors fB, fD and fH will introduce the respective irradiance losses to obtain the GTs value. The GR losses can be included in the analysis by adjusting the ground reflectivity qg.

vertices. However, the Weiler and Atherton (1977) method can determine the intersection between the areas M and S and also allows the subtraction between polygon areas that overlay. Thus, the base provided by the Weiler and Atherton (1977) method is more efficient in this situation, once it allows the finding of the area AS and also performs the subtraction of the overlapped areas, projected by the multiple faces that could compose these objects.

GTs ¼ ð1  fB ÞGB þ ð1  fD ÞGDi þ ð1  fB ÞGDc þ ð1

The calculation method of fB using polygons projection has the advantage of ensuring the same performance for different surface sizes. However, when it is necessary to estimate the shading factor on a monthly or annual basis, a large amount of processing is still required. To optimize simulations for long time periods, the solar positions in the sky dome defined by cS i and aS j throughout the whole sky dome are scanned with angular resolution rM and then the fB values are previously calculated for these different solar positions. In this way, a shading matrix M fB can be built using fBij values calculated element by element for these positions. The M fB resolution rM will be responsible for determining the amount of elements contained in it, thus, the i and j indexes are incremented in the range 0 6 i 6 cS max and 0 6 j 6 aS max . The bounds of the matrix are cS max ¼ ð360  rM Þ=rM and aS max ¼ ð90  rM Þ=rM as shown in (13). 3 2 fBcS ;aS    fBcS ;aSmax fBcS ;aS 0 0 0 1 0 7 6 6 fBc ;a fBcS ;aS    fBcS ;aSmax 7 S1 S0 7 6 1 1 1 7 ð13Þ M fB ¼ 6 7 6 .. .. .. .. 7 6 . . . . 5 4 fBcSmax ;aS fBcSmax ;aS    fBcSmax ;aSmax

 fH ÞGDh þ GR

ð11Þ

4.1. Direct shading factor In Fig. 4, the vertex PS is the result of the vertex P0 projection in the plane l, using the vector m direction. This procedure is easily achieved through the plane l equation (ax + by + cz + d = 0), the point P0(x0, y0, z0) and the vector m ham, bm, cmi line equation components. So, the projection of the face vertices of obstacle O will delimit the polygon S. The intersection between the polygons S and M results in the shaded area AS. The fB value is both defined by Quaschning and Hanitsch (1995) and Cascone et al. (2011), as the ratio between AS and the total analyzed surface area AT. fB ¼

AS AT

ð12Þ

The points of intersection between the polygons M and S can be obtained by some methods, often used in computer graphics. Using the Sutherland and Hodgman (1974) method, it would be possible to determine the AS

Fig. 4. Projection of the obstacle O in the plan l and the shaded area AS defined through the M and S polygons intersection.

4.2. The shading matrix

0

1

As smaller the rM value, as greater the details included in the M fB and the simulation accuracy, nevertheless, the time spent for construct the matrix will also be greater. The set up of rM value should be carried out taking into account the required performance and the available time. To avoid unnecessary operations and speed up the M fB construction, some rules are observed in the fBij calculations. Firstly, the obstacle elevation angles in reference to all analyzed surface vertices are calculated, then, the maximum obstacle elevation angle aOmax is determined. Secondly, the solar positions where h is greater than 90° are ignored. So, the fBij values are calculated in accordance to the rules summarized in (14). ( AS ; if h 6 90 and aS j 6 aOmax fBij ¼ AT ð14Þ 0; otherwise To a given solar positioning, defined by cS and aS, the value of fB is calculated from the M fB values by a bi-linear interpolation, given by (15), so that the accuracy of estimates performed with low-resolution scanning (high rM values) is increased.

20

fðcS ; aS Þ ¼

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

fBcS ;aS i

j

ðcS i  cS iþ1 ÞðaS j  aS jþ1 Þ  aS jþ1 Þ þ

fBcS

ðcS  cS iþ1 ÞðaS

;a iþ1 S j

ðcS iþ1  cS i ÞðaS j  aS jþ1 Þ

ðcS

 cS i ÞðaS  aS jþ1 Þ þ

fBcS ;aS i

jþ1

ðcS i  cS iþ1 ÞðaS jþ1  aS j Þ

 aS j Þ þ

fBcS

ðcS  cS iþ1 ÞðaS

;a iþ1 S jþ1

ðcS iþ1  cS i ÞðaS jþ1  aS j Þ

 ðaS  aS j Þ

ðcS  cS i Þ ð15Þ

To retrieve the fBij values used in bi-linear interpolation (15) from M fB , the matrix can be searched in a straightforward manner until a given position is reached. Or, when it is necessary to know a fB value to any instant in a day, the sun position angles cS and aS corresponding to that date and time are calculated using the Michalsky (1988) method, thus, the key pair ðcS i ; aS j Þ used to address the M fB array position is computed as in (16). 8   > < cS i ¼ INT rcS M   ð16Þ > a : aS j ¼ INT r S M Nevertheless the sun does not visit all the sky dome positions contained in the shading matrix throughout the day and year, the full sky dome scanning process permits a simple numerical approximation to get the diffuse factors by using fBij values. This will be made clear in the next section. 4.3. Isotropic diffuse shading factor In Fig. 5 the hemisphere H is sliced into small surfaces hij bounded at the vertices cS i , cS iþrM , aS j and aS jþrM , for which the central point is ðcS iþr =2 ; aS jþrM =2 Þ. The area Ahij M

of a single slice hij could be approximated by a trapezoid area, with bases ðprM =180ÞcosaS j and ðprM =180ÞcosaS jþrM . This slice could be responsible for the projection of shadows on the surface M in that sun position. Retrieving from M fB a fB value that is related to the central point of hij, when this fB value is null, there will be no GDi blocking in this sky dome slice. But, when this fB value is 1, the totality of Ahij area will contribute to the GDi blocking. Then, the fB indicates where and how much GDi will be blocked in a sky dome slice. This relationship is indicated in (17). In cases where h is greater than 90°, the sun will be behind the surface M and that related area will not contribute to mitigating GDi. Even so, the fB after bi-linear interpolation is not null.   8 rM p < 180 cos aS j þcos aS jþrM ; if h 6 90 Ahij ¼ fBij ð17Þ 2 :  0; if h > 90 During the construction of shading matrix M fB , it is possible to integrate the Ahij values along the whole sky dome. Thus, the fD value related to the surface M is defined by the Eq. (18). PcSmax PaSmax i¼0 j¼0 Ahij ð18Þ fD ¼ ð1 þ cos bÞp In large areas the sky dome seen from different surface points could vary. However, once the fB values will vary along the ascending elevation angles, the fD will result in an average value that accounts for these differences in GDi attenuation through the entire surface. 4.4. Horizon diffuse shading factor The fH value is calculated in the same way as the fD value. However, all the horizon brightening irradiance is considered as impinging from a thin horizon sky band. So, will be assumed that just the Ahij areas inside a thin region defined by 0° 6 aS 6 rM will contribute to GDh attenuation. 4.5. Irradiation and shading factor at different time intervals With fB, fD and fH values calculated, the GTs value can be estimated, for any date and time, by the Eq. (11). The FS value, that accounts for the irradiance losses at that instant, can be calculated with (9). At any time interval, the effective irradiation on tilted plane H TsDt is the result of the GTs integration over a time interval bounded by ti and tf, as shown in (19). This time interval is usually hourly, daily, monthly or yearly. Z tf H TsDt ¼ GTs ðtÞdt ð19Þ ti

Fig. 5. Hypothetical hemisphere H divided into small areas hij defined by cS i , cS iþrM , aS j and aS jþrM .

To get the shading factor F S Dt referencing the losses during a time interval Dt, the solution presented in (20) can be used.

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

21

Fig. 6. Configuration A: real reference cell and its three-dimensional model counterpart.

Fig. 7. Configuration B: real reference cell and its three-dimensional model counterpart.

R tf F SDt ¼ 1  Rtitf ti

GTs ðtÞdt GT ðtÞdt

ð20Þ

5. Results and discussions To test the proposed methodology, a plug-in, hereafter referenced by Solar3DBR, was developed and deployed in Google SketchUp 8 three-dimensional modeling software. The existing classes and methods in the API Ruby of the Google SketchUp was used to perform the reading of the three-dimensional model environment geometry, as well as to read information such as latitude, longitude, date and time. So, this information and the sky condition data, when introduced in the Solar3DBR, allows the estimation of the effective irradiation and the shading factor for a selected surface in a three-dimensional model. To verify the performance of the Solar3DBR in a real situation, an experiment was conducted in the Laborato´rio de Sistemas Fotovoltaicos da Universidade de Sa˜o Paulo (LSF USP), located at 23.56° South latitude and 46.73° West longitude. In this experiment, hourly integrated irradiance data was collected with a 10  10 cm reference cell, and hourly irradiation data was measured using a pyranometer. Both the reference cell and the pyranometer were tilted by 23° and North-oriented. Using this data, the reference cell hourly shading factor FSh was determined during two distinct situations. Situation

A, shown in Fig. 6, aims to create conditions to analyze the shading effect at the beginning and at the end of the day, which is the must common shading situation, so, the shading was intensified with the use of near side shields. In situation B, shown in Fig. 7, the shading was intensified in the middle of the day with the use of an upper shield, this condition is more extreme and unusual, but permits the verification of the prediction accuracy of the diffuse irradiance losses. Besides the shields, other natural obstacles are also included in the tree-dimensional model, Fig. 8 shows the entire site model where it is possible to identify different obstacle types, such as large, thin, near and distant surfaces. Measurements performed in situation A were made on 4/5/2012, a very clear day. In situation B, the measurements were performed on 4/11/2012. In this case, the sky conditions ranged from very clear in the morning to heavily overcast in the afternoon. In each figure the real model is presented, as well as its three-dimensional model counterpart built in Google SketchUp 8. To perform simulations using the Solar3DBR, these general sky conditions was characterized by I and ID data, collected by two others pyranometers, inserted into a comma separated text file that could be readied by Solar3DBR. In Figs. 9 and 10, the FSh values determined by measurements are compared to the FSh values obtained in the simulations performed with the Solar3DBR. To perform these simulations, the shading arrays M fB were constructed with

22

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

Fig. 8. Three-dimensional model view of the site. The arrow points to the reference cell. Besides the shields, other different surface types were also included in the simulation.

1.0

1.0

shading factor

irradiation (kWh/m2)

Reference Cell 0.8

Solar3DBR

0.6 0.4 0.2 0.0 5

6

7

8

9

0.8 0.6 0.4 0.2 0.0

10 11 12 13 14 15 16 17 18 19

5

hour

6

7

8

9

10 11 12 13 14 15 16 17 18 19

hour

Fig. 9. Configuration A: measured and simulated FSh values performed on 04/05/2012.

IT Pyranometer

ITs Reference Cell

IT Solar3DBR

ITs Solar3DBR

Fig. 11. Configuration A: measured and simulated hourly irradiation values on 04/05/2012.

1.0

Reference Cell Solar3DBR

1.0

0.6 0.4 0.2 0.0 5

6

7

8

9

10 11 12 13 14 15 16 17 18 19

hour Fig. 10. Configuration B: measured and simulated FSh values performed on 04/11/2012.

a 1° scanning resolution. The IT values and hourly effective irradiation on tilted plane ITs are also compared in Fig. 11 and in Fig. 12 respectively. To compare the measured and estimated FSh values in each case, the Mean Bias Errors (MBE) and the Root Mean Square Errors (RMSE) were calculated using the data presented in Fig. 9 and in Fig. 10. Thus, 0.046 and 0.014 MBE values and 0.122 and 0.099 RMSE values, were obtained in situation A and B respectively.

irradiation (kWh/m2)

shading factor

0.8

0.8 0.6 0.4 0.2 0.0 5

6

7

8

9

10 11 12 13 14 15 16 17 18 19

hour IT Pyranometer

ITs Reference Cell

IT Solar3DBR

ITs Solar3DBR

Fig. 12. Configuration B: measured and simulated hourly irradiation values on 04/11/2012.

The measured FSh values rely on uncertainties that are introduced in both pyranometer and reference cell measurements mainly due to temperature, solar incidence angle

E.G. Melo et al. / Solar Energy 92 (2013) 15–25 Table 1 Simulated and measured daily values. Configuration

A

Values

Measured 2

HT (kW h/m ) HTs (kW h/m2) FSd

Table 2 FSd, HTs and the time spent in simulations realized with different rM values.

B

6.318 5.089 0.194

Predicted 6.400 5.104 0.202

Measured 5.371 2.290 0.574

Predicted

rM

FSd

HTs (kW h/m2)

Time (s)

5.609 2.581 0.539

1° 2° 5° 10°

0.203 0.203 0.203 0.204

5.104 5.102 5.099 5.095

955 250 39 10

and spectral response effects (Dunn et al., 2012; Gueymard and Myers, 2008). Besides the existence of some inaccuracies in the three-dimensional models, theses factors could have contributed to the points that were found to be the greater divergences, as in the beginning and in the end of these days, where the incidence angle is higher and the irradiance level is lower, so, a minor difference in measurements can insert major errors. The temperature effect is well observed by comparing case A and case B FSh midday hours. At this time, the FSh values in case A, subject to direct radiation, are higher than in case B, subject only to diffuse radiation that results in a lower temperature. The curves shown in Figs. 11 and 12 demonstrate a good accuracy in the IT calculation. The MBE are 0.006 and 0.018 kW h/m2 and the RMSE are 0.034 and 0.035 kW h/ m2 in case A and B respectively. So, the basic fundamental for estimating the effective irradiance was correctly established. Regarding the shading, in case A, the ITs curve presents 0.001 kW h/m2 MBE value and 0.039 kW h/m2 RMSE value. Meanwhile, the MBE value in case B is 0.022 kW h/m2 and the RMSE value is 0.055 kW h/m2, which is just a little greater than the previous values. Despite the FSh values on 4/11/2012 midday hours showing great similarity, the greater MBE and RMSE in case B could be explained by the simplifications introduced in the calculation of the diffuse irradiation losses, since, in the greater part of this day, the upper shield blocks GB, so that, the diffuse component dominates. 1.0 0.8

shading factor

23

Table 1 summarizes the daily global irradiation on tilted plane HT, and the daily effective irradiation on tilted plane HTs, that were obtained by IT and ITs integration throughout each day, and the daily shading factor FSd respective to both configurations. If just the Perez et al. (1987) model had been used, without accounting for the shading irradiance losses, the predicted value at 4/5/2012 (configuration A) could be 25.8% higher. On day 11/5/2012 (configuration B) the difference could be 44.9%. Using the Solar3DBR, the differences between measured and estimated daily values are 0.3% and 12.7% respectively. The accuracy and the required time spent to build the shading matrix M fB were analyzed running configuration A simulations with rM adjusted to 1°, 2°, 5° and 10°. In the Fig. 13, the respective FSh values are presented. Through Fig. 13, it is observed that the estimates of FSh are practically not modified for small rM values as 1°, 2° and 5°. With rM equal to 10°, the values are slightly different from the others at the maximum and minimum points. Table 2 lists the values of FSd, HTs and the time spent to complete the construction of the shading matrix. There is almost no difference in the values of FSd and HTs, which proves the efficiency of the bi-linear interpolation method. The great time difference between the simulations performed with rM adjusted in 1° and the ones whose rM value were 10° is justified by the M fB size, that is 32.400 and 324 respectively. However, the M fB construction with high resolution ensures that, in three-dimensional models with thinner and more distant details from the analyzed surface, the influence of any site element will not be omitted.

0.6

6. Conclusion

0.4

A methodology was proposed to estimate the shading factor and the irradiation on surfaces of three-dimensional models of receiving surfaces in an urban environment. To verify this methodology, a plug-in deployed to Google SketchUp 8 was built and called Solar3DBR. The accurateness of the simulations performed using Solar3DBR were assessed through comparison with actual values obtained in an experiment conducted in the LSF USP. The Solar3DBR toolbar deployed in Google SketchUp 8 and the shading matrix selection box are showed in Fig. 14. Even in a reduced dimension reference cell, shaded through near, distant, thin and larger obstacles, a 0.3%

0.2 0.0 6

7

8

9

10

11

12

13

14

15

16

17

18

hour FSh Reference Cell

FSh 5°

FSh 1°

FSh 10°

FSh 2° Fig. 13. Comparison of simulated FSh values performed with different rM values.

24

E.G. Melo et al. / Solar Energy 92 (2013) 15–25

Fig. 14. Solar3DBR toolbar deployed in Google SketchUp 8 and the shading matrix selection box.

difference between the measured and simulated daily values was obtained in case A. It is better than the 2.15% difference reported by Drif et al. (2008). Assuming that in case B, the upper shield recreates a worse and unusual condition to PV systems, the 12.7% difference in the estimated value is not probably to occur in a real PV application. So, the simplifications introduced to calculate the diffuse irradiance losses thought the Perez et al. (1987) model result in a good precision and simplicity relationship, and permit faster calculations than other methods based on radiance models such as Brunger and Hooper (1993), Perez et al. (1993). Simulations performed with the shading matrix built with different resolutions, showed that the bi-linear interpolation method used to get the fB values from the shading matrix is efficient. So it is possible to perform large time period simulations in reduced time intervals, without significantly affecting the accuracy of the results. This is important because it creates possibilities for using this methodology in a web-based application. Thereby, the proposed methodology has been confirmed as a promising way to estimate the shading factor and irradiation in a three-dimensional model surface in an urban environment. Further long-time measurements will be necessary to verify the accuracy of the methodology under seasonal climatic variations. Comparisons between measured and simulated data from other different surface sizes and orientations have been also planned. This methodology could be applied to other solar energy uses, such as solar thermal power, and it could also be expanded in the future to insert the PV module’s optical losses and the electrical and thermal PV models. Thereby, the Solar3DBR could be used for the development of complete PV systems that will allow PV system designers to research the best surfaces as well as the better positioning for the PV modules installation. Another useful feature is the use of all facilities provided in the Google SketchUp modeling environment in the threedimensional modeling process. Such features are the simple and intuitive modeling tools, visual shading analysis, geolocation, building modeling, image import, among others. This can reduce the time spent in modeling process and

increase the three-dimensional model details, resulting in better simulation results References Aguiar, R., Collares-Pereira, M., 1992. Tag: a time-dependent, autoregressive, gaussian model for generating synthetic hourly radiation. Solar Energy 49, 167–174. Aguiar, R., Collares-Pereira, M., Conde, J., 1988. Simple procedure for generating sequences of daily radiation values using a library of markov transition matrices. Solar Energy 40, 269–279. Bishop, J., 1988. Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits. Solar Cells 25, 73–89. Blanco-Muriel, M., Alarcn-Padilla, D.C., Lpez-Moratalla, T., Lara-Coira, M., 2001. Computing the solar vector. Solar Energy 70, 431–441. Brunger, A.P., Hooper, F.C., 1993. Anisotropic sky radiance model based on narrow field of view measurements of shortwave radiance. Solar Energy 51, 53–64. Cascone, Y., Corrado, V., Serra, V., 2011. Calculation procedure of the shading factor under complex boundary conditions. Solar Energy 85, 2524–2539. Cellura, M., Gangi, A.D., Orioli, A., 2012. A photographic method to estimate the shading effect of obstructions. Solar Energy 86, 886–902. da Silva Jardim, C., Rther, R., Salamoni, I.T., de Souza Viana, T., Rebechi, S.H., Knob, P.J., 2008. The strategic siting and the roofing area requirements of building-integrated photovoltaic solar energy generators in urban areas in brazil. Energy and Buildings 40, 365–370. Decker, B., Jahn, U., 1997. Performance of 170 grid connected pv plants in northern germanyanalysis of yields and optimization potentials. Solar Energy 59, 127–133, Selected Proceeding of ISES 1995: Solar World Congress. Part IV. DGS, 2008. Planning and Installing Photovoltaic Systems: A Guide for Installers, Architects and Engineers, second edition. Earthscan, Berlin. Drif, M., Prez, P., Aguilera, J., Aguilar, J., 2008. A new estimation method of irradiance on a partially shaded pv generator in gridconnected photovoltaic systems. Renewable Energy 33, 2048–2056. Duffie, J., Beckman, W., 2006. Solar Engineering of Thermal Processes, third ed. John Wiley & Sons, New Jersey. Dunn, L., Gostein, M., Emery, K., 2012. Comparison of pyranometers vs. pv reference cells for evaluation of pv array performance. Proceedings of the 38th IEEE Photovoltaic Specialists Conference (PVSC). IEEE, Austin. EPIA, 2012. Global Market Outlook for Photovoltaics Until 2016. Technical report, European Photovoltaic Industry Association, Brussels. Erbs, D., Klein, S., Duffie, J., 1982. Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation. Solar Energy 28, 293–302.

E.G. Melo et al. / Solar Energy 92 (2013) 15–25 Gueymard, C., 1987. An anisotropic solar irradiance model for tilted surfaces and its comparison with selected engineering algorithms. Solar Energy 38, 367–386. Gueymard, C.A., 2008. Rest2: High-performance solar radiation model for cloudless-sky irradiance, illuminance, and photosynthetically active radiation validation with a benchmark dataset. Solar Energy 82, 272– 285. Gueymard, C.A., Myers, D.R., 2008. Modeling Solar Radiation at the Earth’s Surface: Recent Advances. Springer (Chapter 1). Ineichen, P., Perez, R., 2002. A new airmass independent formulation for the linke turbidity coefficient. Solar Energy 73, 151–157. Kurokawa, K., Ikki, O., 2001. The japanese experiences with national pv system programmes. Solar Energy 70, 457–466. Liu, B.Y., Jordan, R.C., 1960. The interrelationship and characteristic distribution of direct, diffuse and total solar radiation. Solar Energy 4, 1–19. Michalsky, J.J., 1988. The astronomical almanac’s algorithm for approximate solar position (19502050). Solar Energy 40, 227–235. Mueller, R., Dagestad, K., Ineichen, P., Schroedter-Homscheidt, M., Cros, S., Dumortier, D., Kuhlemann, R., Olseth, J., Piernavieja, G., Reise, C., Wald, L., Heinemann, D., 2004. Rethinking satellite-based solar irradiance modelling: the solis clear-sky module. Remote Sensing of Environment 91, 160–174. Orioli, A., Gangi, A.D., 2012. An improved photographic method to estimate the shading effect of obstructions. Solar Energy 86, 3470– 3488.

25

Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D., 1987. A new simplified version of the perez diffuse irradiance model for tilted surfaces. Solar Energy 39, 221–231. Perez, R., Seals, R., Michalsky, J., 1993. All-weather model for sky luminance distributionpreliminary configuration and validation. Solar Energy 50, 235–245. Quaschning, V., Hanitsch, R., 1995. Shade calculations in photovoltaic systems, in: ISES Solar World Conference, ISES, Zimbabwe. Quaschning, V., Hanitsch, R., 1998. Irradiance calculation on shaded surfaces. Solar Energy 62, 369–375. Reda, I., Andreas, A., 2004. Solar position algorithm for solar radiation applications. Solar Energy 76, 577–589. Reindl, D., Beckman, W., Duffie, J., 1990. Evaluation of hourly tilted surface radiation models. Solar Energy 45, 9–17. Sutherland, I.E., Hodgman, G.W., 1974. Reentrant polygon clipping. Communications of the ACM 17, 32–42. Ueda, Y., Kurokawa, K., Kitamura, K., Yokota, M., Akanuma, K., Sugihara, H., 2009. Performance analysis of various system configurations on grid-connected residential pv systems. Solar Energy Materials and Solar Cells 93, 945–949, 17th International Photovoltaic Science and Engineering Conference. Utrillas, M., Martinez-Lozano, J., 1994. Performance evaluation of several versions of the perez tilted diffuse irradiance model. Solar Energy 53, 155–162. Weiler, K., Atherton, P., 1977. Hidden surface removal using polygon area sorting. SIGGRAPH Computer Graphics 11, 214–222.