Von der Globalstrahlung zur Ertragsvorhersage

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 und damit den Ertrag 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(γ))​

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