Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rewrite functions in pv_model #689

Draft
wants to merge 22 commits into
base: dev
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
dbfc369
initial implementation of new functions
SimonHuette Dec 20, 2023
55f24cd
Merge branch 'dev' into sh/#688-rewrite-pv-model
SimonHuette Dec 20, 2023
66cd6eb
noticed an error in Test
SimonHuette Dec 20, 2023
e7473d0
Merge remote-tracking branch 'origin/sh/#688-rewrite-pv-model' into s…
SimonHuette Dec 20, 2023
1c92a27
attempt to fix error with wrong unit in the variable timeFrame
SimonHuette Jan 7, 2024
684e9cd
added comment
SimonHuette Jan 17, 2024
eb59fd4
fixed epsilon + testing beam radiation
SimonHuette Jan 22, 2024
bb05389
testing beam and diffuse radiation
SimonHuette Jan 25, 2024
e7e9bcc
further testing
SimonHuette Jan 26, 2024
bd362d9
testing
SimonHuette Jan 29, 2024
8fbf117
noticed error in CalcExtraterrestrialRadiation
SimonHuette Jan 29, 2024
993b0fd
adjusted values for extraterrestrial radiation in test function
SimonHuette Jan 29, 2024
8e8549e
found error in GroovyTest (wrong data)
SimonHuette Feb 7, 2024
88d1062
Merge branch 'dev' into sh/#688-rewrite-pv-model
SimonHuette Mar 5, 2024
761e8eb
spotless apply
SimonHuette Mar 5, 2024
e725f06
added sources in documentation
SimonHuette Mar 30, 2024
b755427
angle in degrees documentation
SimonHuette Mar 30, 2024
ac0828e
updated links in library
SimonHuette Mar 30, 2024
69536ea
A tiny bit of cleanup
sebastian-peter Sep 16, 2024
b5c03d1
Reverting different sloped surface workaround
sebastian-peter Sep 16, 2024
4a54db3
Reverting changed eccentricity correction factor
sebastian-peter Sep 16, 2024
27e382a
Reverting changes in test
sebastian-peter Sep 16, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 61 additions & 17 deletions src/main/scala/edu/ie3/simona/model/participant/PvModel.scala
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,7 @@ final case class PvModel private (
lat,
gammaE,
alphaE,
duration,
)

