solarpower.cpp

Go to the documentation of this file.
00001 //          Solar energy predictor
00002 //
00003 // This program executes a model of the solar energy (kWH/m^2) arriving at a
00004 // solar module that is set at an arbitrary angle, and on an arbitrary day of
00005 // the year. It works by first proposing a Glenn Research Centre model of the
00006 // atmospheric absorption as a function of height, based on a model of air
00007 // density. From the known absorption for sunlight arriving when the sun is
00008 // vertically overhead, a model of the absorption of slant rays is developed.
00009 //
00010 // Using an electrical model of a solar module feeding a battery, it is
00011 // possible to predict the total amount of charge passed into the battery and
00012 // hence the effectiveness of the system. Because the solar module is not
00013 // necessarily working at its maximum power point (MPP), there is an additional
00014 // drop in efficiency that can be demonstrated.
00015 //
00016 // The aim is to predict the amount of energy obtained during the shortest
00017 // clear day for the purposes of understandong the performance of a solar
00018 // module in a battery charging system, and dimensioning a solar power system.
00019 //
00020 /***************************************************************************
00021  *   Copyright (C) 2007 by Ken Sarkies                                     *
00022  *   ksarkies@trinity.asn.au                                               *
00023  *                                                                         *
00024  *   This program is free software; you can redistribute it and/or modify  *
00025  *   it under the terms of the GNU General Public License as published by  *
00026  *   the Free Software Foundation; either version 2 of the License, or     *
00027  *   (at your option) any later version.                                   *
00028  *                                                                         *
00029  *   The program is distributed in the hope that it will be useful,        *
00030  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00031  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00032  *   GNU General Public License for more details.                          *
00033  *                                                                         *
00034  *   You may obtain a copy of the GNU General Public License by writing to *
00035  *   the Free Software Foundation, Inc.,                                   *
00036  *   51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA.             *
00037  ***************************************************************************/
00038 
00039 #include "solarpower.h"
00040 #include <iostream>                                 // Base stream classes
00041 #include <fstream>                                  // File I/O
00042 #include <cstring>
00043 #include <cstdlib>
00044 #include <cmath>
00045 
00046 int main(int argc, char ** argv)
00047 {
00048 // General Constants
00049     const double lossConstant = 0.253617853024622586/pathLoss(1);
00050     const double solarConstant = 1366;          // W/m^2 outer atmosphere
00051     const double solarStandard = 1000;          // W/m^2 used to specify module
00052     const double angleConversion = 3.1415927/180.0;
00053 // Model and Site specific constants
00054     const double maxDeclination = 23.45;        // Northern tropic
00055     const double latitude = -34.929;            // Adelaide
00056     const double declination = maxDeclination;  // Worst case winter
00057     const double moduleAngle = maxDeclination;  // Noon sun facing
00058     const double batteryVoltage = 12;
00059     const double modulePower = 125;
00060 // Energy to battery charge conversion factor
00061     const double energyCharge = modulePower/solarStandard/batteryVoltage;
00062 // Transformed constants
00063     const double rDeclination = declination*angleConversion;
00064     const double cosDeclination = cos(rDeclination);
00065     const double sinDeclination = sin(rDeclination);
00066     const double rModuleAngle = moduleAngle*angleConversion;
00067     const double cosModuleAngle = cos(rModuleAngle);
00068     const double sinModuleAngle = sin(rModuleAngle);
00069     const double rLatitude = latitude*angleConversion;
00070     const double cosLatitude = cos(rLatitude);
00071     const double sinLatitude = sin(rLatitude);
00072 // Model variables
00073     double minuteIncr = 1;                // time integration step size
00074     double minute = 500;                  // Length of integration (minutes)
00075     double cosAngle = 1;
00076     double solarEnergyFixed = 0;
00077     double solarEnergyFollowing = 0;
00078     double solarEnergyFixedCharge = 0;
00079     double solarEnergyFollowingCharge = 0;
00080     double charge12fixed = 0;             // Charge when solar module 12V
00081     double chargeVfixed = 0;              // Charge for solar module at MPP
00082     double charge12following = 0;         // Charge when solar module 12V
00083     double chargeVfollowing = 0;          // Charge for solar module at MPP
00084 
00085     for (double i=1; i<101; i++)
00086     {
00087         double V = OptimalModuleVoltage(i);
00088         double I = moduleCurrent(i,V);
00089         std::cout << i << "," << V << "," << I << ","
00090                   << V*I
00091                   << std::endl;
00092     }
00093     return 0;
00094 
00095 // Compute the power incident on the module during the time interval
00096     while (minute > 0)
00097     {
00098 // Longitudinal angle associated with the movement of the earth at the time,
00099 // relative to a longitudinal axis at noon.
00100         double cosHourAngle = cos(0.25*minute*angleConversion);
00101 // Angle of the sun to the vertical axis at the site and at the time
00102         cosAngle = cosLatitude*cosDeclination*cosHourAngle
00103                  + sinLatitude*sinDeclination;
00104 // Angle of the sun to the module orthogonal axis
00105         double cosIncidence = cosModuleAngle*cosDeclination*cosHourAngle
00106                             + sinModuleAngle*sinDeclination;
00107 // Solar energy received (W/m2) by a fixed module facing the sun at noon
00108         if (cosIncidence > 0)
00109             solarEnergyFixed = solarConstant*cosIncidence*
00110                            exp(-lossConstant*pathLoss(cosAngle));
00111 // Solar energy received by a following module (always facing the sun)
00112         else
00113             solarEnergyFixed = 0;
00114         solarEnergyFollowing = solarConstant*
00115                        exp(-lossConstant*pathLoss(cosAngle));
00116 // Compute the current passed into a 12V battery and a 125W module
00117 // Percentage of solar energy received relative to the standard
00118         double solarEnergyRatioFixed = solarEnergyFixed*100/solarStandard;
00119         double solarEnergyRatioFollowing =
00120                                   solarEnergyFollowing*100/solarStandard;
00121 // Voltage at the Maximum Power Point (MPP) of the module
00122         double Vfixed = OptimalModuleVoltage(solarEnergyRatioFixed);
00123         double Vfollowing = OptimalModuleVoltage(solarEnergyRatioFollowing);
00124 // Current into battery of the module is held at the battery voltage
00125         double current12fixed =
00126                     moduleCurrent(solarEnergyRatioFixed,batteryVoltage);
00127         double current12following =
00128                     moduleCurrent(solarEnergyRatioFollowing,batteryVoltage);
00129 // Current into the battery if the module is at its MPP
00130         double currentVfixed = moduleCurrent(solarEnergyRatioFixed,Vfixed);
00131         double currentVfollowing =
00132                     moduleCurrent(solarEnergyRatioFollowing,Vfollowing);
00133 // Integration of charge for efficiency results
00134 // For the MPP case, an efficient regulator drops the voltage to the battery
00135 // voltage and increases the current keeping the output power the same as the
00136 // input power. In practice such regulators are rarely more than 90% efficient.
00137         charge12fixed += current12fixed;
00138         chargeVfixed += currentVfixed*Vfixed/batteryVoltage;
00139         charge12following += current12following;
00140         chargeVfollowing += currentVfollowing*Vfollowing/batteryVoltage;
00141 // Charge received if module could work at its full power (simplistic model)
00142         solarEnergyFixedCharge += solarEnergyFixed*energyCharge;
00143         solarEnergyFollowingCharge += solarEnergyFollowing*energyCharge;
00144         minute -= minuteIncr;
00145     }
00146     std::cout << "Total (Full Power) charge - Following (AH) "
00147               <<  solarEnergyFollowingCharge/60
00148               << std::endl;
00149     std::cout << "Total charge Real 12V - Following     (AH) "
00150               << charge12following/60
00151               << std::endl;
00152     std::cout << "Total charge Efficient MPP - Following(AH) "
00153               << chargeVfollowing/60
00154               << std::endl;
00155     std::cout << "Total (Full Power) charge - Fixed     (AH) "
00156               <<  solarEnergyFixedCharge/60
00157               <<  std::endl;
00158     std::cout << "Total charge Real 12V - Fixed         (AH) "
00159               << charge12fixed/60
00160               << std::endl;
00161     std::cout << "Total charge Efficient MPP - Fixed    (AH) "
00162               << chargeVfixed/60
00163               << std::endl;
00164     return 0;
00165 }
00166 
00167 //----------------------------------------------------------------------------
00168 // Model for solar module BP3125
00169 // Match three current/voltage points to obtain the characteristic parameters
00170 // of the model pn junction. Model is correct only at these points
00171 // (short-circuit, open-circuit and maximum power).
00172 // Solar current generated is proportional to the short circuit current
00173 // portion only, the second term is the diode loss term.
00174 //
00175 // solarEnergy is the percentage standard incident solar radiation used to
00176 // define the module characteristics (ie 1000 W/m2).
00177 // Voltage is that which is forced on the module by the system.
00178 
00179 double moduleCurrent(const double solarEnergy, const double voltage)
00180 {
00181     const double Isc = 8.02;      // Short Circuit Current
00182     const double I0 = 0.000185;   // Diode dark current
00183     const double Vk = 2.071;      // Model parameter voltage
00184     double b = Isc*solarEnergy*0.01/I0+1;
00185     double current = I0*(b-exp(voltage/Vk));
00186     if (current < 0) current = 0;
00187     return current;
00188 }
00189 //----------------------------------------------------------------------------
00190 // Model for solar module BP3125 with a maximum power point tracker.
00191 //
00192 // solarEnergy is the percentage standard incident solar radiation used to
00193 // define the module characteristics (ie 1000 W/m2).
00194 //
00195 // This uses a simple hill-climbing search for maximum power starting at the
00196 // open-circuit voltage and stepping back to the peak.
00197 
00198 double OptimalModuleVoltage(const double solarEnergy)
00199 {
00200     if (solarEnergy == 0) return 0;
00201     const double Isc = 8.02;      // Short Circuit Current
00202     const double I0 = 0.000185;   // Diode dark current
00203     const double Vk = 2.071;      // Model parameter voltage
00204     double b = Isc*solarEnergy*0.01/I0+1;
00205     double Voc = Vk*log(b);       // Open Circuit voltage
00206     double Vinc = Voc/10;         // Initial increment
00207     double V = Voc;               // Stepping back from here
00208     double powerLast = 0;         // previous power computation
00209     double power;
00210     for (int j = 0; j < 10; j++)
00211     {
00212         bool finished = false;
00213         while (! finished)
00214         {
00215             V -= Vinc;
00216             power = V*I0*(b-exp(V/Vk));
00217             if (power < powerLast) finished = true;
00218             powerLast = power;
00219         }
00220         V += Vinc;
00221         Vinc /= 10;
00222     }
00223 // Buck regulator only
00224 //    if (V < 12) V = 12;       // MPP is below battery voltage
00225 //    if (Voc < 12) V = 0;      // Highest voltage is below battery
00226     return V;
00227 }
00228 //----------------------------------------------------------------------------
00229 // Compute the density of air from a suitable model
00230 //
00231 // Glenn Research Centre http://www.grc.nasa.gov/WWW/K-12/airplane/atmos.html
00232 // The units of the model are given in imperial, so we convert to metric last
00233 //
00234 // Input:  height in metres
00235 // Output: air density in kg/m^3
00236 //
00237 double airDensity(const double height)
00238 {
00239     double Ta;                               // Absolute temperature
00240     double pressure;                         // pressure lbs/ft^2
00241     if (height < 11019)
00242     {
00243         Ta = 288.2 - 0.00649 * height;
00244         pressure = 10331 * pow(0.003471*Ta,5.256);
00245     }
00246     else if (height < 25099)
00247     {
00248         Ta = 216.5;
00249         pressure = 2309.9 * exp(1.73-0.00015748*height);
00250     }
00251     else
00252     {
00253         Ta = 141.5 + 0.00299 * height;
00254         pressure = 253.39 * pow(0.0046*Ta,-11.388);
00255     }
00256     return pressure*0.0341636/Ta;           // density in kg/m^3
00257 }
00258 //----------------------------------------------------------------------------
00259 // Numerical integration of air density over a sloping solar ray path.
00260 //
00261 // The angle of the path to the vertical is phi. To integrate density rho(h)
00262 // as a function of height, over the path length a, we compute the slope of the
00263 // path with respect to the height from ground level to get:
00264 //
00265 //      integral from 0 to infinity rho(h) (da/dh) dh
00266 //
00267 // Input:  cosine of angle of path to vertical phi in degrees
00268 // Output: Path loss. Units are arbitrary as this appears only in ratios.
00269 //
00270 double pathLoss(const double cosPhi)
00271 {
00272     const double R = 6335437;                   // earth radius in metres
00273     double hIncr = 10;                          // Integration increment m
00274     double h = hIncr;                           // height above sea level
00275 // Because of numerical problems near h=0, cosPhi=0 (tangential incidence)
00276 // we integrate over the first height step using constant density equal
00277 // to the average of the step (trapezoidal approximation)
00278     double loss = 0.5*hIncr*(airDensity(h) + airDensity(0))
00279                 * (2*R + h)*h/(R*cosPhi+sqrt(R*R*cosPhi*cosPhi+2*h*R + h*h));
00280 // This starts off the trapezoidal approximation (see notes)
00281 // 99.999% of air mass is below 100km, so we stop iteration there.
00282 // As density contribution falls, increase increment to speed up things.
00283     loss += 0.5*hIncr*airDensity(h)*(R + h)/
00284                 sqrt(R*R*cosPhi*cosPhi+2*h*R + h*h);
00285     while (h < 100000)
00286     {
00287         if (h > 6000) hIncr = 20;
00288         else if (h > 10000) hIncr = 50;
00289         else if (h > 16000) hIncr = 100;
00290         h += hIncr;
00291         loss += hIncr*airDensity(h)*(R + h)/
00292                     sqrt(R*R*cosPhi*cosPhi+2*h*R + h*h);
00293     }
00294     return loss;
00295 }
00296 //----------------------------------------------------------------------------
00297 // Integration of the solar power W/m^2 to kWH/m^2 for a day at given latitude
00298 //
00299 // Integration is done by a simple sum over small elements given in minute
00300 // increments and over half a day. The result is converted to kWH by dividing
00301 // by 60,000 and multiplying by 2 for the second half of the day.
00302 // Represents a module always facing the sun.
00303 // Input: Latitude in degrees, positive north of equator
00304 //        Declination of the sun in degrees
00305 // Output: Total energy per square metre over a day arriving at earth's surface
00306 // Dependencies: pathLoss(cosangle) integral of air density over a slant path
00307 
00308 double dailySolarEnergyFollowing(const double latitude,
00309                                  const double declination)
00310 {
00311     const double angleConversion = 3.1415927/180.0;
00312     const double rDeclination = declination*angleConversion;
00313     const double cosDeclination = cos(rDeclination);
00314     const double sinDeclination = sin(rDeclination);
00315     const double rLatitude = latitude*angleConversion;
00316     const double cosLatitude = cos(rLatitude);
00317     const double sinLatitude = sin(rLatitude);
00318     const double lossConstant = 0.253617853024622586/pathLoss(1);
00319     const double solarConstant = 1366;              // W/m^2 outer atmosphere
00320     double minuteIncr = 1;                          // integration step size
00321     double minute = 0;
00322     double cosAngle = 1;
00323     double solarEnergy = 0;
00324     while (cosAngle > 0)
00325     {
00326         double cosHourAngle = cos(0.25*minute*angleConversion);
00327         cosAngle = cosLatitude*cosDeclination*cosHourAngle
00328                     + sinLatitude*sinDeclination;
00329         solarEnergy += solarConstant*exp(-lossConstant*pathLoss(cosAngle));
00330         minute += minuteIncr;
00331     }
00332     return solarEnergy/30000;
00333 }
00334 //----------------------------------------------------------------------------
00335 // Integration of the solar power W/m^2 to kWH/m^2 for a day at given latitude
00336 //
00337 // Integration is done by a simple sum over small elements given in minute
00338 // increments and over half a day. The result is converted to kWH by dividing
00339 // by 60,000 and multiplying by 2 for the second half of the day.
00340 // Represents a fixed module facing the sun at noon.
00341 // Input: Latitude in degrees, positive north of equator
00342 //        Declination of the sun in degrees
00343 //        Angle of the module to the equatorial plane
00344 // Output: Total energy per square metre over a day arriving at earth's surface
00345 // Dependencies: pathLoss(cosangle) integral of air density over a slant path
00346 
00347 double dailySolarEnergyFixed(const double latitude,
00348                              const double declination,
00349                              const double moduleAngle)
00350 {
00351     const double angleConversion = 3.1415927/180.0;
00352     const double rDeclination = declination*angleConversion;
00353     const double cosDeclination = cos(rDeclination);
00354     const double sinDeclination = sin(rDeclination);
00355     const double rModuleAngle = moduleAngle*angleConversion;
00356     const double cosModuleAngle = cos(rModuleAngle);
00357     const double sinModuleAngle = sin(rModuleAngle);
00358     const double rLatitude = latitude*angleConversion;
00359     const double cosLatitude = cos(rLatitude);
00360     const double sinLatitude = sin(rLatitude);
00361     const double lossConstant = 0.253617853024622586/pathLoss(1);
00362     const double solarConstant = 1366;              // W/m^2 outer atmosphere
00363     double minuteIncr = 1;                          // integration step size
00364     double minute = 0;
00365     double cosAngle = 1;
00366     double solarEnergy = 0;
00367     while (cosAngle > 0)
00368     {
00369         double cosHourAngle = cos(0.25*minute*angleConversion);
00370         cosAngle = cosLatitude*cosDeclination*cosHourAngle
00371                  + sinLatitude*sinDeclination;
00372         double cosIncidence = cosModuleAngle*cosDeclination*cosHourAngle
00373                             + sinModuleAngle*sinDeclination;
00374         if (cosIncidence > 0)
00375             solarEnergy += solarConstant*cosIncidence*
00376                            exp(-lossConstant*pathLoss(cosAngle));
00377         minute += minuteIncr;
00378     }
00379     return solarEnergy/30000;
00380 }
00381 //----------------------------------------------------------------------------
00382 // Length of day in hours for given latitude and solar declination
00383 // There is 15 degrees per hour movement of the sun.
00384 // Input: Latitude in degrees, positive north of equator
00385 //        Declination of the sun in degrees
00386 // Output: length of day sunrise to sunset in hours.
00387 //
00388 double dayLength(const double latitude, const double declination)
00389 {
00390     const double angleConversion = 3.1415927/180.0;
00391     double rLatitude = latitude*angleConversion;
00392     const double rDeclination = declination*angleConversion;
00393     return 2*acos(-tan(rLatitude)*tan(rDeclination))/(15*angleConversion);
00394 }
00395 //----------------------------------------------------------------------------
00396 // Tests air density printout
00397 
00398 //    for (double h=0;h<30000;h+=1000)
00399 //    {
00400 //       std::cout << h << " " << airDensity(h) << std::endl;
00401 //    }
00402 
00403 // Solar power variation sun directly overhead at noon
00404 
00405 //    double cosangle = 1;
00406 //    double hour = 12;
00407 //    while (cosangle > 0)
00408 //    {
00409 //        double solarEnergy =
00410 //              solarConstant*exp(-lossConstant*pathLoss(cosangle));
00411 //        std::cout << hour << "," << solarEnergy << std::endl;
00412 //        hour += 0.1;
00413 //        cosangle = cos(angleConversion*15*(hour-12));
00414 //    }
00415 
00416 // Total daily energy from modules over latitudes
00417 
00418 //    double declination = maxDeclination;
00419 //    for (float n = -60; n < 60; n++)                // Range over latitudes
00420 //    {
00421 //        double latitude = n;
00422 //        std::cout << latitude << ","
00423 //                  << dayLength(latitude,maxDeclination) << ","
00424 //                  << dailySolarEnergyFollowing(latitude,declination) << ","
00425 //                  << dailySolarEnergyFixed(latitude,declination,declination)
00426 //                  << ","
00427 //                  << dailySolarEnergyFixed(latitude,declination,latitude)
00428 //                  << std::endl;
00429 //    }
00430 
00431 // Daily Global Solar Radiation at Armidale
00432 
00433 //    double latitude = -30.5;
00434 //    for (int month = 0; month < 12; month++)
00435 //    {
00436 //        double declination = maxDeclination*sin(2*3.1415927*(month-2.25)/12);
00437 //        std::cout << dailySolarEnergyFixed(latitude,declination,latitude)
00438 //                  << std::endl;
00439 //    }
00440 
00441 //    for (uint v = 0; v < 230; v++) std::cout << (double)v/10 << ","
00442 //         << moduleCurrent(100,(double)v/10) << std::endl;
00443 /*        std::cout << minute << ","
00444                   << solarEnergyFixed*energyCharge << ","
00445                   << solarEnergyFollowing*energyCharge << ","
00446                   << current12fixed << ","
00447                   << current12following << ","
00448                   << Vfixed << "," << currentVfixed << ","
00449                   << Vfollowing << "," << currentVfollowing << ","
00450                   << std::endl;
00451 */

Generated on Mon Sep 3 06:31:15 2007 for Solar Power by  doxygen 1.5.2