Projekt: Optimal opsætning af solpaneler

Projekt: Optimal opsætning af solpaneler#

Indledning#

En husejer nær DTU ønsker at opsætte solpaneler på sit parcelhus med fladt tag. Hun tænker hun skal vende dem mod syd, men er i tvivl om hvilken vinkel panelet bør danne med taget. I dette projekt skal I finde den optimale vinkel for husejerens solpaneler. Projektets mål er at udvikle en model og et Python-program der for angivet placering på jorden (længde- og breddegrad) kan udregne den optimale vinkel for opsætning af solpaneler under en række simplificerende antagelser. Til simulering skal bruges data fra kalenderåret 2024. I starter med at modellere solens bevægelse på himlen ved hjælp af Python-pakken pvlib. I skal herefter opstille en model for solpanelets energiproduktion, der ikke tager højde for skyer og atmosfære. Solen modelleres som en vektorfelt af parallelle vektorer. For at udregne energiproduktionen skal I kunne integrere vektorfeltets flux gennem solpanelfladen. Vi antager først at panelet er fladt, men I skal efterfølgende generalisere jeres model. En mulighed generalisering er følgende: Find en passende flade i bymiljøet omkring DTU, som kan udstyres med tynde, fleksible solpanelsfilm (som https://www.youtube.com/watch?v=TS9ADU0oc50) og udregn energiproduktionen fra sådanne løsninger.

Forberedelse#

Om projektet#

Projektmål#

Det overordnede mål med projektet at lave et Python-program der for angivet placering på jorden (længde- og breddegrad) kan udregne den optimale vinkel for solpaneler.

Optimalitet#

Optimal kan her betyde mange ting:

  • Størst årligt energiproduktion. Det er dette optimalitets mål vi bruger med mindre andet opgives.

  • Krav om en vis daglig minimumslevering (fx om vinteren).

  • Laveste samlede energiomkostninger for boligen. Her skal medregnes forbrugskurver, salg af energi til elnettet, en eventuel elbil, m.m.

  • Osv.

Teori: Model og Antagelser#

Solenergi er omdannelsen af ​​energi fra sollys til elektricitet ved hjælp af solceller (photovoltaic (PV) cells). Solpanelers produktion af strøm afhænger af anlæggets størrelse (antallet af paneler), panelernes individuelle effektivitet, panelernes placering og selvfølgelig vejrforhold.

Solindstråling eller irradians (solar irradiance) er den effekt pr. arealenhed (overfladeeffekttæthed) som modtages fra solen i form af elektromagnetisk stråling i måleinstrumentets bølgelængdeområde. Solindstråling måles i watt pr. kvadratmeter (\(W/m^2\)) i SI-enheder. Solens indstråling har forskellig bølgelængde, fx UV-stråling, der samlet kaldes indstrålingens spektrum.

Den gennemsnitlige årlige solstråling, der ankommer til toppen af ​​jordens atmosfære, er omkring \(1361 W/m^2\). Når solens stråler har passeret atmosfæren, er den samlede stråling (irradians) i Danmark ca. \(S_0 := 1100 W/m^2\) på en klar dag midt på dagen. Spektrummet ændres ligeledes gennem atmosfæren, fx kommer den korteste UV-stråling (under 280 nanometer) sjældent ned til Jordens overflade, men vi vil ignore sådanne spektrale ændringer.

Man kan med god tilnærmelse regne med en lineær sammenhæng mellem solens indstråling og den strømstyrke solpanelet leverer, og vi vil bruge denne tilnærmelse. Da temperaturen i praksis stiger ved øget indstråling, og spændingen falder med øget temperatur, vil effekten dog ikke stige helt lineært. Solens stråler modelleres som et vektorfelt \(\pmb{V}_S\) af parallelle vektorer af længde \(S_0\). Lad os antage at solpanelet er beskrevet ved en flade \(\mathcal{F} = \pmb{r}(\Gamma)\), hvor \(\pmb{r} : \Gamma \to \mathbb{R}^2\), \(\Gamma = [a_1,b_1] \times [a_2,b_2]\), er en parametrisering af fladen. Energiproduktionen for en solcelle afhænger af vinklen mellem fladens normalvektor \(\pmb{n}_\mathcal{F}\) og solstrålerne. Hvis solstrålerne er parallelle med normalvektoren er solens effekt maksimal, mens hvis solstrålerne er vinkelrette på normalvektoren er effekten nul. Solens effekt pr areal i punktet \(\pmb{r}(u,v)\), hvor \((u,v) \in \Gamma\), er projektionen af solens stråler ind på fladens enheds normalvektor:

\[\begin{equation*} \left\langle \pmb{V}, \frac{\pmb{n}_\mathcal{F}(u,v)}{\Vert \pmb{n}_\mathcal{F}(u,v) \Vert} \right\rangle \end{equation*}\]

Den samlede effekt fås at integrere denne størrelse op over held fladen, hvilet svarer til fluxen af \(V\) gennem fladen \(\mathcal{F}\):

\[\begin{equation*} \int_{\Gamma} \langle \pmb{V}, \pmb{n}_\mathcal{F}(u,v) \rangle \mathrm{d}(u,v) \end{equation*}\]

Hvis \(\langle \pmb{V}, \pmb{n}_\mathcal{F}(u,v) \rangle < 0\) betyder det er solen lyser på bagsiden af solpanelet, og effekten skal derfor sættes til nul (effekten kan aldrig blive negativ).

Vi angiver her vores standard-antagelser, som bruges med mindre andet angives:

  • Panelet er fladt og fastmonteret

  • Vi antager at den maksimale irradians/indstråling fra solen på solpanelet er \(S_0 = 1100 W/m^2\). Alle vektorer i solens vektorfelt \(\pmb{V}: \mathbb{R}^3 \to \mathbb{R}^3\) er parallelle og har længden \(S_0\).

  • Solpanelets effekt afhænger lineært af solindstrålingen og fluxen af solens vektorfelt gennem solpanelets (over)flade.

  • Solpanelet effektivitet: Ved en flux på 1000 W per \(m^2\) leverer panelet \(Wp/(L B)\) per \(m^2\), hvor \(Wp\) er panelets peak power watts (STC), \(L\) er længden og \(B\) er bredden.

  • Skydækket og andre atmosfæriske forstyrrelser modelleres ikke, men fastsættes til en gennemsnitlig værdi: vi antager at disse forstyrrelser halverer solens energi gennem hele året: \(S_0 A_0 = 550 W/m^2\), hvor faktoren er \(A_0 = 0.5\).

Der vil i sidste del af projektet være mulighed for at vende tilbage til disse antagelser og undersøge om/hvordan de kan erstattes af bedre antagelser.

Indledende øvelser#

Find via litteratur-søgning den anbefalede vinkel man opstiller solpaneler ved i Danmark. Vinklen siges at være er nul grader hvis panelet ligger fladt ned på jorden (eller taget).

Udvælg en solpanelstype. I kan søge på montører i Danmark og undersøge hvilke paneler de typisk bruger, eller I kan google “solar panel datasheet” eller lignende. Solpanelet bør være et standard panel (der er fladt og altså ikke fx krummer). Find et datablad for det valgte solpanel, og beskriv panelets størrelse (anvend \(L\) længde og \(B\) for bredde) og angiv Wp/Pmax (kaldet max power, peak power watts eller lignende) under standarden STC. Beskriv hvad standarden STC beskriver. Udregn \(Wp/(L B)\) (jvf. listen af standard-antagelser ovenfor). Antag ideelle forhold: solen står vinkelret på solpanelet og solens irradians er \(1100 W/m^2\) igennem en hel time. Hvor mange J og kWh leverer panelet på denne time? Hvor mange Wh er dette per \(m^2\)?

Jeres endelige rapport skal desuden indeholde et indledende afsnit, hvor I beskriver solenergi og solceller. I bestemmer selv hvad afsnittet præcist skal indeholde, men nedenstående kilder kan bruges til at findes information. Det anbefales først at skrive dette afsnit, når I nået længere med projektet.

  1. https://www.pveducation.org/

  2. https://www.acs.org/education/resources/highschool/chemmatters/past-issues/archive-2013-2014/how-a-solar-cell-works.html

Panelet er fladt, så enhedsnormalvektoren \(\pmb{u}_p\) er konstant (og har længde 1) hvor:

\[\begin{equation*} \pmb{u}_p = \frac{\pmb{n}_{\mathcal{F}}(u,v)}{\Vert \pmb{n}_\mathcal{F}(u,v) \Vert} \end{equation*}\]

Angiv en formel eller udtryk for fluxen gennem fladen udtrykt ved \(A_0\), \(\pmb{u}_p, \pmb{V}, L\) og \(B\). Vi definerer fluxen til at være nul hvis vinklen mellem \(\pmb{u}_p\) og \(\pmb{V}\) er større end \(\pi/2\) (90 grader), da vi ikke ønsker negativ flux. Jeres udtryk skal tage højde for dette.

Hint: Da normalvektoren er konstant, kan I slippe af med integraltegnet.

Fluxen bestemmer den øjeblikkelige effekt for panelet for en bestemt retning af solens stråler \(\pmb{V}\) og for en bestemt placering af panelet bestemt ved \(\pmb{u}_p\). For at finde den samlede energiproduktion for en solpanel, har vi har derfor brug for at kunne modellere solens bevægelse i forhold til panelet. Dette modeleres med en solspositionsalgoritme (SPA) fra Python-pakken pvlib, men inden vi kommer så langt har vi brug at udvikle nogle værktøjer i NumPy.

Angiv SI-enheder for \(\pmb{V}\), \(L\), \(B\), \(A_0\), fluxen, energi. Angiv sammenhængen mellem J og kWh.

NumPy#

I projektet har vi brug for at kunne finde minimum, maksimum og/eller nulpunkter as NumPy-arrays, så vi starter med at blive fortrolige med dette. Vi betragter en simpel funktion \(f: [0, 2 \pi] \to \mathbb{R}\) givet ved \(f(x) = \cos(x)\) og ønsker at kunne bestemme minimum, maksimum og nulpunkter af denne. Vi vil i projektet dog ikke have selve funktionen eller dens funktionsudtryk, men kun en (lang) række at funktionsværdier, fx \(f(t_n)\) for \(t_n = 2\pi n/N\) for \(n = 0,1, \dots, N-1\). Vi opfatter her \(\Delta t_n = 2 \pi/N\) som afstanden mellem vores tids-samples. I NumPy:

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
figurstr = (15, 7)

t = np.linspace(0, 2 * np.pi, 1000)
f = np.cos(t)

plt.figure(figsize=figurstr)
plt.plot(t, f, color='b', linestyle='-')
plt.xlabel('t')
plt.ylabel('cos(t)')
plt.show()

# print(t,f)  # uncomment to see the values of t and f
../_images/c645fd7f461e380a6e76414d9eb58a2b6b790d562fc9c938c63b0fb97d7f65bb.png

Vi opfatter t og f som vektorer af hhv variabel- og funktionsværdier. NumPy har mange indbyggede metoder, der kan være nyttige, fx

f.max(), f.min(), f.argmax(), f.argmin(), t[f.argmax()], t[f.argmin()]

(1.0, -0.9999950553174459, 0, 499, 0.0, 3.138447916198812)

Man kan også spørge hvilke funktionsværdier der fx er mindre end \(-0.95\):

idx = f < -0.95
f[idx], t[idx]

# Or in one line
f[f < -0.95], t[f < -0.95]

(array([-0.95192731, -0.95383508, -0.95570513, -0.95753737, -0.95933173,
        -0.96108814, -0.96280654, -0.96448685, -0.96612901, -0.96773295,
        -0.96929861, -0.97082592, -0.97231483, -0.97376528, -0.97517722,
        -0.97655057, -0.9778853 , -0.97918134, -0.98043865, -0.98165717,
        -0.98283687, -0.98397768, -0.98507957, -0.9861425 , -0.98716641,
        -0.98815128, -0.98909706, -0.99000371, -0.9908712 , -0.99169949,
        -0.99248855, -0.99323836, -0.99394887, -0.99462007, -0.99525192,
        -0.9958444 , -0.99639749, -0.99691116, -0.9973854 , -0.99782019,
        -0.9982155 , -0.99857133, -0.99888765, -0.99916446, -0.99940175,
        -0.99959951, -0.99975772, -0.99987639, -0.9999555 , -0.99999506,
        -0.99999506, -0.9999555 , -0.99987639, -0.99975772, -0.99959951,
        -0.99940175, -0.99916446, -0.99888765, -0.99857133, -0.9982155 ,
        -0.99782019, -0.9973854 , -0.99691116, -0.99639749, -0.9958444 ,
        -0.99525192, -0.99462007, -0.99394887, -0.99323836, -0.99248855,
        -0.99169949, -0.9908712 , -0.99000371, -0.98909706, -0.98815128,
        -0.98716641, -0.9861425 , -0.98507957, -0.98397768, -0.98283687,
        -0.98165717, -0.98043865, -0.97918134, -0.9778853 , -0.97655057,
        -0.97517722, -0.97376528, -0.97231483, -0.97082592, -0.96929861,
        -0.96773295, -0.96612901, -0.96448685, -0.96280654, -0.96108814,
        -0.95933173, -0.95753737, -0.95570513, -0.95383508, -0.95192731]),
 array([2.83026365, 2.83655313, 2.8428426 , 2.84913208, 2.85542155,
        2.86171103, 2.8680005 , 2.87428998, 2.88057945, 2.88686892,
        2.8931584 , 2.89944787, 2.90573735, 2.91202682, 2.9183163 ,
        2.92460577, 2.93089525, 2.93718472, 2.9434742 , 2.94976367,
        2.95605315, 2.96234262, 2.9686321 , 2.97492157, 2.98121105,
        2.98750052, 2.99379   , 3.00007947, 3.00636895, 3.01265842,
        3.0189479 , 3.02523737, 3.03152684, 3.03781632, 3.04410579,
        3.05039527, 3.05668474, 3.06297422, 3.06926369, 3.07555317,
        3.08184264, 3.08813212, 3.09442159, 3.10071107, 3.10700054,
        3.11329002, 3.11957949, 3.12586897, 3.13215844, 3.13844792,
        3.14473739, 3.15102687, 3.15731634, 3.16360582, 3.16989529,
        3.17618476, 3.18247424, 3.18876371, 3.19505319, 3.20134266,
        3.20763214, 3.21392161, 3.22021109, 3.22650056, 3.23279004,
        3.23907951, 3.24536899, 3.25165846, 3.25794794, 3.26423741,
        3.27052689, 3.27681636, 3.28310584, 3.28939531, 3.29568479,
        3.30197426, 3.30826374, 3.31455321, 3.32084268, 3.32713216,
        3.33342163, 3.33971111, 3.34600058, 3.35229006, 3.35857953,
        3.36486901, 3.37115848, 3.37744796, 3.38373743, 3.39002691,
        3.39631638, 3.40260586, 3.40889533, 3.41518481, 3.42147428,
        3.42776376, 3.43405323, 3.44034271, 3.44663218, 3.45292166]))

hvor vi også angiver de tilhørende tids-værdier i t-vektoren.

Skriv Python-kode der finder alle funktionsværdier i f i intervallet \([-0.05, 0.05]\) og angiver de tilhørende t-værdier.

Skriv en Python-funktion der kan findes begge fortegnsskift (nulpunkterne) af f.

Tjek at de fundne nulpunkter stemmer rimeligt overens med de eksakte nulpunkter for \(f\). Jeres metode bør være simpel og også virke på “vektorer” np.array([3, 0.5, 0.1,-0.5, -1, -0.5, 0.5, 2]) hvor vi ikke kender den funktion, der ligger bag funktionsværdierne. Specielt bør jeres metode ikke bruge anden viden om \(\cos\) (fx at den kan differentieres) end den information I har i f-vektoren.

I kan finde inspiration i bisektionsmetoden (som I stødte på i Matematik 1a i efteråret), eller I kan finde inspiration i NumPy-funktioner np.where, np.diff, np.sign eller følgende Python-plot:

plt.figure(figsize=figurstr)
plt.plot(t, np.abs(f), color="r")  # linestyle='-', marker='.'
plt.xlabel("t")
plt.ylabel("|cos(t)|")
plt.show()
../_images/4c105789c82cd9469573c442058a57808c4a72686a2c98f4d8c1be00d2a1260a.png

Solspositionsmodel#

Koordinatsystem#

For at kunne regne på fx solens bevægelse og projektionen af solens stråler ind på panelets normalvektor, skal vi have indført (mindst et) passende koordinatsystem. Vi vælger at placere et fast centrum (origo) \((0,0,0)\) i panelets position. Vi forestiller os herefter tangentplanen til jorden i panelets position. Normalvektoren til dette tangentplan vil være vores \(z\)-akse. \(x\)-aksen og \(y\)-aksen “udspænder” således tangentplanen. Vi fastlægger koordinatsystemet ved at \(x\)-aksen peger mod nord og \(y\)-aksen mod øst. Den opmærksomme læser vil have opdaget at dette er et venstredrejet koordinatsystem. Dette er dog traditionen, og koordinatsystemet kaldes det horisontale koordinatsystem https://en.wikipedia.org/wiki/Horizontal_coordinate_system. Vi har altså nu et referencesolpositionssystem, der fokuserer på observatøren (panelet) på en given bredde- og længdegrad på jordens overflade. Lad os opsummere:

  1. Koordinatsystemet har centrum i solpanelet (kaldet observatøren). Det er traditionelt at bruge et ventredrejet koordinatsystem hvor \(x\)-aksen peger mod nord, \(y\)-aksen peger med øst og \(z\)-aksen peger mod zenit-punktet i retning af normalvektoren til tangentplanen. Koordinatsystemet kaldes det horisontale koordinatsystem, da det er orienteret efter observatørens lokale horisont.

  2. Ethvert objekt kan beskrives i det horisontale koordinatsystem enten i kartesiske koordinater \((x,y,z)\) eller sfæriske koordinater \((r,\theta,\phi)\). Her er \(\phi \in [0,2\pi]\) azimut-vinklen der måles fra \(x\)-aksen mod \(y\)-aksen, \(\theta \in [0,\pi]\) er zenit-vinklen fra \(z\)-aksen og \(r\) er radius.

  3. Det horisontale koordinatsystem er fastgjort til et sted på Jorden, ikke stjernerne eller solen. Objekter på stjernehimlen forestilles placeret på himmelsfæren i en fast afstand og deres placering beskrives ved azimut- og zenit-vinklen (radius ignoreres således i denne type beskrivelse). Over tiden vil zenit- og azimut-vinklen for et objekt på himlen (som fx solen) ændres, da objektet ser ud til at drive hen over himlen med Jordens rotation. Da det horisontale koordinatsystem er defineret af observatørens lokale horisont, vil det samme objekt set fra forskellige steder på Jorden på samme tid have forskellige værdier af azimut- og zenit.

Kilde: https://assessingsolar.org/notebooks/solar_position.html

Vi måler enten vinklerne i radianer \((\theta, \phi) \in [0, \pi] \times [0, 2\pi]\) eller grader \((\theta,\phi) \in [0, 180^\circ] \times [0, 360^\circ]\). Der bruges ofte endnu en vinkel, nemlig solens højdevinkel (\(\alpha\)), som er komplementæren til zenitsolvinklen (\(\alpha = 90^\circ - \theta\)), der altså måler vinklen fra horisont-planen mod \(z\)-aksen.

Skriv en Python-funktion def solar_elevation_angle(theta) der givet \(\theta\) i grader udregner \(\alpha\) i grader.

Hvis vinklerne er givet i radianer, kan vi bruge solar_elevation_angle sammen med np.deg2rad og np.rad2deg. Vi behøver derfor ikke lave funktioner for både grader og radianer, da vi let kan genbruge vores funktioner. Hvis vi fx her ønsker at regne i radianer:

theta_in_rad = np.pi / 3  # zenit-vinkel givet i radianer
print(np.pi / 2. - theta_in_rad)  # højdevinkel givet i radianer
# højdevinkel givet i radianer
# np.deg2rad(solar_elevation_angle(np.rad2deg(theta_in_rad)))  # udkommenter denne

0.5235987755982989

Det anbefales i det videre forløb at regne i radianer, men at angive resultater/vikler i grader hvis dette er mere sigende/beskrivende. Bemærk også at vi ikke behøver at bruge solar_elevation_angle i det følgende da både zenit- og elevationsvinklen vil være tilgængelige.

I det horisontale koordinatsystem er ethvert objekt på himmelsfæren fuldstændig bestemt af zenitvinklen (\(\theta\)) og azimutvinklen \(\phi\). Som nævnt ignorer man den radiale koordinat når alle objekter placeres på himmelsfæren. Vi kan dog sagtens medtage den radiale koordinat:

Antag at solen har en fast afstand \(r_s\) til jorden. Find en rimelig værdi for \(r_s\). Angiv et (matematisk) udtryk for hvordan solens \(xyz\)-koordinat kan udregnes ud fra \(r_s\), \(\theta_s\) og \(\phi_s\), hvor \(\theta_s\) og \(\phi_s\) er hhv. zenit og azimut-vinklen for solens placering.

Der placeres et (flat) solpanel i origo af koordinatsystemet. Enhedsnormalen \(\pmb{u}_p \in \mathbb{R}^3\) til solpanelet har zenit-vinkel \(\theta_p\) og azimut-vinkel \(\phi_p\). Vi betrager en normaliseret (dvs enheds-) solvektor \(\pmb{u}_{s} \in \mathbb{R}^3\) givet ved \((\theta_s, \phi_s)\). Solens vektorfelt er således givet ved \(\pmb{V} = S_0 \pmb{u}_{s}\).

Angiv et (matematisk) udtryk for \(\pmb{u}_p\) og for \(\langle \pmb{u}_{s}, \pmb{u}_p \rangle\) ud fra zenit- og azimut-vinklerne. I bør simplificere udtrykket så det indeholder \(\cos(\phi_p-\phi_s)\) og kun 5 trigonometriske funktioner. Vis at \(-1 \le \langle \pmb{n}_{u}, \pmb{u}_p \rangle \le 1\). Forklar man egne ord hvad det betyder når \(\langle \pmb{u}_{s}, \pmb{u}_p \rangle < 0\).

Skriv en Python-funktion def solar_panel_projection(theta_sol, phi_sol, theta_panel, phi_panel) der returnerer \(\langle \pmb{n}_{s}, \pmb{n}_p \rangle\) når det er positivt og ellers returnerer nul.

Kig igen på jeres Python-funktion def solar_panel_projection(theta_sol, phi_sol, theta_panel, phi_panel). Skriv den om så den virker på NumPy-arrays af zenit- og azimut-vinkler. Du kan teste den på følgende de tre situationer, hvor projektionen bør give 0.707107, 0.0 og 0.0 (eller rettere, med numeriske fejl, bør det give array([7.07106781e-01, 6.12323400e-17, 0.0])). Forklar solpanelets orientering og solens placering i de tre situationer.

theta_sol = np.array([np.pi / 4, np.pi / 2, 0.0])
phi_sol = np.array([np.pi, np.pi / 2, 0.0])
theta_panel = np.array([0.0, np.pi / 2, np.pi])
phi_panel = np.array([np.pi, 0.0, 0.0])

# solar_panel_projection(theta_sol, phi_sol, theta_panel, phi_panel)

Solpositionsmodellering ved Pvlib#

I Python kan solpositionsvinklerne, benævnt \((\theta_s, \phi_s)\), nemt beregnes på ethvert sted ved hjælp af solar positionsalgoritmen (SPA) med pakken pvlib, der som standard er implementeret med National Renewable Energy Laboratory’s SPA-algoritme [Reda og Andreas, 2003, https://www.nrel.gov/docs/fy08osti/34302.pdf]). Vi følger https://assessingsolar.org/notebooks/solar_position.html.

import pandas as pd
import pvlib
from pvlib.location import Location

Vi skal først have defineret observatørens/panelets geografiske placering. Det gøres via objektet pvlib.location.Location i biblioteket pvlib, hvor vi skal angive bl.a. breddegrad, længdegrad, tidszone og højde. Til simulering bruges data for \((\theta_s, \phi_s)\) fra fx kalenderåret 2024, men her i denne indledende øvelse nøjes vi med data for nuværende måned, april 2024:

tidszone = "Europe/Copenhagen"
start_dato = "2024-04-01"
slut_dato = "2024-04-30"
delta_tid = "Min"  # "Min", "H",

# Definition of Location object. Coordinates and elevation of Amager, Copenhagen (Denmark)
site = Location(
    55.660439, 12.604980, tidszone, 10, "Amager (DK)"
)  # latitude, longitude, time_zone, altitude, name

# Definition of a time range of simulation
times = pd.date_range(
    start_dato + " 00:00:00", slut_dato + " 23:59:59", inclusive="left", freq=delta_tid, tz=tidszone
)

Vælg placering/lokation for jeres solpanel, fx DTU. Ret overstående GPS koordinater (målt i DecimalDegrees), højde og navn så det passer med den valgte lokation.

Vi kan nu finde solpositionen ud fra det horisontale koordinatsystem placering i site for det angivne tidsinterval ved følgende kald:

# Estimate Solar Position with the 'Location' object
solpos = site.get_solarposition(times)

# Visualize the resulting DataFrame
solpos.head()
apparent_zenith zenith apparent_elevation elevation azimuth equation_of_time
2024-04-01 00:00:00+02:00 117.855212 117.855212 -27.855212 -27.855212 339.197342 -3.866396
2024-04-01 00:01:00+02:00 117.904726 117.904726 -27.904726 -27.904726 339.473686 -3.866190
2024-04-01 00:02:00+02:00 117.953602 117.953602 -27.953602 -27.953602 339.750301 -3.865984
2024-04-01 00:03:00+02:00 118.001838 118.001838 -28.001838 -28.001838 340.027184 -3.865777
2024-04-01 00:04:00+02:00 118.049434 118.049434 -28.049434 -28.049434 340.304332 -3.865571

Vi ser at DataFramen indeholder solpositionen for hvert minut i april 2024. Tids-samplingen \(\Delta t\) kan styres ved delta_tid = "Min" (minute) sat ovenfor. Når vi senere skal udregne energiproduktionen over hele 2024, kan det være tilstrækkeligt med at kende solpositionen for hver time (for hele året 2024), DataFramen bliver nemlig meget stor, hvis man bruger delta_tid = "Min" for et helt år. Dette klares ved delta_tid = "H" (hour). Bemærk at delta_tid = "M" (month) sætter \(\Delta t\) til en måned (hvilket er for stort til vores behov).

Efter at solvinklerne er estimeret ved hjælp af pvlib, kan de visualiseres, fx for den 1. april:

import matplotlib.dates as mdates

valgt_dato = "2024-04-01"

# Plots for solar zenith and solar azimuth angles
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30, 10))
fig.suptitle("Solar Position Estimation in " + site.name + valgt_dato)