// === Diffuse Radiation Parameters ===//
Expand Down Expand Up @@ -337,7 +338,7 @@ final case class PvModel private (
val e0 = 1.000110 +
0.034221 * cos(jInRad) +
0.001280 * sin(jInRad) +
0.000719 * cos(2d * jInRad) +
0.00719 * cos(2d * jInRad) +
Copy link
Member

@sebastian-peter sebastian-peter Mar 5, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to Iqbal, this is 0.000719 with three zeroes. The primary source, Spencer, also states the same. Duffie and Zheng as well. Please check again.

Furthermore, could you add the missing sources I just mentioned to the list of sources in documentation as well?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

image

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is Duffie p. 9, where I noticed the 2 zeros only. Can you maybe send me the PDFs of the other books, so I can look it up aswell?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interestingly, the fifth edition of Duffie seems to differ from all earlier editions in this point. Given that the cited sources stayed the same and I didn't find any other explanation (and the layout seems broken in this edition for the first time), I assume that the change did not happen on purpose and the number with three zeros is still correct. There should be some more investigations here, though.

Copy link
Member

@sebastian-peter sebastian-peter Sep 16, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It'd be great if we added commentary to the readthedocs documentation regarding this issue, i.e. why we're using the fourth edition and not the fifth. (Extraterrestrial Radiation)

0.000077 * sin(2d * jInRad)

// solar constant in W/m2
Expand Down Expand Up @@ -409,32 +410,27 @@ final case class PvModel private (
omegaSS: Angle,
omegaSR: Angle,
): Option[(Angle, Angle)] = {
val thetaGInRad = thetaG.toRadians
val omegaSSInRad = omegaSS.toRadians
val omegaSRInRad = omegaSR.toRadians

val omegaOneHour = toRadians(15d)
val omegaHalfHour = omegaOneHour / 2d

val omega1InRad = omega.toRadians // requested hour
val omega2InRad = omega1InRad + omegaOneHour // requested hour plus 1 hour

// (thetaG < 90°): sun is visible
// (thetaG > 90°), otherwise: sun is behind the surface -> no direct radiation
if (
thetaGInRad < toRadians(90)
// omega1 and omega2: sun has risen and has not set yet
&& omega2InRad > omegaSRInRad + omegaHalfHour
&& omega1InRad < omegaSSInRad - omegaHalfHour
// requested time is between sunrise and sunset (+/- one hour)
omega1InRad > omegaSRInRad - omegaOneHour
&& omega1InRad < omegaSSInRad
) {

val (finalOmega1, finalOmega2) =
if (omega1InRad < omegaSRInRad) {
// requested time earlier than sunrise
(omegaSRInRad, omegaSRInRad + omegaOneHour)
} else if (omega2InRad > omegaSSInRad) {
// sunset earlier than requested time
(omegaSSInRad - omegaOneHour, omegaSSInRad)
(omegaSRInRad, omega2InRad)
} else if (omega1InRad > omegaSSInRad - omegaOneHour) {
// requested time is less than one hour before sunset
(omega1InRad, omegaSSInRad)
} else {
(omega1InRad, omega2InRad)
}
Expand Down Expand Up @@ -464,13 +460,23 @@ final case class PvModel private (
* @return
* the beam radiation on the sloped surface
*/

def calculateTimeFrame(omegas: Option[(Angle, Angle)]): Double = {
omegas match {
case Some((omega1, omega2)) =>
(omega2 - omega1).toDegrees / 15

case None => 0d
}
}
private def calcBeamRadiationOnSlopedSurface(
eBeamH: Irradiation,
omegas: Option[(Angle, Angle)],
delta: Angle,
latitude: Angle,
gammaE: Angle,
alphaE: Angle,
duration: Time,
): Irradiation = {

omegas match {
Expand All @@ -482,6 +488,9 @@ final case class PvModel private (

val omega1InRad = omega1.toRadians
val omega2InRad = omega2.toRadians
// variable that accounts for cases when the integration interval is shorter than 15° (1 hour equivalent), when the time is close to sunrise or sunset
val timeFrame =
(omega2 - omega1).toDegrees / 15d / duration.toHours // original term: (omega2 - omega1).toRadians * 180 / Math.PI / 15, since a one hour difference equals 15°

val a = ((sin(deltaInRad) * sin(latInRad) * cos(gammaEInRad)
- sin(deltaInRad) * cos(latInRad) * sin(gammaEInRad) * cos(
Expand All @@ -503,7 +512,7 @@ final case class PvModel private (

// in rare cases (close to sunrise) r can become negative (although very small)
val r = max(a / b, 0d)
eBeamH * r
eBeamH * r * timeFrame
case None => WattHoursPerSquareMeter(0d)
}
}
Expand Down Expand Up @@ -533,6 +542,41 @@ final case class PvModel private (
* @return
* the diffuse radiation on the sloped surface
*/

private def calcEpsilon(
eDifH: Irradiation,
eBeamH: Irradiation,
thetaZ: Angle,
): Double = {
val thetaZInRad = thetaZ.toRadians

((eDifH + eBeamH / cos(thetaZInRad)) / eDifH + (5.535d * 1.0e-6) * pow(
thetaZ.toDegrees,
3,
)) /
(1d + (5.535d * 1.0e-6) * pow(thetaZ.toDegrees, 3))
}

private def calcEpsilonOld(
eDifH: Irradiation,
eBeamH: Irradiation,
thetaZ: Angle,
): Double = {
val thetaZInRad = thetaZ.toRadians

((eDifH + eBeamH) / eDifH + (5.535d * 1.0e-6) * pow(thetaZ.toRadians, 3)) /
(1d + (5.535d * 1.0e-6) * pow(thetaZ.toRadians, 3))
}

private def firstFraction(
eDifH: Irradiation,
eBeamH: Irradiation,
thetaZ: Angle,
): Double = {

(eDifH + eBeamH) / eDifH
}

private def calcDiffuseRadiationOnSlopedSurfacePerez(
eDifH: Irradiation,
eBeamH: Irradiation,
Expand All @@ -555,12 +599,12 @@ final case class PvModel private (

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 +
var epsilon = ((eDifH + eBeamH / cos(thetaZInRad)) / eDifH +
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both Duffie and the original Perez use the formula as before. Where did you find the formula with cos thetaZ again?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

image

Copy link
Contributor Author

@SimonHuette SimonHuette Mar 30, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is my reference (Duffie p. 95f). Note that "I_b,n" is not the beam radiation on a horizontal surface, but rather the beam radiation on a surface normal to the direction of the sun. This can be calculated by dividing the beam radiation on a horizontal surface "I_b" (=eBeamH in code) by cos(ThetaZ). This is also demonstrated in Duffie: example 2.16.2

Copy link
Member

@sebastian-peter sebastian-peter Sep 16, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we need some more investigation here, probably focused on the weather data specification. I_b,n is certainly the "normal beam incidence radiation", which means the radiation in direction of the beam, i.e. on a surface that is exactly facing the sun, as you explained. I_b is certainly the radiationon a horizontal plane (i.e. on earth's surface).
All of this is illustrated by figure 1.8.1 in Duffie (4th ed).
The current source we use, Myers, is not quite clear on this. In Perez 1990 it is clearly stated, though: The "normal incidence direct irradiance" is asked for (formula (1), p. 273).

The question now becomes which type of radiation the weather data sources provide. For ERA5, solar radiation always seems to be specified on a horizontal plane. I have not found the specs for ICON and COSMO yet.

So it seems likely to me that you're right with your assumptions. What should be done now is adapting the formula in documentation (replace E_beamH with E_beamH/cos(theta_Z)). Also it be nice to add a hint for fig. 1.8.1 in Duffie to the source.

(5.535d * 1.0e-6) * pow(
thetaZInRad,
thetaZ.toDegrees,
SimonHuette marked this conversation as resolved.
Show resolved Hide resolved
3,
)) / (1d + (5.535d * 1.0e-6) * pow(
thetaZInRad,
thetaZ.toDegrees,
3,
))

Expand Down
40 changes: 24 additions & 16 deletions src/test/groovy/edu/ie3/simona/model/participant/PvModelTest.groovy
Original file line number Diff line number Diff line change
Expand Up @@ -300,9 +300,9 @@ class PvModelTest extends Specification {

where:
j || I0Sol
0d || 1414.91335d // Jan 1st
2.943629280897834d || 1322.494291080537598d // Jun 21st
4.52733626243351d || 1355.944773587800003d // Sep 21st
0d || 1423.7592070000003d // Jan 1st
2.943629280897834d || 1330.655828592125d // Jun 21st
4.52733626243351d || 1347.6978765254157d // Sep 21st
}

def "Calculate the angle of incidence thetaG"() {
Expand All @@ -322,7 +322,7 @@ class PvModelTest extends Specification {
Angle omegaRad = Sq.create(Math.toRadians(omegaDeg), Radians$.MODULE$)
//Inclination Angle of the surface
Angle gammaERad = Sq.create(Math.toRadians(gammaEDeg), Radians$.MODULE$)
//Sun's azimuth
//Surface azimuth
Angle alphaERad = Sq.create(Math.toRadians(alphaEDeg), Radians$.MODULE$)

when:
Expand Down Expand Up @@ -449,19 +449,22 @@ class PvModelTest extends Specification {

expect:
"- should calculate the beam contribution,"

pvModel.calcBeamRadiationOnSlopedSurface(eBeamH, omegas, delta, latitudeInRad, gammaE, alphaE) =~ Sq.create(eBeamSSol, WattHoursPerSquareMeter$.MODULE$)
def calculatedsunsetangle = pvModel.calcSunsetAngleOmegaSS(latitudeInRad, delta)
def calculateAngleDifference = (omegas.get()._1() - omegas.get()._2())
def timeframe = pvModel.calculateTimeFrame(omegas)
def beamradiation = pvModel.calcBeamRadiationOnSlopedSurface(eBeamH, omegas, delta, latitudeInRad, gammaE, alphaE)
beamradiation =~ Sq.create(eBeamSSol, WattHoursPerSquareMeter$.MODULE$)


where: "the following parameters are given"
latitudeInDeg | slope | azimuth | deltaIn | omegaIn | thetaGIn || eBeamSSol
40d | 0d | 0d | -11.6d | -37.5d | 37.0d || 67.777778d // flat surface => eBeamS = eBeamH
40d | 60d | 0d | -11.6d | -37.5d | 37.0d || 112.84217113154841369d // 2011-02-20T09:00:00
40d | 60d | 0d | -11.6d | -78.0d | 75.0d || 210.97937494450755d // sunrise
40d | 60d | 0d | -11.6d | 62.0d | 76.0d || 199.16566536224116d // sunset
//40d | 0d | 0d | -11.6d | -37.5d | 37.0d || 67.777778d // flat surface => eBeamS = eBeamH
//40d | 60d | 0d | -11.6d | -37.5d | 37.0d || 112.84217113154841369d // 2011-02-20T09:00:00
//40d | 60d | 0d | -11.6d | -78.0d | 75.0d || 210.97937494450755d // sunrise
//40d | 60d | 0d | -11.6d | 62.0d | 76.0d || 199.16566536224116d // sunset
40d | 60d | 0d | -11.6d | 69.0d | 89.9d || 245.77637766673405d // sunset, cut off
40d | 60d | 0d | -11.6d | 75.0d | 89.9d || 0d // no sun
40d | 60d | -90.0d | -11.6d | 60.0d | 91.0d || 0d // no direct beam
//40d | 60d | 0d | -11.6d | 75.0d | 89.9d || 0d // no sun
//40d | 60d | -90.0d | -11.6d | 60.0d | 91.0d || 0d // no direct beam
}

def "Calculate the estimate diffuse radiation eDifS"() {
Expand All @@ -473,9 +476,9 @@ class PvModelTest extends Specification {
// 0.244 MJ/m^2 = 67.777778 Wh/m^2
//Beam Radiation on horizontal surface
Irradiation eBeamH = Sq.create(67.777778d, WattHoursPerSquareMeter$.MODULE$)
// 0.769 MJ/m^2 = 213,61111 Wh/m^2
// 0.796 MJ/m^2 = 221,111288 Wh/m^2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a typo we made at some point, so this correction should find its way into the new PvModelSpec as well.

//Diffuse beam Radiation on horizontal surface
Irradiation eDifH = Sq.create(213.61111d, WattHoursPerSquareMeter$.MODULE$)
Irradiation eDifH = Sq.create(221.111288d, WattHoursPerSquareMeter$.MODULE$)
//Incidence Angle
Angle thetaG = Sq.create(Math.toRadians(thetaGIn), Radians$.MODULE$)
//Zenith Angle
Expand All @@ -487,11 +490,16 @@ class PvModelTest extends Specification {
"- should calculate the beam diffusion"
// == 0,7792781569074354 MJ/m^2

pvModel.calcDiffuseRadiationOnSlopedSurfacePerez(eDifH, eBeamH, airMass, I0Quantity, thetaZ, thetaG, gammaE) =~ Sq.create(eDifSSol, WattHoursPerSquareMeter$.MODULE$)
def epsilon = pvModel.calcEpsilon(eDifH, eBeamH, thetaZ) // epsilon(Duffie) = 1,28451252
def epsilonOld = pvModel.calcEpsilonOld(eDifH, eBeamH, thetaZ)
def firstFraction = pvModel.firstFraction(eDifH, eBeamH, thetaZ)

def diffuseradiation = pvModel.calcDiffuseRadiationOnSlopedSurfacePerez(eDifH, eBeamH, airMass, I0Quantity, thetaZ, thetaG, gammaE)
diffuseradiation =~ Sq.create(eDifSSol, WattHoursPerSquareMeter$.MODULE$)

where: "the following parameters are given"
thetaGIn | thetaZIn | slope | airMass | I0 || eDifSSol
37.0 | 62.2 | 60 | 2.13873080095658d | 1399.0077631849722d || 216.46615469650982d
37.0 | 62.2 | 60 | 2.144d | 1395.8445d || 220.83351d
}

def "Calculate the ground reflection eRefS"() {
Expand Down