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.

Sign up for GitHub

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

Jump to bottom

rewrite functions in pv_model #689

Draft
wants to merge 18 commits into
base: dev
Choose a base branch
from
+95 −37
Draft

rewrite functions in pv_model #689

Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
18 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
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
73 changes: 50 additions & 23 deletions src/main/scala/edu/ie3/simona/model/participant/PvModel.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,8 @@ final case class PvModel private (
delta,
lat,
gammaE,
alphaE
alphaE,
duration
)

// === Diffuse Radiation Parameters ===//
Expand Down Expand Up @@ -333,7 +334,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 •
edited
Loading

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.

0.000077 * sin(2d * jInRad)

// solar constant in W/m2
Expand Down Expand Up @@ -405,32 +406,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 @@ -460,13 +456,24 @@ 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
alphaE: Angle,
duration: Time
): Irradiation = {

omegas match {
Expand All @@ -478,6 +485,8 @@ 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°
Copy link
Member

Choose a reason for hiding this comment

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

I'm not fully convinced by your proposal yet.

The original problem is that the ratio between a and b can tend to infinity because the radiation on the horizontal surface goes to 0, while the radiation on tilted surface can stay fairly big if tilted towards sunset.

The solution proposed with R_{b,ave} = a/b is to calculate some sort of average over a time period between omega1 and omega2. This way, the extreme values calculated when using just one point in time close to sunrise or sunset gets reduced. This works because the value further away from sunrise or -set pulls down the average again.

All this was just how it's supposed to work. If now the time slice between omega1 and omega2 is very small, the "relief" we get from the second value is quite reduced. Now it seems like some sort of "crutch on a crutch" or "workaround around a workaround" to try to calculate our way back to the time interval of one hour. This is because calculating an average from a time interval is already a workaround (which could be debatable, I guess). And furthermore this would only work if the change in radiation is proportional to the change in angle omega, which doesn't seem clear to me.

Curious about your thoughts.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I am not quite sure what you mean by "calculate our way back to the time interval of one hour". If the difference between omega1 and omega2 is very small, we are not looking at a full hour, but rather a few minutes. The problem is that the available weather data contains the radiation energy of one full hour in Wh/m^2, which would obviously result in a false result if we multiply it by R_b,ave. This is why there is the linear factor "timeFrame" between 1 and 0, which should solve this problem.

Also, I partially get what you mean by "workaround around a workaround", but doesn't that also apply for the current implemenation aswell? There, the problem is solved by artificially keeping the difference between omega1 and omega2 fixed, which is inaccurate at times close to sunrise or sunset.

Besides, can you please explain why the change in radiation must be proportional to the change in omega?

Furthermore, if you already looked at the integration test, did you notice anything odd? Were there any very high peaks?


val a = ((sin(deltaInRad) * sin(latInRad) * cos(gammaEInRad)
- sin(deltaInRad) * cos(latInRad) * sin(gammaEInRad) * cos(
Expand All @@ -499,7 +508,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 @@ -529,6 +538,27 @@ 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 @@ -551,12 +581,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 •
edited
Loading

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

(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 Expand Up @@ -593,10 +623,7 @@ final case class PvModel private (

// calculate the f_ij components based on the epsilon bin
val f11 = -0.0161 * pow(x, 3) + 0.1840 * pow(x, 2) - 0.3806 * x + 0.2324
val f12 = 0.0134 * pow(x, 4) - 0.1938 * pow(x, 3) + 0.8410 * pow(
x,
2
) - 1.4018 * x + 1.3579
val f12 = 0.0134 * pow(x, 4) - 0.1938 * pow(x, 3) + 0.8410 * pow(x, 2) - 1.4018 * x + 1.3579
val f13 = 0.0032 * pow(x, 3) - 0.0280 * pow(x, 2) - 0.0056 * x - 0.0385
val f21 = -0.0048 * pow(x, 3) + 0.0536 * pow(x, 2) - 0.1049 * x + 0.0034
val f22 = 0.0012 * pow(x, 3) - 0.0067 * pow(x, 2) + 0.0091 * x - 0.0269
Expand Down
40 changes: 24 additions & 16 deletions src/test/groovy/edu/ie3/simona/model/participant/PvModelTest.groovy
View file
Open in desktop
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
//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