# plot for solar zenith angle
ax1.plot(solpos.loc[valgt_dato].zenith)
ax1.set_ylabel("Solar zenith angle (degree)")
ax1.set_xlabel("Time (hour)")
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H"))

# plot for solar azimuth angle
ax2.plot(solpos.loc[valgt_dato].azimuth)
ax2.set_ylabel("Solar azimuth angle (degree)")
ax2.set_xlabel("Time (hour)")
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H"))
../_images/39122bcf09c3b23e618a506ff5a5166bb9c096399f4c46eb11030d112465ae80.png

Bemærk at \(x\)-aksen bruger UTC-tidszone og der derfor skal \(+2\) for at få dansk tid. De plottede vektorer kan udskrives:

valgt_dato = "2024-04-01"
print(solpos.loc[valgt_dato].zenith)
print(solpos.loc[valgt_dato].elevation)
print(solpos.loc[valgt_dato].azimuth)

2024-04-01 00:00:00+02:00    117.855212
2024-04-01 00:01:00+02:00    117.904726
2024-04-01 00:02:00+02:00    117.953602
2024-04-01 00:03:00+02:00    118.001838
2024-04-01 00:04:00+02:00    118.049434
                                ...    
2024-04-01 23:55:00+02:00    117.237633
2024-04-01 23:56:00+02:00    117.289913
2024-04-01 23:57:00+02:00    117.341566
2024-04-01 23:58:00+02:00    117.392590
2024-04-01 23:59:00+02:00    117.442983
Freq: min, Name: zenith, Length: 1440, dtype: float64
2024-04-01 00:00:00+02:00   -27.855212
2024-04-01 00:01:00+02:00   -27.904726
2024-04-01 00:02:00+02:00   -27.953602
2024-04-01 00:03:00+02:00   -28.001838
2024-04-01 00:04:00+02:00   -28.049434
                               ...    
