MOKSLAS – LIETUVOS ATEITIS SCIENCE – FUTURE OF LITHUANIA

The model of solar radiation, which takes into account direct, diffused, and reflected components of solar energy, has been presented. Model is associated with geographical coordinates and local time of every day of the year. It is shown that using analytic equations for modelling the direct component, it is possible to adopt it for embedded systems with low computational power and use in solar tracking applications. Reflected and diffused components are especially useful in determining the performance of photovoltaic modules in certain location and surroundings. The statistical method for cloud layer simulation based on local meteorological data is offered. The presented method can’t be used for prediction of weather conditions but it provides patterns of solar radiation in time comparable to those measured with pyranometer. Cloud layer simulation together with total solar radiation model is a useful tool for development and analysis of maximum power point tracking controllers for PV modules. keywords: solar radiation, solar energy sources, cloud layer simulation, photovoltaic (PV) modules, maximum power point tracking (MPPT).


Introduction
Constantly increasing prices for fossil fuels, ecological problems and international regulations (Directive 2009/28/ EC) draws significant attention to development of new and efficient ways to utilise renewable energy such as wind and solar power.
While new technologies for gathering renewable energy are introduced constantly, price and efficiency of sources still remain a big issue. Solar energy deserves special attention, because by far it is the largest energy resource (Quaschning 2016), and it is easily available for domestic users for electricity and heat generation. Photovoltaic modules are most commonly used to produce electrical energy from solar radiation. Due to high price per Watt and low efficiency (up to 20%) of PV modules, lot of efforts are made not only in development of modules but also in research of sophisticated control systems, ensuring minimal energy losses such as maximum power point tracking (MPPT) controllers (Kohata et al. 2010). Having mathematical model of solar irradiance allows estimation of performance of different control techniques before investing to hardware, so increasing the speed and decreasing the price of development. There are many papers available offering different approaches to terrestrial solar irradiance modelling 290 based on data which can be received from local hydro meteorological service. The aim of this model is to imitate different weather conditions, cloud motion due to wind, monthly and annual solar energy prediction for given area.

Direct solar radiation
As described in previous work (Vasarevičius, Martavičius 2011), direct solar insolation EH ( ) I n falling to the horizontal ground plane can be derived geometrically from Earth's orbit around the Sun, Earth's rotation around its own axis and heel of axis. Also solar beam attenuation in atmosphere has to be assumed: here: n -number of the day in a year (January 1 -day 1, December 31 -365), S0 1367 I = W/m 2 is a solar constant (Zekai 2008). ( ) n ε -eccentricity factor of Earth's orbit, Z θ -Zenith angle, τ -attenuation of solar energy flux in atmosphere. First three terms in expression are common to any location on Earth while the attenuationτ can be highly dependent on climatic conditions, air humidity e.t.c. in specific location. This term also estimates a height above sea level. There are several approaches proposed for calculation of attenuation in atmosphere. For current model two methods have been tested.
Simple equation presented in (Grace 2006) for estimation of attenuation in atmosphere is: here: kH is optical depth of atmosphere and can be derived from solar radiation data using Labert-Beer law. A better performance of the model achieved using empirical equation proposed in (Hottel 1976): here: 1 1 a r a = , and * k k r k = -are the factors characterizing the attenuation of solar radiation in the atmosphere, h -height above sea level in kilometers. Correction coefficients 0 r , 1 r and k are offered for four climatic zones in (Hottel 1976), but usually need to be calculated according local radiation data using Levenberg-Marquardt algorithm (LMA) (Vieira et al. 2002).
The Zenith angle Z θ describes the height of the sun above the horizon. Zenith angle equals the angle formed by the sun's rays falling to land surface and the direction of zenith. This angle defines Earth's rotation around its axis and geographic location using the following expression: here: β -altitude (elevation) angle, LA φ -latitude angle of location A, deg., M MN MA ∆φ = φ − φ -longitude difference angle between local noon MN φ and location A MA φ , rad, δ -distance from the sun to equator (declination) angle, rad (Masters 2004): Declination angle δ defines the heel of the Earth axis which is a key factor in determining the change of seasons.
To make model local time dependent, longitude difference angle between local noon meridian and the location A M ∆φ is replaced with the following expression: here: CT t -local clock time, ð / 12 w = -an angular velocity of the Earth's rotation around its axis, rad/h, -solar noon time correction due to Earth's elliptical orbit around the sun on n-th day, h (Vasarevičius, Martavičius 2011). Angle M ∆φ is sometimes referred as hour angle and can be used to calculate azimuth angle S CT ( ) t α : Described direct solar radiation model was implemented as Matlab/Simulink block. To test calculations of Zenith and azimuth angles polar solar path diagram was generated ( Fig. 1). Radial axis represents the altitude angle β . When 90 β =°, the sun is at zenith and when 0 β =°, the sun is right on the horizon line. The angular axis of this plot shows azimuth angle S α . This plot allows determining the position of the sun at the sky at any time of any day of the year.
Also Equations (7) -(11) were used to create microcontroller program for sun tracking. Results of real life testing during one month (July) have shown mean error 1.3° in Azimuth angle calculations and 1.6° in Zenith angle calculation. Although these errors can't be assumed to be decisive, because there was no suitable hardware available to properly calibrate suns direction sensor.

