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

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 derzei tnicht 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

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