2024-04-01 23:55:00+02:00   -27.237633
2024-04-01 23:56:00+02:00   -27.289913
2024-04-01 23:57:00+02:00   -27.341566
2024-04-01 23:58:00+02:00   -27.392590
2024-04-01 23:59:00+02:00   -27.442983
Freq: min, Name: elevation, Length: 1440, dtype: float64
2024-04-01 00:00:00+02:00    339.197342
2024-04-01 00:01:00+02:00    339.473686
2024-04-01 00:02:00+02:00    339.750301
2024-04-01 00:03:00+02:00    340.027184
2024-04-01 00:04:00+02:00    340.304332
                                ...    
2024-04-01 23:55:00+02:00    337.993232
2024-04-01 23:56:00+02:00    338.267217
2024-04-01 23:57:00+02:00    338.541481
2024-04-01 23:58:00+02:00    338.816021
2024-04-01 23:59:00+02:00    339.090835
Freq: min, Name: azimuth, Length: 1440, dtype: float64

Her har vi valgt at plotte solvinklerne for 1. april.

Plot solens zenit-, azimut- og elevationsvinkel, dvs. \(\theta_s, \phi_s, \alpha_s\), for hele dagen den 20. april 2024 som funktion af tiden.

Anbefaling: I det følgende har vi brug af at arbejde med vektorer af fx zenit- og azimut-vikler, specielt at kunne finde maksimumsværdier, nulpunkter, integrere, osv. Det anbefales derfor at arbejde med dataen fra solpos som NumPy-arrays. Dette kan gøres ved fx:

