pvlib-python icon indicating copy to clipboard operation
pvlib-python copied to clipboard

IAM that supports AR coating like Fresnel

Open mikofski opened this issue 2 years ago • 4 comments

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])

mikofski avatar Jul 22 '22 23:07 mikofski

+1. This reference might be relevant.

cwhanse avatar Jul 23 '22 00:07 cwhanse

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

adriesse avatar Jul 24 '22 14:07 adriesse

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

mikofski avatar Jul 24 '22 16:07 mikofski

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

adriesse avatar Jul 24 '22 19:07 adriesse

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.

kdebrab avatar Dec 16 '22 12:12 kdebrab

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

adriesse avatar Dec 16 '22 21:12 adriesse

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

kdebrab avatar Dec 17 '22 15:12 kdebrab

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

adriesse avatar Dec 17 '22 17:12 adriesse

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?

mikofski avatar Dec 17 '22 17:12 mikofski

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 avatar Dec 17 '22 18:12 kdebrab

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.

kdebrab avatar Dec 17 '22 20:12 kdebrab

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.

adriesse avatar Dec 20 '22 10:12 adriesse

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.

kandersolar avatar Dec 28 '22 20:12 kandersolar

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!

adriesse avatar Jan 18 '23 21:01 adriesse

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.

cwhanse avatar Jan 18 '23 21:01 cwhanse

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

kdebrab avatar Feb 01 '23 08:02 kdebrab

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

adriesse avatar Feb 01 '23 09:02 adriesse