is vf_row_sky correct?

[mikofski]

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 146 in 7e88d21

y = 0.5 * (cosd(psi_t_shaded) + cst)

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 is psi_t_shaded so finally this angle from the slant height all the way up to the top of the next row is surface_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°).

[mikofski]

ditto:

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 153 in 7e88d21

y = 0.5 * (cosd(psi_t_unshaded) + cst)

and

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 251 in 7e88d21

return 0.5 * (cosd(psi_t_shaded) - cst)

[mikofski]

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.

[cwhanse]

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.

[kandersolar]

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$.

[mikofski]

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!

[cwhanse]

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.

[mikofski]

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:

[mikofski]

Hi @cwhanse did you mean to reopen this?

[cwhanse]

Yes. I think we need to better understand these two different calculations of view factors before we consider this question resolved.

[mikofski]

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.

[mikofski]

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.

[mikofski]

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.

[mikofski]

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

[cwhanse]

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.

[mikofski]

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:

• B-71 is correct for the ground -> sky
• B-1 is correct for module -> sky with no obstructions
• A-10 integrated from zero to 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)

[cwhanse]

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.

[kandersolar]

The orientation of each differential matters.

D'oh! Of course, thanks :)

[mikofski]

icymi: @jdnewmil