Von der Globalstrahlung zur Ertragsvorhersage
Berechnung der Einstrahlung
Gegeben seien die folgenden Größen
- Sg = Intensität der Globalstrahlung (diffuse und direkte Strahlung) auf einer horizontalen Fläche
- α = Azimutwinkel des Sonnenstandes, 180° = Südrichtung
- β = Höhenwinkel des Sonnenstandes, 0° = waagerecht, 90° = senkrecht
- n = Tag des Jahres, 1..365
- φ = Azimutwinkel der PV-Anlage, 0° = Nordrichtung, 180° = Südrichtung, gemessen im Uhrzeigersinn
- γ = Neigungswinkel der PV-Anlage, 0° = horizontal, 90° = senkrecht
Die Aufgabe besteht nun darin, aus diesen Daten die Globalstrahlung auf der geneigten PV-Anlage zu berechnen.
Als Zwischengrößen benötigen wir
- Si = Intensität der direkten Strahlung auf einer horizontalen Fläche
- Sf = Intensität der diffusen Strahlung auf einer horizontalen Fläche
womit also gilt
Sg = Sf + Si
Klarheitsindex und diffuser Anteil
Im nächsten Schritt schätzen wir aus der Globalstrahlung einen Klarheitsindex Kt ab, indem wir zunächst mit der Solarkonstante
I0 = 1367 W/m²
die außeratmosphärische normale Einstrahlung I0n modellieren. Diese beträgt am n-ten Tag des Jahres ungefähr
I0n = I0 * ( 1 + 0.033 * cos( 360° * n/365))
Horizontal ergibt sich damit
I0h = I0n*sin(β)
für β > 0. Der Klarheitsindex ergibt sich als
Kt = Sg/I0h
Um Instabilitäten bei kleinen Sonnenwinkeln zu vermeiden, sollte dieser Wert auf das Intervall [0,0.85] geclippt werden.
Der diffuse Anteil nach Erbs ist dann
Sf = kd * Sg
wobei der Wert kd sich wie folgt aus Kt ergibt:
- Für Kt <= 0.22 => kd = 1-0.09*Kt
- Für 0.22 < Kt <= 0.80 => kd = 0.9511−0.1604*Kt+4.388*Kt²−16.638*Kt³+12.336*Kt⁴
- Für 0.80 <= Kt => kd = 0.165
Der direkte Strahlungsanteil auf der horizontalen Fläche ist also
Si = (1 - kd) * Sg
Direkte Einstahlung auf die PV-Fläche
Die direkte Normalstrahlung auf horizontaler Fläche ist
DNI = Si * 1/sin(β)
Der Kosinus des Einfallswinkels auf die Modulfläche ergibt sich als
cos(θ)=sin(β)*cos(γ)+cos(β)*sin(γ)*cos(α − φ)
Dieser Wert wird negativ, wenn sich die Sonne hinter der PV-Fläche befindet.
- Wenn es ältere Module mit undurchsichtiger Rückseite sind, entfällt in diesem Fall die direkte Einstrahlung. Als Faktor für die weitere Rechnung ist also zu verwenden
F = max(0,cos(θ))
- Wenn es sich um eine Anlage mit modernen Glas-Glas-Modulen handelt, ist die direkte Einstrahlung auch dann noch vorhanden. Als Faktor für die weitere Rechnung ist bei gleicher Empfindlichkeit aus beiden Richtungen also zu verwenden
F = abs(cos(θ))
- Im Zweifelsfall ist hier ein Korrekturfaktor anzubringen, wenn cos(θ) < 0
Die direkte Strahlung auf die PV-Anlage ergibt sich damit als
B = Si * F/sin(β)
für β > 0. Zur besseren Absicherung gegen nummerische Instabilität sollte das aber tatsächlich nur verwendet werden wenn β > 3°, ansonsten hart auf Null gesetzt werden.
Diffuse Einstrahlung auf die PV-Fläche
Das so genannte isotrope Himmelsmodell ergibt für die diffuse Einstrahlung auf die PV-Fläche
D = Sf * 0.5*(1+cos(γ))
Besser geeignet ist ein Modell nach Hay-Davies, in welchem die diffuse Einstrahlung in einen Anteil in Sonnennähe und einen isotropen Anteil zerlegt wird. Dafür benötigen wir den Zenitwinkel des Sonnenstandes, dieser ist gleich 90°-β, um einen geometrischen Faktor r zu berechnen
r = F/cos(90°-β) = F/sin(β)
Darin ist F der im vorigen Abschnitt berechnete Faktor, also in der Regel der Kosinus des Einfallswinkels der direkten Strahlung auf die PV-Anlage.
Desweiteren bestimmen wir einen Anisotropiefaktor A
A = DNI/I0n = Si/I0n * 1/sin(β) = Si/I0h = Si/Sg * Kt = (1 - kd)* Kt
Damit ergibt sich dann als diffuse Einstrahlung auf die PV-Fläche nach Hay-Davies
D = Sf * (A*r + (1−A) * 0.5*(1+cos(γ)))
Bodenreflexion
Bei geeigneter Anlage kann man noch die Bodenreflexion hinzufügen, mit einer Albedo
- ρ = 0.2 für "normalen" Boden
- ρ = 0.6 .. 0.8 bei Schnee
ergibt sich die Einstrahlung aus Bodenreflexion als
R = ρ * Sg * 0.5*(1-cos(γ))
Gesamtergebnis für die Einstrahlung
Die Einstrahlung auf die PV-Anlage ist also
I = B + D + R
Ertragsvorhersage
Bis zur Ertragsvorhersage sind noch mehrere Schritte zu unternehmen. Denn Solarmodule weisen unterschiedliche Empfindlichkeiten für die verschiedenen Strahlungsarten auf.
Perl-Programm
Ein einfaches Perl-Programm zur Berechnung folgt hier. Es bekommt 7 Parameter übergeben:
- sg => Sg = Intensität der Globalstrahlung (diffuse und direkte Strahlung) auf einer horizontalen Fläche
- sunazi => α = Azimutwinkel des Sonnenstandes, 180° = Südrichtung
- sunele => β = Höhenwinkel des Sonnenstandes, 0° = waagerecht, 90° = senkrecht
- day => n = Tag des Jahres, 1..365
- strazi => φ = Azimutwinkel der PV-Anlage, 0° = Nordrichtung, 180° = Südrichtung, gemessen im Uhrzeigersinn
- strtilt => γ = Neigungswinkel der PV-Anlage, 0° = horizontal, 90° = senkrecht
- model => 0: Berechnung mit isotropem Himmelsmodell, 1: Berechnung nach Hay-Davies
Rückgabewert ist die Einstrahlung auf die PV-Anlage in W/m²
sub pfx_PVforecast($$$$$$$) {
my ($sg, $sunazi, $sunele, $day, $strazi, $strtilt, $model) = @_;
my $rho = 0.2; # Bodenalbedo
my $pi = 4 * atan2(1,1);
my $deg = $pi / 180.0;
my $Isc = 1367.0; # Solarkonstante W/m²
return 0 if !defined($sg) || $sg <= 0;
return 0 if !defined($sunele) || $sunele <= 0;
$model = 0 if !defined($model);
my $sunazi_r = $sunazi * $deg;
my $sunele_r = $sunele * $deg;
my $strazi_r = $strazi * $deg;
my $strtilt_r = $strtilt * $deg;
my $sin_ele = sin($sunele_r);
my $cos_ele = cos($sunele_r);
return 0 if $sin_ele <= 0.0;
my $I0n = $Isc * (1.0 + 0.033 * cos(2.0 * $pi * $day / 365.0));
my $I0h = $I0n * $sin_ele;
return 0 if $I0h <= 0.0;
my $Kt = $sg / $I0h;
$Kt = 0.0 if $Kt < 0.0;
$Kt = 1.0 if $Kt > 1.0;
my $kd;
if ($Kt <= 0.22) {
$kd = 1.0 - 0.09 * $Kt;
}
elsif ($Kt <= 0.80) {
$kd = 0.9511
- 0.1604 * $Kt
+ 4.388 * $Kt**2
- 16.638 * $Kt**3
+ 12.336 * $Kt**4;
}
else {
$kd = 0.165;
}
$kd = 0.0 if $kd < 0.0;
$kd = 1.0 if $kd > 1.0;
my $sf = $kd * $sg; # diffuse horizontale Strahlung = DHI
my $si = $sg - $sf; # direkte horizontale Strahlung = BHI
$si = 0.0 if $si < 0.0;
my $cos_thetai =
$sin_ele * cos($strtilt_r)
+ $cos_ele * sin($strtilt_r) * cos($sunazi_r - $strazi_r);
my $cos_thetai_pos = $cos_thetai;
$cos_thetai_pos = 0.0 if $cos_thetai_pos < 0.0;
my $dni = ($sin_ele > 0.01) ? ($si / $sin_ele) : 0.0;
my $B_tilt = $dni * $cos_thetai_pos;
$B_tilt = 0.0 if $B_tilt < 0.0;
my $D_tilt;
if ($model == 1) {
my $Ai = ($sg > 0.0) ? (($si / $sg) * $Kt) : 0.0;
$Ai = 0.0 if $Ai < 0.0;
$Ai = 1.0 if $Ai > 1.0;
my $Rb = ($sin_ele > 0.01) ? ($cos_thetai_pos / $sin_ele) : 0.0;
$Rb = 0.0 if $Rb < 0.0;
$D_tilt = $sf * (
$Ai * $Rb
+ (1.0 - $Ai) * (1.0 + cos($strtilt_r)) / 2.0
);
}
else {
$D_tilt = $sf * (1.0 + cos($strtilt_r)) / 2.0;
}
$D_tilt = 0.0 if $D_tilt < 0.0;
my $R_tilt = $rho * $sg * (1.0 - cos($strtilt_r)) / 2.0;
$R_tilt = 0.0 if $R_tilt < 0.0;
my $G_tilt = $B_tilt + $D_tilt + $R_tilt;
$G_tilt = 0.0 if $G_tilt < 0.0;
return $G_tilt;
}