diff --git a/docs/readthedocs/_static/bibliography/bibtexAll.bib b/docs/readthedocs/_static/bibliography/bibtexAll.bib index 23ecc48185..8761d828b4 100644 --- a/docs/readthedocs/_static/bibliography/bibtexAll.bib +++ b/docs/readthedocs/_static/bibliography/bibtexAll.bib @@ -12,12 +12,6 @@ @Article{Maleki.2017 DOI = {10.3390/en10010134} } -@MISC{Itaca_Sun, -author = {Itacanet}, -title = {The Sun As A Source Of Energy}, -howpublished={\url{https://www.itacanet.org/the-sun-as-a-source-of-energy/part-3-calculating-solar-angles/}} -} - @article{Spencer.1971, added-at = {2018-06-18T21:23:34.000+0200}, author = {Spencer, J. W.}, @@ -73,7 +67,7 @@ @book{Quaschning.2013 author = {Quaschning, Volker}, year = {2013}, title = {Regenerative Energiesysteme: Technologie; Berechnung; Simulation}, - url = {http://ebooks.ciando.com/book/index.cfm/bok_id/471802}, + url = {https://www.hanser-elibrary.com/doi/book/10.3139/9783446435711}, price = {39.99 EUR}, address = {M{\"u}nchen}, edition = {8., aktualisierte und erw. Aufl.}, @@ -92,7 +86,8 @@ @Inbook{Schoenberg.1929 pages={1--280}, abstract={Die ersten Versuche, Messungen der Lichtst{\"a}rke der Himmelsk{\"o}rper auszuf{\"u}hren, um auf diesem Wege Aufschl{\"u}sse {\"u}ber ihre Beschaffenheit zu erhalten, fallen in die Zeit jenes gewaltigen Aufschwungs, welchen die Optik durch die Arbeiten von Newton und Huygens am Ende des siebzehnten und zu Beginn des achtzehnten Jahrhunderts erfahren hatte. Schon Huygens1 selbst versuchte die Helligkeit des Sirius mit derjenigen der Sonne zu vergleichen, indem er das Sonnenlicht durch eine kleine {\"O}ffnung im verschlossenen Ende eines langen Rohres abschw{\"a}chte. Der schwedische Physiker Celsius2 suchte nach einem Gesetze der Lichtabnahme beleuchteter Fl{\"a}chen, doch waren seine Schlu{\ss}folgerungen, ebenso wie diejenigen von Huygens, infolge der Unzul{\"a}nglichkeit der angewandten Methoden nicht stichhaltig. Au{\ss}er der Formel der Lichtabnahme mit dem Quadrate der Entfernung besa{\ss} die Physik noch keinerlei R{\"u}stzeug an strengen Definitionen und Gesetzen und keinerlei Apparate zur Messung der Lichtst{\"a}rken. Buffon3 versuchte den Verlust zu bestimmen, den das Sonnenlicht bei Reflexion an gl{\"a}sernen Spiegeln erleidet. Aber erst die gro{\ss} angelegten Arbeiten von Pierre Bouguer4 5 (1698--1758) und Johann Heinrich Lambert (1728--1777) schufen die Grundlagen der Photometrie. In systematischer experimenteller Arbeit, die durch sinnreiche Theorien erl{\"a}utert wird, behandelt Bouguer ihre Hauptprobleme: die Absorption des Lichtes bei der Reflexion und beim Durchgang durch feste und fl{\"u}ssige K{\"o}rper sowie die diffuse Reflexion an matten Fl{\"a}chen und beim Durchgange durch tr{\"u}be Medien.}, isbn={978-3-642-90703-6}, -doi={10.1007/978-3} +doi={10.1007/978-3}, +url = {https://link.springer.com/chapter/10.1007/978-3-642-90703-6_1} } @MISC{WikiAirMass, @@ -122,7 +117,7 @@ @book{Iqbal.1983 author = {Iqbal, Muhammad}, year = {1983}, title = {An Introduction To Solar Radiation}, - url = {http://site.ebrary.com/lib/alltitles/docDetail.action?docID=10678925}, + url = {https://books.google.de/books/about/An_Introduction_To_Solar_Radiation.html?id=3_qWce_mbPsC&redir_esc=y}, address = {Oxford}, publisher = {{Elsevier Science}}, isbn = {9780123737502} diff --git a/docs/readthedocs/models/pv_model.md b/docs/readthedocs/models/pv_model.md index b36df95f00..da64e2067a 100644 --- a/docs/readthedocs/models/pv_model.md +++ b/docs/readthedocs/models/pv_model.md @@ -127,7 +127,7 @@ $$ **References:** * {cite:cts}`Maleki.2017` -* {cite:cts}`Itaca_Sun` +* {cite:cts}`Duffie.2013` p. 17 (formula 1.6.10) ### Solar Altitude Angle @@ -146,7 +146,7 @@ $$ **References:** * {cite:cts}`Maleki.2017` p. 5 -* {cite:cts}`Itaca_Sun` +* {cite:cts}`Duffie.2013` p. 15 (formula 1.6.5) with $\sin (\alpha_s) = \cos (\theta_z)$ ### Zenith Angle @@ -198,11 +198,11 @@ $$ Calculating the air mass ratio by dividing the radius of the earth with approx. effective height of the atmosphere (each in kilometer) $$ -airmassratio = (\frac{6371 km}{9 km}) = 707.8\overline{8} +\mathrm{airmassratio} = (\frac{6371 km}{9 km}) = 707.8\overline{8} $$ $$ -airmass = \sqrt{(707.8\overline{8} \cdot \cos({\theta_z}))^2 +2 \cdot 707.8\overline{8} +1)} - 707.8\overline{8} \cdot \cos{(\theta_z)}) +\mathrm{airmass} = \sqrt{(707.8\overline{8} \cdot \cos({\theta_z}))^2 +2 \cdot 707.8\overline{8} +1)} - 707.8\overline{8} \cdot \cos{(\theta_z)}) $$ **References:** @@ -234,6 +234,8 @@ $$ * {cite:cts}`Zheng.2017` p. 53, formula 2.3b * {cite:cts}`Iqbal.1983` +* {cite:cts}`Spencer.1971` +* {cite:cts}`Duffie.2013` (the fifth ed. seems to have a typo in formula (1.4.1b): factor $0.000719$ is missing a zero) ### Beam Radiation on Sloped Surface @@ -274,7 +276,7 @@ b = (\cos(\phi) \cdot \cos(\delta)) \cdot (\sin(\omega_{2}) - \sin(\omega_{1})) $$ $$ -E_{beam,S} = E_{beam,H} \cdot \frac{a}{b} +E_{\mathrm{beam},S} = E_{\mathrm{beam},H} \cdot \frac{a}{b} $$ **Please note:** $\frac{1}{180}\pi$ is omitted from these formulas, as we are already working with data in *radians*. @@ -286,7 +288,7 @@ $$ **$\omega_1$** = hour angle $\omega$\ **$\omega_2$** = hour angle $\omega$ + 1 hour\ **$\alpha_e$** = surface azimuth angle\ -**$E_{beam,H}$** = beam radiation (horizontal surface) +**$E_{\mathrm{beam},H}$** = beam radiation (horizontal surface) **Reference:** @@ -300,13 +302,15 @@ The diffuse radiation is computed using the Perez model, which divides the radia A cloud index is defined by $$ -\epsilon = \frac{\frac{E_{dif,H} + E_{beam,H}}{E_{dif,H}} + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}{1 + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3} +\epsilon = \frac{\frac{E_{\mathrm{dif},H} + E_{\mathrm{beam},N}}{E_{\mathrm{dif},H}} + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}{1 + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3} $$ +with angle $\theta_z$ values in **degrees** ({cite:cts}`Duffie.2013` p. 94) and $E_{\mathrm{beam},N} = \frac{E_{\mathrm{beam},H}}{\cos (\theta_z)}$ ({cite:cts}`Duffie.2013` p. 95). + Calculating a brightness index $$ -\Delta = m \cdot \frac{E_{dif,H}}{I_{0}} +\Delta = m \cdot \frac{E_{\mathrm{dif},H}}{I_{0}} $$ **Perez Fij coefficients (Myers 2017):** @@ -384,7 +388,7 @@ $$ the diffuse radiation can be calculated: $$ -E_{dif,S} = E_{dif,H} \cdot (\frac{1}{2} \cdot (1 + +E_{\mathrm{dif},S} = E_{\mathrm{dif},H} \cdot (\frac{1}{2} \cdot (1 + cos(\gamma_{e})) \cdot (1- F_{1}) + \frac{a}{b} \cdot F_{1} + F_{2} \cdot \sin(\gamma_{e})) $$ @@ -396,25 +400,27 @@ $$ **$\gamma_{e}$** = slope angle of the surface\ **$I_{0}$** = Extraterrestrial Radiation\ **$m$** = air mass\ -**$E_{beam,H}$** = beam radiation (horizontal surface)\ -**$E_{dif,H}$** = diffuse radiation (horizontal surface) +**$E_{\mathrm{beam},H}$** = beam radiation (horizontal surface)\ +**$E_{\mathrm{beam},N}$** = beam radiation (normal incidence, thus radiation on a plane normal to the direction of the beam)\ +**$E_{\mathrm{dif},H}$** = diffuse radiation (horizontal surface) **References:** * {cite:cts}`Perez.1987` * {cite:cts}`Perez.1990` * {cite:cts}`Myers.2017` p. 96f +* {cite:cts}`Duffie.2013` p. 95f ### Reflected Radiation on Sloped Surface $$ -E_{ref,S} = E_{Ges,H} \cdot \frac{\rho}{2} \cdot (1- +E_{\mathrm{ref},S} = E_{\mathrm{Ges},H} \cdot \frac{\rho}{2} \cdot (1- \cos(\gamma_{e})) $$ *with*\ -**$E_{Ges,H}$** = total horizontal radiation ($E_{beam,H} + E_{dif,H})$\ +**$E_{\mathrm{Ges},H}$** = total horizontal radiation ($E_{\mathrm{beam},H} + E_{\mathrm{dif},H})$\ **$\gamma_e$** = slope angle of the surface\ **$\rho$** = albedo @@ -428,13 +434,13 @@ $$ Received energy is calculated as the sum of all three types of irradiation. $$ -E_{total} = E_{beam,S} + E_{dif,S} + E_{ref,S} +E_{\mathrm{total}} = E_{\mathrm{beam},S} + E_{\mathrm{dif},S} + E_{\mathrm{ref},S} $$ *with*\ -**$E_{beam,S}$** = Beam radiation\ -**$E_{dif,S}$** = Diffuse radiation\ -**$E_{ref,S}$** = Reflected radiation +**$E_{\mathrm{beam},S}$** = Beam radiation\ +**$E_{\mathrm{dif},S}$** = Diffuse radiation\ +**$E_{\mathrm{ref},S}$** = Reflected radiation A generator correction factor (depending on month surface slope $\gamma_{e}$) and a temperature correction factor (depending on month) multiplied on top. diff --git a/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala b/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala index 5d83478da0..98f2162ed2 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala @@ -24,7 +24,6 @@ import tech.units.indriya.unit.Units._ import java.time.ZonedDateTime import java.util.UUID -import java.util.stream.IntStream import scala.math._ final case class PvModel private ( @@ -545,20 +544,21 @@ final case class PvModel private ( val gammaEInRad = gammaE.toRadians // == brightness index beta ==// - val beta = eDifH * airMass / extraterrestrialRadiationI0 + val delta = eDifH * airMass / extraterrestrialRadiationI0 // == cloud index epsilon ==// - // if we have no clouds, the epsilon bin is 8, as epsilon bin for an epsilon in [6.2, inf.[ = 8 - var x = 8 + val x = if (eDifH.value.doubleValue > 0) { + // if we have diffuse radiation on horizontal surface we have to consider + // the clearness parameter epsilon, which then gives us an epsilon bin x - if (eDifH.value.doubleValue > 0) { - // if we have diffuse radiation on horizontal surface we have to check if we have another epsilon due to clouds get the epsilon - var epsilon = ((eDifH + eBeamH) / eDifH + + // Beam radiation is required on a plane normal to the beam direction (normal incidence), + // thus dividing by cos theta_z + var epsilon = ((eDifH + eBeamH / cos(thetaZInRad)) / eDifH + (5.535d * 1.0e-6) * pow( - thetaZInRad, + thetaZ.toDegrees, 3, )) / (1d + (5.535d * 1.0e-6) * pow( - thetaZInRad, + thetaZ.toDegrees, 3, )) @@ -579,18 +579,23 @@ final case class PvModel private ( // get the corresponding bin val finalEpsilon = epsilon - x = IntStream - .range(0, discreteSkyClearnessCategories.length) - .filter((i: Int) => - (finalEpsilon - discreteSkyClearnessCategories(i)( - 0 - ) >= 0) && (finalEpsilon - discreteSkyClearnessCategories( - i - )(1) < 0) - ) - .findFirst - .getAsInt + 1 + discreteSkyClearnessCategories.indices + .find { i => + (finalEpsilon - + discreteSkyClearnessCategories(i)(0) >= 0) && + (finalEpsilon - + discreteSkyClearnessCategories(i)(1) < 0) + } + .map(_ + 1) + .getOrElse(8) + } else { + // epsilon in [6.2, inf.[ + 8 } + } else { + // if we have no clouds, the epsilon bin is 8, + // as the epsilon bin for an epsilon in [6.2, inf.[ is 8 + 8 } // calculate the f_ij components based on the epsilon bin @@ -604,19 +609,17 @@ final case class PvModel private ( val f22 = 0.0012 * pow(x, 3) - 0.0067 * pow(x, 2) + 0.0091 * x - 0.0269 val f23 = 0.0052 * pow(x, 3) - 0.0971 * pow(x, 2) + 0.2856 * x - 0.1389 - // calculate circuumsolar brightness coefficient f1 and horizon brightness coefficient f2 - val f1 = max(0, f11 + f12 * beta + f13 * thetaZInRad) - val f2 = f21 + f22 * beta + f23 * thetaZInRad + // calculate circumsolar brightness coefficient f1 and horizon brightness coefficient f2 + val f1 = max(0, f11 + f12 * delta + f13 * thetaZInRad) + val f2 = f21 + f22 * delta + f23 * thetaZInRad val aPerez = max(0, cos(thetaGInRad)) val bPerez = max(cos(1.4835298641951802), cos(thetaZInRad)) // finally calculate the diffuse radiation on an inclined surface eDifH * ( - ((1 + cos( - gammaEInRad - )) / 2) * (1 - f1) + (f1 * (aPerez / bPerez)) + (f2 * sin( - gammaEInRad - )) + ((1 + cos(gammaEInRad)) / 2) * (1 - f1) + + (f1 * (aPerez / bPerez)) + + (f2 * sin(gammaEInRad)) ) } diff --git a/src/test/scala/edu/ie3/simona/agent/em/EmAgentIT.scala b/src/test/scala/edu/ie3/simona/agent/em/EmAgentIT.scala index 8667aaa84e..630f17575d 100644 --- a/src/test/scala/edu/ie3/simona/agent/em/EmAgentIT.scala +++ b/src/test/scala/edu/ie3/simona/agent/em/EmAgentIT.scala @@ -259,11 +259,11 @@ class EmAgentIT scheduler.expectMessage(Completion(storageAgent)) /* TICK 0 - LOAD: 0.000269 MW - PV: -0.005685 MW + LOAD: 0.269 kW + PV: -5.617 kW STORAGE: SOC 0 % -> charge with 5 kW - -> remaining -0.0004161 MW + -> remaining -0.348 kW */ emAgentActivation ! Activation(0) @@ -272,7 +272,7 @@ class EmAgentIT 0, weatherService.ref.toClassic, WeatherData( - WattsPerSquareMeter(400d), + WattsPerSquareMeter(540d), WattsPerSquareMeter(200d), Celsius(0d), MetersPerSecond(0d), @@ -285,19 +285,19 @@ class EmAgentIT emResult.getInputModel shouldBe emInput.getUuid emResult.getTime shouldBe 0L.toDateTime emResult.getP should equalWithTolerance( - -0.000416087825.asMegaWatt + -0.00034885012.asMegaWatt ) - emResult.getQ should equalWithTolerance(0.0000882855367.asMegaVar) + emResult.getQ should equalWithTolerance(0.00008828554.asMegaVar) } scheduler.expectMessage(Completion(emAgentActivation, Some(7200))) /* TICK 7200 - LOAD: 0.000269 MW (unchanged) - PV: -0.003797 MW + LOAD: 0.269 kW (unchanged) + PV: -3.651 kW STORAGE: SOC 63.3 % - -> charge with 3.5282 kW - -> remaining 0 MW + -> charge with 3.382 kW + -> remaining 0 kW */ emAgentActivation ! Activation(7200) @@ -322,24 +322,24 @@ class EmAgentIT emResult.getQ should equalWithTolerance(0.0000882855367.asMegaVar) } - scheduler.expectMessage(Completion(emAgentActivation, Some(13107))) + scheduler.expectMessage(Completion(emAgentActivation, Some(13362))) - /* TICK 13107 - LOAD: 0.000269 MW (unchanged) - PV: -0.003797 MW (unchanged) + /* TICK 13362 + LOAD: 0.269 kW (unchanged) + PV: -3.651 kW (unchanged) STORAGE: SOC 100 % -> charge with 0 kW - -> remaining -0.003528 MW + -> remaining -3.382 kW */ - emAgentActivation ! Activation(13107) + emAgentActivation ! Activation(13362) resultListener.expectMessageType[ParticipantResultEvent] match { case ParticipantResultEvent(emResult: EmResult) => emResult.getInputModel shouldBe emInput.getUuid - emResult.getTime shouldBe 13107L.toDateTime + emResult.getTime shouldBe 13362L.toDateTime emResult.getP should equalWithTolerance( - -0.0035281545552.asMegaWatt + -0.003382375474.asMegaWatt ) emResult.getQ should equalWithTolerance(0.0000882855367.asMegaVar) } @@ -347,11 +347,11 @@ class EmAgentIT scheduler.expectMessage(Completion(emAgentActivation, Some(14400))) /* TICK 14400 - LOAD: 0.000269 MW (unchanged) - PV: -0.000066 MW + LOAD: 0.269 kW (unchanged) + PV: -0.066 kW STORAGE: SOC 100 % -> charge with -0.202956 kW - -> remaining 0 MW + -> remaining 0 kW */ // send weather data before activation, which can happen @@ -558,11 +558,11 @@ class EmAgentIT val weatherDependentAgents = Seq(pvAgent, heatPumpAgent) /* TICK 0 - LOAD: 0.000269 MW - PV: -0.005685 MW + LOAD: 0.269 kW + PV: -5.617 kW Heat pump: off, can be turned on or stay off -> set point ~3.5 kW (bigger than 50 % rated apparent power): turned on - -> remaining -0.000566 MW + -> remaining -0.499 kW */ emAgentActivation ! Activation(0) @@ -572,7 +572,7 @@ class EmAgentIT 0, weatherService.ref.toClassic, WeatherData( - WattsPerSquareMeter(400d), + WattsPerSquareMeter(540d), WattsPerSquareMeter(200d), Celsius(0d), MetersPerSecond(0d), @@ -586,7 +586,7 @@ class EmAgentIT emResult.getInputModel shouldBe emInput.getUuid emResult.getTime shouldBe 0.toDateTime emResult.getP should equalWithTolerance( - -0.000566087824.asMegaWatt + -0.000498850118.asMegaWatt ) emResult.getQ should equalWithTolerance(0.001073120041.asMegaVar) } @@ -594,11 +594,11 @@ class EmAgentIT scheduler.expectMessage(Completion(emAgentActivation, Some(7200))) /* TICK 7200 - LOAD: 0.000269 MW (unchanged) - PV: -0.003797 MW + LOAD: 0.269 kW (unchanged) + PV: -3.651 kW Heat pump: running (turned on from last request), can also be turned off -> set point ~3.5 kW (bigger than 50 % rated apparent power): stays turned on with unchanged state - -> remaining 0 MW + -> remaining 1.468 kW */ emAgentActivation ! Activation(7200) @@ -621,17 +621,18 @@ class EmAgentIT case ParticipantResultEvent(emResult: EmResult) => emResult.getInputModel shouldBe emInput.getUuid emResult.getTime shouldBe 7200.toDateTime - emResult.getP should equalWithTolerance(0.00132184544484.asMegaWatt) + emResult.getP should equalWithTolerance(0.001467624526.asMegaWatt) emResult.getQ should equalWithTolerance(0.001073120041.asMegaVar) } scheduler.expectMessage(Completion(emAgentActivation, Some(14400))) /* TICK 14400 - LOAD: 0.000269 MW (unchanged) - PV: -0.000066 MW + LOAD: 0.269 kW (unchanged) + PV: -0.066 kW Heat pump: Is still running, can still be turned off -> flex signal is 0 MW: Heat pump is turned off + -> remaining 0.203 kW */ emAgentActivation ! Activation(14400) @@ -662,10 +663,11 @@ class EmAgentIT scheduler.expectMessage(Completion(emAgentActivation, Some(21600))) /* TICK 21600 - LOAD: 0.000269 MW (unchanged) - PV: -0.000032 MW + LOAD: 0.269 kW (unchanged) + PV: -0.026 kW Heat pump: Is not running, can run or stay off -> flex signal is 0 MW: Heat pump is turned off + -> remaining 0.242 kW */ emAgentActivation ! Activation(21600) @@ -675,6 +677,7 @@ class EmAgentIT 21600, weatherService.ref.toClassic, WeatherData( + // Same irradiation, but different angle of the sun WattsPerSquareMeter(5d), WattsPerSquareMeter(5d), Celsius(0d), @@ -688,17 +691,18 @@ class EmAgentIT case ParticipantResultEvent(emResult: EmResult) => emResult.getInputModel shouldBe emInput.getUuid emResult.getTime shouldBe 21600.toDateTime - emResult.getP should equalWithTolerance(0.0002367679996.asMegaWatt) + emResult.getP should equalWithTolerance(0.000242284024.asMegaWatt) emResult.getQ should equalWithTolerance(0.000088285537.asMegaVar) } scheduler.expectMessage(Completion(emAgentActivation, Some(28665))) /* TICK 28666 - LOAD: 0.000269 MW (unchanged) - PV: -0.000032 MW (unchanged) + LOAD: 0.269 kW (unchanged) + PV: -0.026 kW (unchanged) Heat pump: Is turned on again and cannot be turned off -> flex signal is no control -> 0.00485 MW + -> remaining 5.092 kW */ emAgentActivation ! Activation(28665) @@ -707,7 +711,7 @@ class EmAgentIT case ParticipantResultEvent(emResult: EmResult) => emResult.getInputModel shouldBe emInput.getUuid emResult.getTime shouldBe 28665.toDateTime - emResult.getP should equalWithTolerance(0.0050867679996.asMegaWatt) + emResult.getP should equalWithTolerance(0.005092284024.asMegaWatt) emResult.getQ should equalWithTolerance(0.001073120040.asMegaVar) } diff --git a/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala b/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala index 157ff0a041..c7aef19b1c 100644 --- a/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala +++ b/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala @@ -13,18 +13,11 @@ import edu.ie3.datamodel.models.input.{NodeInput, OperatorInput} import edu.ie3.datamodel.models.voltagelevels.GermanVoltageLevelUtils import edu.ie3.simona.test.common.{DefaultTestData, UnitSpec} import edu.ie3.util.quantities.PowerSystemUnits._ -import edu.ie3.util.scala.quantities.{ - ApparentPower, - Irradiation, - Kilovoltamperes, - Megavars, - ReactivePower, - WattHoursPerSquareMeter, -} +import edu.ie3.util.scala.quantities._ import org.locationtech.jts.geom.{Coordinate, GeometryFactory, Point} import org.scalatest.GivenWhenThen import squants.Each -import squants.energy.{Kilowatts, Power} +import squants.energy.Kilowatts import squants.space.{Angle, Degrees, Radians} import tech.units.indriya.quantity.Quantities.getQuantity import tech.units.indriya.unit.Units._ @@ -97,7 +90,6 @@ class PvModelSpec extends UnitSpec with GivenWhenThen with DefaultTestData { private implicit val apparentPowerTolerance: ApparentPower = Kilovoltamperes( 1e-10 ) - private implicit val powerTolerance: Power = Kilowatts(1e-10) private implicit val reactivePowerTolerance: ReactivePower = Megavars(1e-10) "A PV Model" should { @@ -763,23 +755,37 @@ class PvModelSpec extends UnitSpec with GivenWhenThen with DefaultTestData { } } - "calculate the estimate diffuse radiation eDifS correctly" in { + "calculate the estimated diffuse radiation eDifS correctly" in { + def megaJoule2WattHours(megajoule: Double): Double = { + megajoule / (3.6 / 1000) + } + val testCases = Table( ("thetaGDeg", "thetaZDeg", "gammaEDeg", "airMass", "I0", "eDifSSol"), - (37.0d, 62.2d, 60d, 2.13873080095658d, 1399.0077631849722d, - 216.46615469650982d), + // Reference Duffie 4th ed., p.95 + // I_0 = 5.025 MJ * 277.778 Wh/MJ = 1395.83445 Wh + // eDifSSol is 0.79607 MJ (0.444 + 0.348 + 0.003) if one only calculates the relevant terms + // from I_T on p. 96, but Duffie uses fixed f values, so the inaccuracy is fine (approx. 4.5 Wh/m^2 or 0.016 MJ/m^2) + ( + 37.0d, + 62.2d, + 60d, + 2.144d, + megaJoule2WattHours(5.025), + megaJoule2WattHours(0.812140993078191252), + ), ) forAll(testCases) { (thetaGDeg, thetaZDeg, gammaEDeg, airMass, I0, eDifSSol) => - // Reference p.95 - // https://www.sku.ac.ir/Datafiles/BookLibrary/45/John%20A.%20Duffie,%20William%20A.%20Beckman(auth.)-Solar%20Engineering%20of%20Thermal%20Processes,%20Fourth%20Edition%20(2013).pdf + // Reference Duffie 4th ed., p.95 Given("using the input data") // Beam Radiation on horizontal surface val eBeamH = - 67.777778d // 1 MJ/m^2 = 277,778 Wh/m^2 -> 0.244 MJ/m^2 = 67.777778 Wh/m^2 + megaJoule2WattHours(0.244) // 0.244 MJ/m^2 = 67.777778 Wh/m^2 // Diffuse Radiation on a horizontal surface - val eDifH = 213.61111d // 0.769 MJ/m^2 = 213,61111 Wh/m^2 + val eDifH = + megaJoule2WattHours(0.796) // 0.796 MJ/m^2 = 221.111111 Wh/m^2 When("the diffuse radiation is calculated") val eDifSCalc = pvModel.calcDiffuseRadiationOnSlopedSurfacePerez(