Diffused and reflected components
To estimate total amount of solar energy falling on to a PV modules plane, diffused and reflected solar irradiance components must be assumed (Grace 2006).
Diffused radiation component D CT ( , ) I n t appears due to scattering of light in Earth's atmosphere. Precise calculation of this component is rather difficult due to complex processes, light reflections and scattering in atmosphere. To simplify calculations the sky is assumed to be isotropicthe light reaches horizontal plane with equal intensity from all directions (Masters 2004). This assumption allows usage of following expression: here: D τ -clear sky diffusion factor, V α -slope of PV modules active plane, E CT ( , ) I n t -direct solar insolation, at time CT t . In (Liu, Jordan 1960) empirical relationship is derived between clear sky diffusion factor D τ and atmosphere attenuation τ : The last component of solar energy R CT ( , ) I n t is reflected from the surroundings in front of the panel (Masters 2004). Its power is mostly dependent on reflectance of the surface, commonly called albedo -ρ. In some cases reflected component can be deterministic factor in choosing the location for PV panel. For example, albedo of the snow can be up to 0.85. Albedo of water represents special case: pond of still water in noon can have albedo as low as 0.02, but in morning and evening, when β is small, its value can reach 0.8. If area is windy and there is a high probability of waves on surface of pond, its albedo can be 0.5 which is anyway higher than grass (~0.2) or forest (~0.15) (Zekai 2008). Assuming that Solar radiation, falling on to surroundings of PV module is reflected equally in all directions, the ref- Having all three components, total clear sky solar irradiance C CT ( , ) I n t on inclined plane can be calculated: I n t I n t I n t I n t here: θ -is the angle between normal to the surface of PV module and incident solar beam. This angle is calculated according (Dusabe et al. 2009): here: H α -angle between horizontal projection of normal to the surface of PV module and a direction of south. It is necessary to check Equation (16) for a negative result and set it to zero in such a case.
Equation (15) provides total radiation to the surface of PV module during a cloudless day. Unfortunately according statistical solar radiation data, only one or two such days can happen during a month. Fast moving small clouds are very common to this region. So to have a model imitating local conditions it is mandatory to assume the influence of the cloud layer.