np.array(solpos.loc[valgt_dato].elevation)

array([-27.85521243, -27.90472593, -27.95360189, ..., -27.34156572,
       -27.39258993, -27.44298348])

Plot solens elevationsvinkel og find ud af hvornår på dagen solen står højest den 20. april 2024. Forklar hvad det betyder når \(\alpha_s < 0\) eller \(\theta_s > 90^\circ\).

Find tidspunktet for solopgang og solnedgang på DTU den 20. april 2024. Sammenlign med “kendte” værdier fx fra DMI. Hint: Hvis I ønsker præcise værdier skal I bruge apparent_elevation (apparent sun elevation accounting for atmospheric refraction) i stedet for elevation. I behøver ikke tage højde for jordens krumning.

Find solens højeste punkt på himlen (i grader) på sommersolhverv på DTU, og hvornår på dagen det sker? Hint: Du bliver nødt til at ændre på start og slut dato for solpos-objektet.

Lav en Python-funktion som kan beregne solens højeste punkt \(\alpha_{max}\) på himlen (i grader) på en given dato (year-month-day) i en given lokation (fx by) angivet ved en breddegrad og længdegrad. Hint: Svaret bør ikke afhænge af længdegraden, da solens højeste punkt på himlen kun afhænger af breddegraden.

I har i en tidligere opgave fundet et udtryk for solens \(xyz\)-koordinat ud fra \(r_s\), \(\theta_s\) og \(\phi_s\).

