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

IAM that supports AR coating like Fresnel #1501

Closed
mikofski opened this issue Jul 22, 2022 · 17 comments · Fixed by #1616
Closed

IAM that supports AR coating like Fresnel #1501

mikofski opened this issue Jul 22, 2022 · 17 comments · Fixed by #1616

Comments

@mikofski
Copy link
Member

mikofski commented Jul 22, 2022

Problem

Currently pvlib supports the DeSoto physical model (similar to normal glass), ASHRAE, Martin & Ruiz, and SAPM polynomial, but it doesn't have a pure Fresnel model that allows additional interfaces like an AR coating.

  • DeSoto physical model is most similar to the Fresnel for normal glass but only has one interface, so is limited to IAM curves below it only, while an AR coating would have a greater ρ
  • Martin & Ruiz could be used to approximate an AR coated glass if the correct a_r were known. The default of a_r=0.16 is slightly above the normal glass Fresnel IAM, but an a_r=0.14 seems to match an AR coating with index of refraction of 1.2 most closely.

pvlib_iam

Proposal

a new method in pvl.iam.fresnel_ar(aoi, n_ar=1.2, n_air=1.0, n_glass=1.56) that implements the Fresnel equation

Alternative

Suggest readers to use Martin & Ruiz with a_r=0.14 instead of default.

additional content

PVsyst has switched to Fresnel equations. We can duplicate their methods ignoring additional reflections and the encapsulant layer:
Fresnel-v-ASHRAE

import numpy as np
import matplotlib.pyplot as plt
plt.ion()


# constants
n_glass = 1.56
n_air = 1.0
theta_inc = np.linspace(0, 88, 100)


def snell(theta_1, n1, n2):
    """Snell's equation"""
    sintheta_2 = n1/n2 * np.sin(np.radians(theta_1))
    return sintheta_2, np.degrees(np.arcsin(sintheta_2))


def refl_s(theta_1, theta_2, n1, n2):
    """Fresnel's equation"""
    n1_costheta_1 = n1*np.cos(np.radians(theta_1))
    n2_costheta_2 = n2*np.cos(np.radians(theta_2))
    return np.abs((n1_costheta_1 - n2_costheta_2)/(n1_costheta_1 + n2_costheta_2))**2


def refl_p(theta_1, theta_2, n1, n2):
    """Fresnel's equation"""
    n1_costheta_2 = n1*np.cos(np.radians(theta_2))
    n2_costheta_1 = n2*np.cos(np.radians(theta_1))
    return np.abs((n1_costheta_2 - n2_costheta_1)/(n1_costheta_2 + n2_costheta_1))**2


def refl_eff(rs, rp):
    """effective reflectivity"""
    return (rs+rp)/2


def trans(refl):
    """transmissivity"""
    return 1-refl


def refl0(n1, n2):
    """reflectivity at normal incidence"""
    return np.abs((n1-n2)/(n1+n2))**2


def fresnel(theta_inc, n1=n_air, n2=n_glass):
    """calculate IAM using Fresnel's Law"""
    _, theta_tr = snell(theta_inc, n1, n2)
    rs = refl_s(theta_inc, theta_tr, n1, n2)
    rp = refl_p(theta_inc, theta_tr, n1, n2)
    reff = refl_eff(rs, rp)
    r0 = refl0(n1, n2)
    return trans(reff)/trans(r0)


def ashrae(theta_inc, b0=0.05):
    """ASHRAE equation"""
    return 1 - b0*(1/np.cos(np.radians(theta_inc)) - 1)


def fresnel_ar(theta_inc, n_ar, n1=n_air, n2=n_glass):
    """calculate IAM using Fresnel's law with AR"""
    # use fresnel() for n2=n_ar
    _, theta_ar = snell(theta_inc, n1, n_ar)
    rs_ar1 = refl_s(theta_inc, theta_ar, n1, n_ar)
    rp_ar1 = refl_p(theta_inc, theta_ar, n1, n_ar)
    r0_ar1 = refl0(n1, n_ar)
    # repeat with fresnel() with n1=n_ar
    _, theta_tr = snell(theta_ar, n_ar, n2)
    rs = refl_s(theta_ar, theta_tr, n_ar, n2)
    rp = refl_p(theta_ar, theta_tr, n_ar, n2)
    # note that combined reflectivity is product of transmissivity!
    # so... rho12 = 1 - (1-rho1)(1-rho2) 
    reff = refl_eff(1-(1-rs_ar1)*(1-rs), 1-(1-rp_ar1)*(1-rp))
    r0 = 1-(1-refl0(n_ar, n2))*(1-r0_ar1)
    return trans(reff)/trans(r0)


