Von der Globalstrahlung zur Ertragsvorhersage

Dieser Artikel erläutert die grundlegenden Berechnungen, um aus einer Angabe der Globalstrahlung auf einer horizontalen Fläche die Strahlungskomponenten der direkten Strahlung B, der diffusen Strahlung D und der durch Bodenreflexion hinzukommenden Strahlung R zu erhalten.

Berechnung der Einstrahlung

Gegeben seien die folgenden Größen

• Sg = Intensität der Globalstrahlung (diffuse und direkte Strahlung) auf einer horizontalen Fläche in W/m². Achtung: Die Vorhersagedaten des DWD sind in der Einheit kJ/m²h angegeben und müssen durch 3.6 dividiert werden, um die Intensität in W/m² zu erhalten
• γs = 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(γs − γ)

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(β))​

Hay-Davies

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. Allerdings muss man ggf. berücksichtigen, dass direkte Strahlung nur von "vorne" kommt, diffuse Strahlung aber ggf. auch von der anderen Seite.

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(β)))

Perez

Das Modell von Perez ist deutlich komplizierter, weil es die diffuse Einstrahlung in einen Anteil in Sonnennähe und einen isotropen Anteil zerlegt, aber zusätzlich eine Horizontaufhellung berücksichtigt. Im Endeffekt wird die diffuse Einstrahlung auf die PV-Fläche nach Perez berechnet als

D = Sf * (C1 ​+ (1−C2)*0.5*(1+cos(β)) + C3*sin(β))

Allerdings ist die Berechnung der Koeffizienten C1..C3 sehr kompliziert, erfordert u.a. eine Wertetabelle. Außerdem ist für PV-Anlagen am Horizont mit Bebauung etc. zu rechnen, so dass die Verwendung dieses Modells derzeit nicht sinnvoll scheint.

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

Ig = B + D + R

Einfluss der Moduloberfläche

Bis zur Ertragsvorhersage sind noch mehrere Schritte zu unternehmen. Denn Solarmodule weisen unterschiedliche Empfindlichkeiten für die verschiedenen Strahlungsarten auf. Formal wird dies über so genannte 'Incident Angle Modifier' IAM berücksichtigt, die jeweils Funktionen des oben berechneten Einfallswinkels θ sind. Das gebräuchlichste Modell dafür stammt von der American Society of Heating, Refrigerating and Air-Conditioning Engineers ASHRAE und lautet

IAM(θ) = 1 - b0*(1/cos(θ) - 1)

mit einem modulabhängigen Parameter b0 = 0.05 .. 0.2 je nach Glasoberfläche der Module. Bei θ=0 ist dieser Faktor also 1 und fällt bei größeren Einfallswinkeln ab. Wichtig: IAM darf nicht negativ werden.

Diesen IAM kann man für die direkte Strahlung verwenden. Bei diffuser Strahlung und Bodenreflexion kommen prinzipiell alle Einfallswinkel vor, so dass man zwar dieselbe Funktion verwenden kann, allerdings mit jeweils einem effektiven Einfallswinkel. Gebräuchliche Näherungen für diese effektiven Einfallswinkel sind

θd = 59.7° − 0.1388*β + 0.001497*β²
θr = 90° − 0.5788*β + 0.002693*β²

mit β in Grad.

Die tatsächlich in elektrische Energie umsetzbare Strahlung ist also

Ieff = IAM(θ)*B + IAM(θd)*D + IAM(θr)*R

Eine brauchbare Näherung ergibt sich bereits, wenn sowohl IAM(θd) als auch IAM(θr) auf 1 gesetzt werden.

Berechnung des Ertrags

Die elektrische Leistung aus der PV-Anlage ergibt sich, indem die berechnet effektive Einstrahlung Ieff mit der Modulfläche Am und dem Wirkungsgrad η multipliziert wird. Dieser Wirkungsgrad hängt stark von der Temperatur Tm der Module ab, er sinkt gegenüber dem Referenzwert η0 bei 25 °C um jeweils ca. 1% pro Kelvin ab.

Die zu erwartende elektrische Leistung ist also

P = Ieff * Am * η0 * (1 - 0.01*(Tc - 25))

Perl-Programm

Ein einfaches Perl-Programm zur Berechnung folgt hier. Es bekommt 8 Parameter übergeben:

• sg => Sg = Intensität der Globalstrahlung (diffuse und direkte Strahlung) auf einer horizontalen Fläche. Achtung: Die vom DWD angegebenen Daten haben die Einheit kJ/m²h. Um diese in W/m² umzurechnen, müssen sie durch 3.6 dividiert werden.
• sunazi => γs = 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
• b0 => b0 Materialfaktor für IAM, 0.05 .. 0.2 (Startwert für eine Vorhersage sollte 0.05 sein)

Rückgabewert ist die effektive Einstrahlung auf die PV-Anlage in W/m²

sub pfx_PVforecast($$$$$$$$) {
  my ($sg, $sunazi, $sunele, $day, $strazi, $strtilt, $model, $b0) = @_;

  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);
  $b0    = 0.05 if !defined($b0);

  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;

  # IAM nur für Direktstrahlung (ASHRAE)
  my $IAMb = 1.0;
  if ($B_tilt > 0.0 && $cos_thetai_pos > 0.0) {
    $IAMb = 1.0 - $b0 * (1.0 / $cos_thetai_pos - 1.0);
    $IAMb = 0.0 if $IAMb < 0.0;
    $IAMb = 1.0 if $IAMb > 1.0;
  }

  my $B_eff = $B_tilt * $IAMb;

  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;
  }

[[Kategorie:Energieerzeugungsmessung]]
  $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_eff + $D_tilt + $R_tilt;
  $G_tilt = 0.0 if $G_tilt < 0.0;

  return $G_tilt;
}