Skriv en Python-funktion (til brug med NumPy arrays) der omregner fra solens zenit og azimuth til solens position angivet i \(xyz\)-koordinaten. Husk om I regner i radianer eller grader. Her kan np.deg2rad()-funktionen være nyttig. Det er fint at bruge en cirka værdi for \(r_{s}\) men man kan finde en mere korrekt værdi ved: pvlib.solarposition.nrel_earthsun_distance(times) * 149597870700, hvor 149597870700 er antal meter på en astronomisk enhed AU.

Skriv en Python funktion der omregner fra solens position på himlen i et \(xyz\) koordinater til zenit og azimuth (i grader eller radianer). Her kan np.arctan2(y, x) og np.rad2deg() være nyttige.

Effekt- og energiberegninger#

Anbefaling: Det anbefales at I regner alting i radianer. Husk at man kan bruge np.deg2rad eller np.rad2deg. Hvis I har vinklerne i grader, skal man altså bruge np.rad2deg.

Vi betragter den 20. april. solpos indeholder solpositionsdata for hver minut gennem hele dagen. Da der er 1440 minutter på en dag er solpositionsvinkler over denne dag beskrevet ved en vektor af denne længde, fx:

solpos.loc[valgt_dato].zenith

2024-04-01 00:00:00+02:00    117.855212
2024-04-01 00:01:00+02:00    117.904726
2024-04-01 00:02:00+02:00    117.953602
2024-04-01 00:03:00+02:00    118.001838
2024-04-01 00:04:00+02:00    118.049434
                                ...    