# plot Fresnel for normal glass and ASHRAE
plt.plot(theta_inc, fresnel(theta_inc))
plt.plot(theta_inc, ashrae(theta_inc))

# calculate IAM for AR with n=1.1 and plot
iam_ar11 = fresnel_ar(theta_inc, n_ar=1.1)
plt.plot(theta_inc, iam_ar11)

# repeat for AR with n=1.2
iam_ar12 = fresnel_ar(theta_inc, n_ar=1.2)
plt.plot(theta_inc, iam_ar12)

# make plot pretty
plt.legend(['Fresnel, normal glass', 'ASHRAE, $b_0=0.05$', 'Fresnel $n_{AR}=1.1$', 'Fresnel $n_{AR}=1.2$'])
plt.title("IAM correction, Fresnel vs. ASHRAE, using basic eqn's")
plt.ylabel('IAM')
plt.xlabel(r'incidence angle $\theta_{inc} [\degree]$')
plt.grid()
plt.ylim([0.55,1.05])
@cwhanse
Copy link
Member

cwhanse commented Jul 23, 2022

+1. This reference might be relevant.

@adriesse
Copy link
Member

I seem to recall from somewhere that PVsyst actually interpolates from a fixed set of pre-calculated values when simulating.

@mikofski
Copy link
Member Author

PVsyst allows a user specified custom IAM v AOI lookup table in the module PAN file, but that presupposes there exist qualified IAM measurements either from a lab or the manufacturer. Otherwise they use Fresnel as of v6.67. See https://www.pvsyst.com/help/iam_loss.htm

@adriesse
Copy link
Member

Yes, what I meant is that they use the Fresnel equations to populate the table for interpolation. At least this is my recollection.

@kdebrab
Copy link
Contributor

kdebrab commented Dec 16, 2022

In #1616 I submitted a proposal to extend the ìam.physical() function with an extra optional argument n_ar to support AR coating.

Main differences with the code proposed above by @mikofski :

  • instead of adding new functions fresnel and fresnel_ar, the existing function physical (which is equivalent with fresnel) is extended with support for AR coating through the optional argument n_ar.
  • internal transmissions are no longer ignored but taken into account.