Cloud layer simulation
There are many methods proposed for cloud layer modeling in papers (Iziomon, Mayer 2002;Taylor, Ellingson 2008;Davies 1989), but most of them are developed for meteorological purposes and usually estimate total average values for some period of time. These models are also derived from long time statistical data in the region of interest. The goal of this work is to develop mathematical model providing following functions: − simulate different environmental condition in a real time with desired time step, − estimated total solar energy output during the month and annually should not significantly differ from local average data. For the first function presented models are not suitable. Also there is no detailed statistical data for the region of interest. These limitations require different approach to model development.
Only available data from local hydro meteorological service is monthly averages of cloud cover. The amount of clouds in the sky is estimated in 10 point scale: 0very clear sky, 10 -all the sky is covered with clouds. Annual average for Lithuania is 7.3 points. The most sunny days are observed in June and July (average ~4 points), and October and November are the most cloudy months (~8 points).
In work of Soubdhan et al. (2009) a method to classify daily clearness index t k distributions using mixture Dirichlet distributions. Solar radiation data sampled with frequency of 1 Hz was used for this work. In (Chia, Hutchinson 1991) beta distribution is offered as a model for frequency distribution of daily cloud duration.
While there is only limited solar radiation data available for local region, another type of asymmetric probability density function (PDF) -gamma distribution is chosen empirically for present work to define the probability of sky coverage by clouds for a given day. In a similar manner as in (Soubdhan et al. 2009) and (Chia, Hutchinson 1991) the peak of PDF is chosen in known average for a given month.
Two modifications are performed on gamma distribution: it is mirrored when average percentage of sky covered by clouds for a given month is higher than 50% and expanded beyond limits of 1-10 to leave a small probability of clear day during very cloudy month and bad weather during relatively clear month. Such PDF is described by following equation: here: variable 0 10 x ≤ ≤ , i = 1,2,…,12 -the number of month, where 1 is January, i α and i β -shape coefficients derived for each month, i ψ -normalization coefficient, i p -average part of sky, covered by clouds for a selected month, i Ã( ) α -gamma function. As seen from (17), cloud cover during each month can be described using set of three members: i α , i β and i ψ . Due to widening of function it does not integrate to 1, so normalization coefficient i ψ has to be calculated for each set of members i α and i ψ : (18) Figure 2 shows probability density functions for four months: February, May, August and November. For any given day n, one random value of expected cloud coverage p is generated according these functions. To define the duration of the single cloud, random normally distributed values with the mean p are generated, which control variable pulse generator created in Matlab/Simulink environment. As a result pulses of variable duration are generated with average duty ratio equal to p .
For analysis of MPPT systems the rate of change of solar radiation near the edges of the clouds is important. Usually that rate depends on the wind speed -at higher wind clouds move faster and solar radiation falling on PV modules surface changes faster. This process can be simulated by passing the output of pulse generator through a Gaussian filter with following impulse response ( ) h t : here: B(ν) -pass-band of the filter, dependent on the wind speed ν, T S -sampling period of input signal, t -time.
To simulate low wind speed, narrow pass-band B(v) of the filter is set, giving slow rise and fall rates of output signal. For high wind speed -wide filter pass-band is set and output signal is changing rapidly.
Low power white noise with filtered high frequency is added to the resulting signal to simulate fluctuations of solar radiation due to such effects as warm air rising from the surface of PV module etc.
Resulting signal can be denoted as transparency ratio of the sky -( ) t γ . When the sky is absolutely clear -γ = 1, when direct path of the solar beam is fully covered by cloud -γ = 0. The albedo of low clouds can be as high as 0.6 (Zekai 2008), so the light reflected from the clouds also reach surface of PV module and in some cases total instantaneous radiation can become higher than of clear sky (Soubdhan et al. 2009).
There are some methods offered in Maxwell et al. (1986) and Suckling, Hay (1977) to evaluate multiple reflections from clouds. In present work a different approach is being tested -to derive an average relation between direct beam and diffused radiation.
For a cloudy day Expression (13) cannot be used to determine the diffuse factor of the sky. In this case sky is anisotropic and diffusion of radiation is much more complicated process.
In Figure 3 normalized data sets of measured diffused and direct beam radiation (hourly averages) during June in Lithuania are presented. Each series represent data pairs of one hour during 30 days. 12 hours of day time are analyzed. Early morning and late evening data is excluded and should be handled separately.
In Figure 3 it is obvious, that when the sky is clearmostly direct radiation reaches the surface of PV module and in very cloudy conditions only diffused portion remains. A line represents approximated relationship derived using LMA: ( ) Replacing sky diffusion factor D τ in Equation (12) with CT ( ( )) t ψ γ , according to (14) and (15) The presented model is implemented in Matlab/ Simulik environment to produce real time solar radiation values on a given day of the year in any given location in Lithuania. In Figure 4 simulation results of 3 solar radiation components during one day of June in Lithuania are presented: direct radiation to horizontal plane EH I for clear day as a reference and total T I and diffused D I solar radiation components for a cloudy day. As seen in figure, when a day is perfectly clear, the diffused component is comparatively weak. During cloudy day the direct component can be totally attenuated and the diffused component increases significantly due to reflections from clouds. The effect of multiple reflections described in (Suckling, Hay 1977) is achieved by splitting D I to multiple elements with different weight coefficients (integrating to one) and delaying or advancing each element by random value. Due to this some instantaneous values of T I can be higher than those obtained during a clear day. Such situation can be seen in charts presented in (Soubdhan et al. 2009).
To increase a precision of the presented model for real time cloud simulation additional parameter -Zenith angle Z θ has to be assumed in Equation (21) and solar radiation data of at least one year has to be processed to perform surface matching algorithm for derivation of relationship Summarising it can be stated, that the developed solar irradiance model can generate continuous solar radiation data patterns, comparable to those naturally occurring. Has been shown, that the developed model is suitable for development and performance analysis of control systems for solar thermal collectors, PV panels (especially MPPT tracking algorithms) and other energy sources, utilising solar energy.

Conclusions
Mathematical model of solar radiation can significantly increase the speed and reduce the price of development of control systems for solar energy sources. The paper presents such a model assuming four elements: direct, diffused and reflected solar radiation components and influence of the cloud layer on total solar radiation falling on to the inclined plane. Model of direct radiation can be used to calculate azimuth and altitude angles of the sun at any location, defined by geographical coordinates at any time. This makes it useful for solar tracking applications and optimal angle calculations of PV modules and solar thermal collectors. Described model assumes reflectivity of surrounding environment, which in some cases can have significant influence on performance of solar energy source. The method is presented for simulation of the cloud layer using gamma distribution and simple local statistical meteorological data. As shown, model generates solar radiation patterns in time comparable with real data, so making the model suitable for analysis and development of MPPT algorithms.
Further investigation of the proposed method for cloud simulation is needed to increase its accuracy and will be provided when significant statistical data will be available for a region of interest.