2024-04-01 23:55:00+02:00    117.237633
2024-04-01 23:56:00+02:00    117.289913
2024-04-01 23:57:00+02:00    117.341566
2024-04-01 23:58:00+02:00    117.392590
2024-04-01 23:59:00+02:00    117.442983
Freq: min, Name: zenith, Length: 1440, dtype: float64

Vi betragter kun \(\theta_p \in [0, \pi/2]\), da \(\theta_p > \pi/2\) svarer til at vi begynder at vende panelet med bagsiden op ad.

Lav en Python-funktion som kan udregne fluxen af solens vektorfelt gennem solpanelets flade for hvert minut gennem dagen. I bør bruge solar_panel_projection(theta_sol, phi_sol, theta_panel, phi_panel). Husk under alle omstændigheder kun at medtage sol-zenitvinkler \(\theta_s \in [0, \pi/2]\) (hvorfor?) så panelets flux er nul hvis \(\theta_s\)-værdierne (i en vektor som solpos.loc[valgt_dato].zenith) er over \(\pi/2\) dvs. 90 grader.

For at finde energiproduktionen fra solpanelet skal vi integrere fluxen (dvs effekten) op over den betragtede tidsperiode. Vi regner altid i SI-enheder, men I bør angive endelig resultater (fx den samlede energiproduktion) i relavante enheder (fx både i Joules og kWh). Til at integrere kan vi bruge Trapez-metoden kendt fra Matematik 1b. I må gerne bruge jeres egen implementering af Trapez-metoden, men vi vælger her dog at bruge Simpson’s regel https://da.wikipedia.org/wiki/Simpsons_regel fra SciPy-pakken. Energiproduktionen kan angives per \(m^2\) panel. Husk at medregne solpanelets effektivitet (mht fluxen) som beskrevet i standard-antagelserne.

from scipy import integrate
# flux = np.array(...)  # fra forrige opgave

# husk at tage højde for panelets effektivitet, jvf standard antagelserne

# dx=60 since there are 60 s between time samples

# Udkommenter
# integral_value = integrate.simps(..., dx=60)
# integral_value

flux er her vektoren (np.array) der indeholder fluxen udregnet for hvert minut dagen igennem. Parameteren dx=60 fortæller SciPy at fluxen er samplet hvert minut (1 minut er 60 sekunder i SI-enhed). Hvis I senere vælger at bruge delta_tid = "H" når I skal regne energiproduktionen for et helt år, så skal I huske at fortælle integrate.sipms om den ændrede tids-sampling, nemlig dx=3600. Man kan dog sagtens regne med delta_tid = "Min" hele vejen igennem uden at det bliver for tungt at regne på.

Peg solpanelet mod syd, dvs azimut-vinkel \(\phi_p = 180^\circ\). Udregn energiproduktionen for den 20.april for hver heltals vinkel \(\theta_p\) mellem 0 og 90 grader.

Optimal vinkel#

Vi skal nu endelig betragte energiproduktionen for hele 2024. Kald et nyt solpos-objekter med det relevante tidsinterval.

Peg solpanelet mod syd, dvs azimut-vinkel \(\phi_p = 180^\circ\). Udregn energiproduktionen for hele 2024 for hver heltals vinkel \(\theta_p\) mellem 0 og 90 grader.

Find den optimale vinkel \(\theta_p\) og angiv energiproduktionen. Hvor meget mindre bliver energiproduktionen hvis \(\phi_p\) fx er \(175^\circ\) eller lignende?

Lav en realistisk opsætning af \(X\) antal solpaneler, hvor I vælger \(X\) efter en typisk opsætning på et parcelhus. Solpaneler opsættes efter den optimale vinkel. Udregn energiproduktionen for hver dag og plot dette som funktion af tiden (angivet i dage).