Main differences with the previous code of iam.physical() (besides support for AR coating):

  • the Fresnel equations for the reflectance for parallel and perpendicular polarized light has been adapted from the tan formulation (as used by Desoto) to the much more common cos formulation (as used everywhere else). As far as I can tell they are mathematically fully equivalent (in combination with Snell's law), but the cos formulation doesn't suffer from 0/0 division for aoi equal to zero. That's why no special handling is needed anymore for angles near zero.
  • the result for 90° is no longer exactly 0, but (on my laptop) 3.47262923e-16, which I think is fully fine. But for that I had to include a_tol in the existing test_physical unit test.
  • (probably not relevant) though aoi is normally expected to be between 0° and 180°, there were apparently also tests for negative aoi. The new code is now also consistent between positive and negative angles for angles outside the range -360°, 360° (though I didn't include tests for it as I doubt whether we should really support that).

Main differences with PVSyst:

  • I don't know the reason but there is still a (very minor) difference between the values in PVSyst compared to the values returned by this function.

@adriesse
Copy link
Member

adriesse commented Dec 16, 2022

  • the result for 90° is no longer exactly 0, but (on my laptop) 3.47262923e-16, which I think is fully fine. But for that I had to include a_tol in the existing test_physical unit test.

For fun you could try using scipy.special.cosdg()

  • I don't know the reason but there is still a (very minor) difference between the values in PVSyst compared to the values returned by this function.

I don't think they do the internal reflections.

@kdebrab
Copy link
Contributor

kdebrab commented Dec 17, 2022

  • the result for 90° is no longer exactly 0, but (on my laptop) 3.47262923e-16, which I think is fully fine. But for that I had to include a_tol in the existing test_physical unit test.

For fun you could try using scipy.special.cosdg()

I didn't know about scipy.special.cosdg(). Using that function one gets indeed exactly 0 for aoi of 90°.

I still don't think it matters as:

  • 3.47262923e-16 is really very close to zero
  • for aoi > 90° we already get exactly 0
  • timing of sunset, sunrise has its uncertainties and depends on atmospheric conditions
  • the sun moves by 1° every 4 minutes
  • the diameter of the sun is 0.53°.

Nevertheless, maybe pvlib should in general use these scipy.special functions instead of the pvlib.tools functions.

  • I don't know the reason but there is still a (very minor) difference between the values in PVSyst compared to the values returned by this function.

I don't think they do the internal reflections.

Results are actually closer to PVSyst when taking internal reflections into account compared to when not taking them into account. Compare e.g. the results for n_ar = 1.29 (and L=0):

aoi pvsyst with reflections without reflections
0 1 1 1
30 0.999 0.999 0.999
50 0.987 0.987 0.987
60 0.962 0.963 0.962
70 0.892 0.892 0.890
75 0.816 0.814 0.812
80 0.681 0.679 0.675
85 0.440 0.437 0.434
90 0 0 0

So, I'd guess also PVSyst takes internal reflections into account.

@adriesse
Copy link
Member

Are all the relevant parameter values documented in PVsyst? I can't find them on their web site.

@mikofski
Copy link
Member Author

I’m -1 on scipy.special.cosdg see previous thread.

great to see closer agreement with PVsyst, not that it’s the reference, but still it’s de facto. Their online help says they consider reflections, but in my pseudo-code I ignored it. What is L parameter?

@kdebrab
Copy link
Contributor

kdebrab commented Dec 17, 2022

great to see closer agreement with PVsyst, not that it’s the reference, but still it’s de facto. Their online help says they consider reflections, but in my pseudo-code I ignored it. What is L parameter?

L is the thickness of the glass. It's an optional parameter of iam.physical, which is by default 2 mm (though most PV panels are actually 3.2 mm thick). It has a minor impact on the IAM values. You considered it zero as well in your calculation above. I don't know what PVSyst considers but taking it into account reduces the IAM values and thus causes an increase of the deviations with regard to PVSyst.

@kdebrab
Copy link
Contributor

kdebrab commented Dec 17, 2022

Are all the relevant parameter values documented in PVsyst? I can't find them on their web site.

The refractive index for the glass and AR coating can be selected in the tool and its default values (n_glass=1.526, n_ar=1.29) for the Incidence Angle mode 'Fresnel, Glass with AF' can be seen in some threads on the PVsyst forum, e.g. here and here. I couldn't find any values for the other two parameters (glazing thickness L and glazing extinction coefficient K), but that's possibly because PVsyst ignores the absorption in the IAM calculation. That's why I assumed L=0 for the comparison table (and the glazing extinction coefficient doesn't matter any more when the glazing thickness is zero). Moreover, assuming a non-zero glazing thickness would only lower the values, and we need higher values to get closer to the PVsyst values.

@adriesse
Copy link
Member

Well, from my perspective if it is clear you implement these formulas correctly for a two-layer window then it doesn't need to match PVsyst or claim that it is doing the same as PVsyst.

@kandersolar
Copy link
Member

PVWatts v5 and SAM both have the option of a second snell & fresnel iteration for AR coatings, although I don't think either of them consider repeated internal reflections. I point these out as potentially useful inclusions in the references and possible influences on the name of the new function.

@adriesse
Copy link
Member

adriesse commented Jan 18, 2023

L is the thickness of the glass. It's an optional parameter of iam.physical, which is by default 2 mm (though most PV panels are actually 3.2 mm thick).

This is one default that has always bothered me out of principle. For the incorrect value there is an academic paper that can be referenced, thus giving someone else the responsibility, but for the correct value we haven't got one!

@cwhanse
Copy link
Member

cwhanse commented Jan 18, 2023

L is the thickness of the glass. It's an optional parameter of iam.physical, which is by default 2 mm (though most PV panels are actually 3.2 mm thick).

This is one default that has always bothered me out of principle. For the incorrect value there is an academic paper that can be referenced, thus giving someone else the responsibility, but for the correct value we haven't got one!

I'd be OK changing it, assuming explanatory comments are added.

The only reason to stick with L=0.002 m is to stay consistent with the reference. The basis for that value: it was what W. De Soto used in his Master's thesis (same title as the paper), with no further discussion provided. So we shouldn't regard 2mm as the result of any concerted research.

@kdebrab
Copy link
Contributor

kdebrab commented Feb 1, 2023

L is the thickness of the glass. It's an optional parameter of iam.physical, which is by default 2 mm (though most PV panels are actually 3.2 mm thick).

This is one default that has always bothered me out of principle. For the incorrect value there is an academic paper that can be referenced, thus giving someone else the responsibility, but for the correct value we haven't got one!

FYI, the yearly 'International Technology Roadmap for Photovoltaics (ITRPV)' (donwloadable from https://www.vdma.org/international-technology-roadmap-photovoltaic) has a figure about the "world market share of front glass thickness in modules". Below you find the relevant figure in the latest report from November 2022.

image

@adriesse
Copy link
Member

adriesse commented Feb 1, 2023

Good reference to have. I suspect that the increasing fraction of thinner front glass is related to an increasing fraction of glass-glass modules.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants