Closed mikofski closed 1 year ago
FYI: these formulas are sort of referenced in doi: 10.1109/PVSC40753.2019.8980572 in eqn's (10) & (11), except that they are multipliers meant to undo the VF the previous transposition step.
I disagree and believe that line of code and its description are correct.
The implementation of the infinite_sheds doesn't use the multiplier approach indicated by the paper. The reference frame for all view factor calculations is the ground. The function _vf_row_sky_integ
first calculates the view factor to the visible sky from each point x
along the slant length of the module, which is given by 1/2 (cos(psi_t) + cos(tilt)), where psi_t is the angle from x to the top of the facing row of modules. Then the view factors are integrated over x to return an "average" view factor value.
I've been chewing on this for a while and am not quite sure where the flaw is in Mark's math. I suspect it has something to do with unwittingly adding or deleting a segment from the sky dome, making it not a hemisphere anymore and causing the view factor to be taken as a fraction of the wrong thing, but I don't have any math to back that up. I'd be interested to see a better explanation if anyone has one.
Anyway, regarding the calculation in the code: start with the usual $(1 + \cos(\beta))/2$ view factor for an isolated row and subtract off the portion of the sky blocked by the facing row. That second VF is the sky VF of an upside-down row tilted at psi_t
, so it can be calculated with $(1 + \cos(180-\psi_t))/2 = (1 - \cos(\psi_t))/2$. Subtracting the latter from the former yields $(\cos(\beta) + \cos(\psi_t))/2$.
I suspect it has something to do with unwittingly adding or deleting a segment from the sky dome, making it not a hemisphere anymore and causing the view factor to be taken as a fraction of the wrong thing
~Yes, I think this was my mistake. In my analogy, isotropic like seems to be coming from below the horizon.~ BTW: I wasn't proposing the "multiplier" approach - but it does make me ~think that~ wonder if [edited] the equations in the reference are incorrect ~because the sky dome isn't entirely visible if the reference frame is rotated to the slant height reference frame of the module surface~. Thanks for the insight!
Although I don't see an issue that needs correcting, I think it is worth documenting how that line of code is derived. In the attached diagram, we want the view factor for the angle of view corresponding to the red arc. Measuring angles from the horizontal (positive to the right), the view factor is 1/2 (cos(tilt) - cos(180 - psi)) = 1/2(cos(tilt) + cos(psi)). The second angle (180-psi) results because pvlib uses psi = shading.masking_angle
which always returns a positive value, i.e., psi measured clockwise as indicated in this diagram.
However, if angles are measured from slanted line (positive upwards from the ground), I get the formula above: view factor = 1/2(cos(0) - cos(180 - (psi + tilt))) = 1/2 (1 + cos(psi + tilt)). I think the bridging the gap between these two equations must lie in the direction (sign) of the various angles, and some clever way to separate the '1' into cosine terms.
angles are measured from slanted line (positive upwards from the ground)
That was my mistake I think, but now I think that the way the angles are measured should not be arbitrary. Measuring the angles the way I did (from the slanted line) requires rotating the reference frame from horizontal, which I believe is equivalent to a module surface tilted further away from the isotropic sky dome by psi_t
.
If I just rely on Marion as the reference then we have:
vf = 0.5 * (cos(beta) - cos(180-psi_t)) = 0.5 * (cos(beta) + cos(psi_t))
according to trig identities.
Also the existing formulation (cos(psi_t) + cos(beta))/2
differs but is numerically closer to the Passias shade factor:
# for backside beta = 160-deg:
(1-pvl.shading.sky_diffuse_passias(np.arange(10)))*(1 + np.cos(tilt))/2
# array([0.03015369, 0.03015139, 0.03014451, 0.03013303, 0.03011696,
# 0.03009632, 0.0300711 , 0.03004131, 0.03000696, 0.02996807])
I think it's important to realize that these two formulas do not evaluate the same:
Hi @cwhanse did you mean to reopen this?
Yes. I think we need to better understand these two different calculations of view factors before we consider this question resolved.
Okay, I think I see your point. I was going to say this is given by B-71 from Siegel's and Howell's "Thermal Radiation" 7th edition which is the fraction of radiation from a:
Differential element of any length to surface generated by a line of infinite length parallel to the plane of the element, and moved parallel to itself. Plane of element does not intersect surface.
But the problem is the orientation of dA1
is pointing up, and the plane it lies in never intersects with A2
so this is appropriate for the view factor of radiation leaving the sky and incident on the ground, which using trig identities and reciprocity derives the equation that's in Marion's reference.
However, a differential surface dA1
on the module is inclined at an angle, and the plane it lies in definitely will intersect A2
so I think a more appropriate view factor would actually be B-1 if you reorient yourself so that dA1
inclined by phi1
and A2
is horizontal.
Some other interesting papers I found are Joseph Applebaum's "View Factors of Flat Collectors, Including Photovoltaics, Visible to Partial Sky" in which they show several view factors using different methods. One is the view factor of an fixed tilt collector on an inclined plane.
By the way, this problem can be solved using A-10 by integrating from zero to Y_psi = sin(psi+phi)/sin(psi)
but the integral is nasty. It also looks like superposition could be used if there were view factor from dA1
to the part of A2
starting where psi
intersects it. I'll try to add a figure or use SymPy maybe. My calculus is very rusty.
Testing both formulas against integrating A-10 in this gist, I find better agreement with the formulation in PR #1666 and the original reference doi: 10.1109/PVSC40753.2019.8980572 in eqn's (10) & (11). I believe the current pvlib formulation would only be valid for the ground's view of the sky as in the Marion reference.
@cwhanse I think the reason the existing formulation doesn't work is because the normal of the differential surface dA1
is pointing upward instead of normal to the panel surface and that affects the angle from the normal to a line connecting the two surfaces:
@kanderso-nrel I think the reason that superposition, Fd1->a = Fd1->2 - Fd1*->b
as you suggested, doesn't work because they are not the same differential dA1
as you can tell by the normal from each surface, which is different for the differential rotated by psi
than for the differential rotated by phi
. The orientation of each differential matters.
I think we should consider merging the PR. I think we can use Siegel & Howell as a reference.
Also I should’ve clarified that I did not derive the numerator expression. I took it from work done by Jeff Newmiller as explained in the wiki bifacial notes. I remember it took me awhile to understand it, as I also originally thought it should be a difference of two cosines, vs a cosine of a difference. Evidentially it didn’t stick because I forgot. That was pre-pandemic
I believe the current pvlib formulation would only be valid for the ground's view of the sky
Yes. Herein lies the conceptual mistake I made in the implementation. B-71 in Siegel and Howell is the correct formula, but the angles must be relative to the plane containing the receiver, not to a reference plane such as the ground surface. I concur that #1666 is the valid correction to that line of code.
Thanks @cwhanse for running this to ground. I will copy some of this into the wiki bifacial notes for reference. Just for clarification I think:
Y_psi
is correct for module -> sky with obstruction above the horizon which approaches B-1 as Y_psi -> inf
where Y_psi = sin(psi+phi)/sin(psi)
B-71 is correct for module -> sky, but in its cosine form since we want to measure angles relative to the module's plane rather than its normal. I wouldn't point to B-1 since that assumes the radiating surface is a flat plane, as opposed to a surface formed by moving a parallel line to form an arbitrary surface. There's an equivalency between B-71 and B-1, but that's one more conceptual step; B-71 and B-1 yield the same formula when the angles are determined relative to the module's plane.
The orientation of each differential matters.
D'oh! Of course, thanks :)
icymi: @jdnewmil
https://github.com/pvlib/pvlib-python/blob/7e88d212c786d0ad334dce6fcafaf29339ff60ab/pvlib/bifacial/infinite_sheds.py#L146
I think this should be:
$$\frac{1 + \cos \left( \text{surface tilt} + \psi_{t}\ \text{shaded} \right)}{2}$$
because in the reference frame of the module surface the angle pointing along the slant height to the sky is actually zero, $cos(0) = 1$, and the angle above the slant height to a horizontal line would be the
surface_tilt
itself, then the angle from the horizontal to the top of the next row ispsi_t_shaded
so finally this angle from the slant height all the way up to the top of the next row issurface_tilt + psi_t_shaded
:For example, this is why if
psi_t_shaded
is zero, then the view factor should collapse to the isotropic view factor $(1+\cos(\beta))/2$ as given on the PVPMC website modeling reference for POA sky diffuse.The actual value difference between the two formulas can be quite small when
psi_t_shaded
is close to zero (eg less than 5°), but it's significant when as the masking angle is larger (eg greater than 5°).