Udvidelser#

I skal nu vælge at arbejde videre med mindst en af følgende udvidelser:

Solpanelsfilm på en krum flade#

En solpanelsfilm (eng: thin solar/power film) er en tynd film som man kan påklistre forskellige bygningsflader og som virker som et solpanel. Tag på ekskursion i DTUs nærområde, find en ikke-plan flade i bymiljøet som egner sig til påmontering af solpanelsfilm. Men en “ikke-plan” flade menes en krum flade hvor fladens normalvektor ikke er konstant. Den kan være taget på et busskur, øverste etage på Jægersborg vandtårn, osv. Det bør være en flade som I kan parametrisere. Opstil en parametrisering for den valgte bygningsflade. Find et datablad for en solpanelsfilm som (med rimelighed) kan bruges på den valgte bygningsflade og angiv de relevante data. Generaliser jeres model og Python-kode så der tages højde for at solpanelets normalvektor ikke er konstant. Udregn den årlige energiproduktion for solpanelsfilmen på bygningsdelen.

Da normalvektoren ikke længere er konstant over hele fladen, skal I nu integrere projektionen af \(\pmb{V}\) ind på fladens enhedsnormelvektor op over hele fladen. Til dette kan I bruge SciPy. Husk at projektionen skal sætte til nul, når den er negativ. I bør vælge en flade hvor skyggen for fladen selv ikke bliver for kompliceret at håndtere.

Medregn skygge fra et træ eller bygning#

I det horizontale koordinatsystem placeres fx en trækrone i en given afstand, fx 10 m. Betragt sfæren \(r_{træ}=10m\) med centrum i Origo (panelets position) og beskriv træet form på denne sfæren. En simpel model er som følger: Antag at der kun er to muligheder: enten \(0\%\) skygge eller \(100 \%\) skygge. Udregn overfladearealet af hele sfæren som funktion af \(r_{træ}\), antag en størrelse af trækronen (fx dens omtrentlige diameter) og udregn hvor man zenit og azimut-grader træet cirka skygger. Dette kunne være at der er \(100 \%\) skygge når \(\theta \in [70,80]\) og \(\phi \in [150,160]\), og energioptaget skal derfor fraregnes når solen befinder sig i dette interval. Træets form bliver noget unaturlig når man bruger akseparalle områder i \((\theta,\phi)\)-planen. Hvilken form svarer det til på sfæren \(r_{træ}=10m\)? Diskuter hvor stor effekt skygger fra træer og bygninger kan have. Er det muligt at beskrive en kuglerund trækrone? Hvordan vil denne form se ud på sfæren \(r_{træ}\)?

I kan alternativt modellere skygge fra en nabobygning. Ideen er den samme men skygge-intervallet for \(\theta\) bør gå helt ned til jorden, dvs \(\theta = 90\). Afhængig af bygningsstørrelse og afstand til solpanel, kan vinkel-intervallerne meget vel blive større, fx \(\theta \in [65,90]\). I kan finde inspiration med dette værktøj: https://www.findmyshadow.com/

Medregn energipriser og forbrugsmønstre#

I stedet for at maksimere den årlige energiproduktion for den valgte solpanelsopsætning, kan det være relevant af minimere parcelhuset årlige energiomkostninger. Find et typisk årligt energiforbrug for et standard parcelhus’, gerne med data time for time. Find tilsvarende tal med timepriserne for energi. Energiforbruget og energipriserne er typisk høje i tidsintervallet 17-20 og højere om vinteren end sommeren. Hvis I ikke kan finde relevante data, så må I antage et forbrug og et priser, fx at priserne er dobbelt så høje i tidsrummet 17-20 end resten af dagen. Find de optimale vinkler \(\theta_p\) og \(\phi_s\), der minimere energiomkostningerne ud fra det angivne forbrug. I bør specielt være opmærksomme på \(\phi_s\), da det er meget muligt at det kan betale sig at dreje panelerne mod vest for at få mest energi ud af aftensolen hvor energipriserne og forbruget er højt.

Priserne time-for-time for 2020 kan findes i: denne fil.

Mere om optimering#

  1. Er sydvendt panel optimalt? Det må det jo være, men lad os undersøge det matematisk. I opgaven ovenfor antog vi \(\phi_p = 180^\circ\). I skal nu droppe denne antagelser og udregne panelets energiproduktion som funktion af både \(\phi_p\) og \(\theta_p\).

  2. Et panel med motor: Antag at solpanelets vinkel kan justeres enten hver dag, hver måned, eller hver kvartal. Hvis I fx vælger hver måned, så skal den optimale vinkel for hver måned angives, og panelet årlige effekt skal derefter udregnes. Sammenlign den årlige effekt med effekten for en tilsvarende fastmonteret panel. Diskuter om det kan betale sig at have solpanelanlæg med paneler hvis vinkel på denne måde kan justeres.

Kilder#

  1. https://www.pveducation.org/

  2. https://www.acs.org/education/resources/highschool/chemmatters/past-issues/archive-2013-2014/how-a-solar-cell-works.html

  3. https://assessingsolar.org/intro.html

  4. https://en.wikipedia.org/wiki/Horizontal_coordinate_system

  5. https://en.wikipedia.org/wiki/Solar_irradiance

Contents