Detailed illumination model for bifacial solar cells

Klaus Jäger; Klaus Jäger; Peter Tillmann; Peter Tillmann; Christiane Becker

doi:10.1364/OE.383570

1. Introduction

Solar energy is the most abundant of all sustainable energy technologies [1]. Further, the cost of photovoltaics (PV) has decreased rapidly in the last decade [2]. Therefore, many studies suggest that the largest fraction of global energy needs in a sustainable carbon-free economy can be covered by photovoltaic solar energy [3,4].

Currently, most photovoltaic modules use monofacial solar cells, which only can utilize light that hits the cell at the front. In contrast to that, bifacial solar cells can also convert light impinging at the back and therefore increase the energy yield significantly. Because of novel silicon solar cell technologies using passivated emitter rear contacts (PERC/PERT/PERL) or IBC contacting schemes, bifacial solar cell operation becomes feasible and indeed, they are expected to have 60% market-share by 2029 [5–7].

To estimate the energy yield of bifacial modules it is vital to develop detailed illumination models, especially because the largest fraction of light that hits the module at the back arises from the ground. Several different illumination models were developed in recent years [8]. Sometimes, modules are treated as isolated and shadowing of module rows in front or behind is not considered [9], which overestimates the irradiation onto the ground. Already several authors developed detailed shadowing models of the ground [10–12], however, often variations of the irradiance along the module are neglected [12]. Calcabrini and coworkers developed a simplified model for front illumination using a sky view factor, which quantifies the landscape around the PV module, and a sun coverage factor, which they define ‘at a location with a raised horizon as the ratio between time that the sun is behind the module or blocked by the skyline per year and the annual sunshine duration at the same location with a clear horizon’ [13].

In a recent publication we minimized the levelized cost of electricity (LCOE) of large PV fields of bifacial modules using a Bayesian optimization algorithm [14]. In the present manuscript, we detail the illumination model, which we used in the preceding manuscript. The model, which is described in section 2 considers direct light and diffuse light from the sky and contains a detailed illumination model of the ground, where the model discriminates between illumination of the ground arising from diffuse and direct light. The model calculates the irradiance components on arbitrarily many positions along the module, which is relevant for finding the minimal irradiance, which determines the PV module performance for many PV modules. Finally, we discuss several examples in section 3.

2. Illumination model

With the illumination model we calculate the irradiance onto a solar module, which is placed somewhere in a big PV-field. We assume this field to be so big that effects from the boundaries can be neglected. Further, we assume the modules to be homogeneous: we neglect effects from the module boundaries or module space in between solar cells. Hence, we can treat this problem as 2-dimensional with periodic boundary conditions, as illustrated in Fig. 1. As input we use the direct normal irradiance (DNI) and the diffuse horizontal irradiance (DHI), which we obtain from freely available weather databases, and the solar position. In [14], DNI and DHI were adapted according to the Perez model to account for the diffuse irradiance more accurately [15].

Fig. 1. Illustrating the irradiance components, which reach the module on the front and back sides; the front components are summarized in Table 1. Here, we assume w.l.o.g. that the PV system is located on the northern hemisphere and oriented towards South, $\phi _m = 180^\circ$. For better visibility, the point $P_m$, on which the irradiance components are evaluated, is not depicted on module $\#0$ but on $\#-1$.

Download Full Size | PDF

Table 1. Irradiance components constituting the illumination of a solar module with $s=\overline {BP_m}$, where $P_m$ is defined in Fig. 1 . These components have to be considered for front and back sides – hence eight components in total.

View Table | View all tables in this article

For the front and back sides of the module, we have to consider four irradiance components each – hence eight components in total. These components are summed up in Table 1 and depicted in Fig. 1. In this model we assume the solar modules to be completely black, meaning they do not reflect any light that could reach another module.

The total irradiance (or intensity; given in watt per m$^2$) on front is given by

(1)$$I_f(s) = I_{\textrm{dir, }f}^{\textrm{sky}}(s) + I_{\textrm{diff, }f}^{\textrm{sky}}(s) + I_{\textrm{dir, }f}^{\textrm{gr.}}(s) + I_{\textrm{diff, }f}^{\textrm{gr.}}(s),$$

and similar for the back side with a subscript $b$ instead of $f$.

In section 2.1 we will derive expressions for the four components impinging onto the solar module from the sky. Expressions for the components impinging from the ground are derived in section 2.2. Table 2 shows the input parameters required for calculating the components.

Table 2. Input parameters for the illumination model.

View Table | View all tables in this article

The irradiance $I$ is based on the more fundamental notion radiance $L$ (in watt per steradian per m$^2$), which is connected to the irradiance via

(2)$$I(\mathbf{r}) = \int_\Omega L(\mathbf{r},\theta,\phi) \cos\theta\,\,\textrm{d}\Omega,$$

where $\mathbf {r}$ defines the coordinates on the surface on which $I$ and $L$ are evaluated [16]. $\theta$ and $\phi$ denote the polar angle and azimuth, where $\theta$ is defined such that it denotes the angle between the evaluated ray and the normal to the surface at $\mathbf {r}$. The term $\cos \theta$ arises from Lambert’s cosine law, which accounts for the decreasing irradiance of a beam of light which strikes a surface under increasing angle, because the radiant power is distributed across a larger area. The solid angle element is given by $\,\textrm {d}\Omega = \sin \theta \,\,\textrm {d}\theta \,\,\textrm {d}\phi$.

For the following derivation we define dimensionless geometrical distribution functions as

(3)$$\begin{array}{ccc} \iota_{\textrm{dir}}(s) := \frac{I_{\textrm{dir}}(s)}{\mathrm{DNI}} &\textrm{and} & \iota_{\textrm{diff}}(s) := \frac{I_{\textrm{diff}}(s)}{\mathrm{DHI}} \end{array}$$

for the components arising from direct sunlight and diffuse skylight, respectively. Here we omitted the superscripts “sky” and “gr.” as well as “$f$” and “$b$”. DNI is the direct normal irradiance and DHI denotes the diffuse horizontal irradiance; see also Table 2.

2.1 Irradiance components from the sky

The direct irradiance components from the sky are given by

(4a)$$\begin{aligned} I_{\textrm{dir, }f}^{\textrm{sky}}(s) &= \begin{cases} \textrm{DNI}\cdot\phantom{|}\cos\sigma_{mS}\phantom{|}, & \textrm{if } \cos\sigma_{mS}\;>\;0,\\ 0, & \textrm{if } \cos\sigma_{mS}\leq 0, \end{cases} \end{aligned}$$

(4b)$$\begin{aligned} I_{\textrm{dir, }b}^{\textrm{sky}}(s) &= \begin{cases} \textrm{DNI}\cdot|\cos\sigma_{mS}|, & \textrm{if } \cos\sigma_{mS}\;<\;0,\\ 0, & \textrm{if } \cos\sigma_{mS}\geq 0. \end{cases} \end{aligned}$$

where the cosine between sunlight and module normal is calculated with

(5)$$\cos\sigma_{mS} = \mathbf{n}_m\cdot\mathbf{n}_S = \begin{pmatrix} \cos\phi_m\sin\theta_m \\ \sin\phi_m\sin\theta_m \\ \cos\theta_m\end{pmatrix}\cdot \begin{pmatrix} \cos(-\phi_S) \sin\theta_S \\ \sin(-\phi_S) \sin\theta_S \\ \cos\theta_S \end{pmatrix}.$$

The $-\phi _S$ arises from the fact that the solar azimuth $\phi _S$ has a negative (clockwise) orientation in the right-handed coordinate system defined in Fig. 1. Hence, it increases from $0$ at North ($x$-direction) to $\frac {\pi }{2}$ at East ($-y$-direction) and $\pi$ at South. Note that $\cos \theta _{mS}\;>\;0$ when direct sunlight hits the front of the module and $<\;0$ when it hits the back. For modules facing South ($\theta _m=180^\circ$), the latter can occur in the mornings and the evenings in the half year between the vernal equinox and the autumnal equinox. Further, our algorithm checks whether the position $s$ on the module is shadowed by the row of modules in front. In that case, $I_{\textrm {dir, }f}^{\textrm {sky}}(s) = 0$.

For the diffuse irradiance from the sky $I_{\textrm {diff, }f}^{\textrm {sky}}(s)$ on the module front we have to consider the angular range from which light can reach the module. This concept is also known as view factors (VF) in literature [8]. Generally speaking, diffuse light can reach the module front from an opening that spans from the top of the front row $D_{-1}$ to the plane of the module itself, marked by the topmost point $D_0$. In the example in Fig. 1, point $P_m$ on module $\#0$ is illuminated by light from directions within $\sphericalangle D_{-1}P_mD_0$. With respect to the module normal, this range is constrained by the angles $\alpha _{s1} = -\pi /2$ and $\alpha _{s2}(s)$. However, diffuse light does not only reach the module from directions within the $xz$-plane but from a spherical wedge as explained in the appendix. Under the assumption that diffuse light arrives isotropically from the upper hemisphere and using Eqs. (22) and (24) we find

(6a)$$I_{\textrm{diff, }f}^{\textrm{sky}}(s) = \textrm{DHI}\cdot \textrm{SVF}_f^{\textrm{sky}} = \frac{\textrm{DHI}}{2}\left[\sin\alpha_{s2}(s)-\sin\alpha_{s1}\right] = \frac{\textrm{DHI}}{2}\left[\sin\alpha_{s2}(s)+1\right].$$

For the diffuse irradiance from the sky $I_{\textrm {diff, }b}^{\textrm {sky}}(s)$, the SVF for $P_m$ is restricted by $\sphericalangle D_0P_mD_1$. With respect to the module normal, this range is constrained by the angles $\epsilon _{s1}(s)$ and $\epsilon _{s2}=\pi /2$.

(6b)$$I_{\textrm{diff, }b}^{\textrm{sky}}(s) = \textrm{DHI}\cdot\textrm{SVF}_b^{\textrm{sky}} = \frac{\textrm{DHI}}{2}\left[\sin\epsilon_{s2}-\sin\epsilon_{s1}(s)\right] = \frac{\textrm{DHI}}{2}\left[1-\sin\epsilon_{s1}(s)\right]$$

2.2 Irradiance components from the ground

As illustrated in Fig. 2, the ground receives direct and diffuse light from the sky, where we introduce the irradiances $G_{\textrm {dir}}(x_g)$ and $G_{\textrm {diff}}(x_g)$. We assume the ground to be a Lambertian reflector with albedo $A$. Under this assumption the ground emits isotropically with the radiance $L^{\textrm {gr.}}(x_g)$, which is connected to the irradiance $G(x_g)$ via

(7)$$L^{\textrm{gr.}}(x_g) = \frac{A}{\pi}G(x_g),$$

where we omitted the subscripts “diff” and “dir”. Similar to Eq. (3) we define geometrical distribution functions on the ground:

(8)$$\begin{array}{ccc} \gamma_{\textrm{dir}} := \frac{G_{\textrm{dir}}}{\mathrm{DNI}}, &\textrm{and}& \gamma_{\textrm{diff}} := \frac{G_{\textrm{diff}}}{\mathrm{DHI}}=\textrm{SVF}_g. \end{array}$$

Fig. 2. Illustrating the irradiance components, which reach the ground under a large field of PV modules. Here, we assume w.l.o.g. that the PV system is located on the northern hemisphere and oriented towards South.

Download Full Size | PDF

The direct geometrical distribution function is position dependent and takes values of 0 for shaded areas (no direct sunlight) and the projection of the direct irradiance on the ground otherwise:

(9)$$\gamma_{\textrm{dir}}(x_g) = \begin{cases} \small \cos\theta_S, & \textrm{if direct sunlight hits ground at }x_g,\\ 0, & \textrm{if a module blocks direct light at }x_g, \end{cases}$$

hence $\gamma _{\textrm {dir}}$ is independent of the azimuth of the Sun $\phi _S$. In Fig. 2 the fractions of the ground, which can be reached by direct sunlight, are marked in orange. Depending on the geometrical module parameters and the position of the Sun, the illuminated area may lay completely within the unit cell as in the example in Fig. 2, it may extend from one unit cell into the next as in Fig. 3 or no direct light can reach the ground. The latter can occur when the module spacing $d$ decreases or when the zenith angle $\theta _S$ is large, which means that the Sun has a low altitude.

Fig. 3. An example for (a) a module configuration and (b) the corresponding diffuse and direct geometrical distribution functions $\gamma$ on the ground, as defined in Eq. (8). The following parameters were used: $\ell = 1.96$ m, $d=7.30$ m, $h=0.50$ m and $\theta _m=52$°. The solar position for the direct component was $a_S=31.9$° and $\phi _S=144.1$° (Berlin, 20 June 2019, 11:52 am CEST). The number of positions $x_g$ on the ground, on which the functions are evaluated, is controlled by the parameter $N_g$, which was 101 in this example. The unit cell is represented as shaded area.

Download Full Size | PDF

The diffuse geometrical distribution function describes the fraction of the sky that is visible at a certain location on the ground. For the example shown in Fig. 2, light can reach $P_g$ from three ranges: $\sphericalangle D_{-1}P_gB_0$, $\sphericalangle D_0P_gD_1$ and $\sphericalangle B_1P_gD_2$. Using Eq. (22), we find

(10)$$\begin{aligned} \textstyle \gamma_{\textrm{diff}}(x_g) = \frac{\textrm{1}}{2}[ & \left(\sin\delta_{-1}-\sin\beta_0\right)\\ + & \left(\sin\delta_0-\sin\delta_1\right)\\ + & \left(\sin\beta_1-\sin\delta_2\right)]. \end{aligned}$$

If $|\delta _{-1}|\leqslant |\beta _0|$ no light reaches the ground from between modules $\#-1$ and $\#0$. If $x_g\;>\;x'$, we have to use the ranges $\sphericalangle D_0P_gB_1$ and $\sphericalangle D_1P_gD_2$ instead of the ones mentioned above.

Now we derive expressions for the irradiance components, which illuminate the module from the ground: $I_{\textrm {dir, }f}^{\textrm {gr.}}$, $I_{\textrm {diff, }f}^{\textrm {gr.}}$, $I_{\textrm {dir, }b}^{\textrm {gr.}}$ and $I_{\textrm {diff, }b}^{\textrm {gr.}}$. The following derivation is valid for both components originating either from direct sunlight or diffuse skylight. Hence, we omit the subscripts “dir” and “diff” in the following. As illustrated in Fig. 1, the relevant angular range is given by $\sphericalangle B_{-1}P_mB_0$ for the front and $\sphericalangle B_0P_mB_1$ for the back, respectively. With respect to the module normal, this range is constrained by the angles $\alpha _{g1}(s)$ and $\alpha _{g2} = \pi /2$ for the front and by the angles $\epsilon _{g1}=-\pi /2$ and $\epsilon _{g2}(s)$ for the back. The diffuse irradiance components from the ground can be calculated using Eq. (21) in the appendix,

(11a)$$\begin{aligned} I_f^{\textrm{gr.}}(s) &= \frac{\pi}{2}\int_{\alpha_{g1}(s)}^{\alpha_{g2}}L^{\textrm{gr.}}\left[x^g_f\left(s,\alpha\right)\right]\cos\alpha\,\,\textrm{d}\alpha, \end{aligned}$$

(11b)$$\begin{aligned} I_b^{\textrm{gr.}}(s) &= \frac{\pi}{2}\int_{\epsilon_{g1}}^{\epsilon_{g2}(s)}L^{\textrm{gr.}}\left[x^g_b\left(s,\epsilon\right)\right]\cos\epsilon\,\,\textrm{d}\epsilon. \end{aligned}$$

As illustrated in Fig. 2, $x_g(s,\alpha )$ and $x_g(s,\epsilon )$ are the coordinates of the respective points $P_g(s,\alpha )$ and $P_g(s,\epsilon )$ on the ground, which are defined such that the angles between the line $\overline {P_gP_m}$ and the module normal $\mathbf {n}_m$ are equal to $\alpha$ or $\epsilon$, respectively. Using $\phi _m=180^\circ$, we find

(12a)$$\begin{aligned} x_g(s,\alpha) &= s\cdot\cos\theta_m + \left(h+s\cdot\sin\theta_m\right)\tan\left(\theta_m+\alpha\right), \end{aligned}$$

(12b)$$\begin{aligned} x_g(s,\epsilon) &= s\cdot\cos\theta_m + \left(h+s\cdot\sin\theta_m\right)\tan\left(\theta_m-\epsilon\right). \end{aligned}$$

Applying the definitions from Eqs. (3) and (8) together with Eq. (7) to Eqs. (11) yields

(13a)$$\begin{aligned}\iota_f^{\textrm{gr.}}(s) &= \frac{A}{2}\int_{\alpha_{g1}(s)}^{\alpha_{g2}}\gamma\left[x_g\left(s,\alpha\right)\right]\cos\alpha\,\,\textrm{d}\alpha, \end{aligned}$$

(13b)$$\begin{aligned}\iota_b^{\textrm{gr.}}(s) &= \frac{A}{2}\int_{\epsilon_{g1}}^{\epsilon_{g2}(s)}\gamma\left[x_g\left(s,\epsilon\right)\right]\cos\epsilon\,\,\textrm{d}\epsilon. \end{aligned}$$

2.3 Implementation details

The implementation is written in Python and utilizes the NumPy scientific computing package fast tensor operations. The code was implemented with the fast evaluation of time-series data for energy yield calculations in mind. For a fixed-array geometry and changing solar positions only the geometric distribution functions for direct light $\iota _{\textrm {dir}}$ have to be calculated separately for every timestamp while the different $\iota _{\textrm {diff}}$ are time independent:

(14)$$\begin{array}{ccc}I_{\textrm{dir}}(s,t) := \mathrm{DNI}(t) \cdot \iota_{\textrm{dir}}(s,t) &\textrm{and}& I_{\textrm{diff}}(s,t) := \mathrm{DHI}(t) \cdot \iota_{\textrm{diff}}(s) \end{array}$$

The irradiance of direct and diffuse light from the sky, which directly hit the solar modules, can be calculated with simple geometric functions. These numeric functions are evaluated very fast with vectorized code.

To account for reflected light from the ground the integral in Eq. (13) would need to be solved for every timestep. To avoid this, the ground is discretized into a finite amount of elements $N_g$. Within each element the radiance is assumed to be constant. This allows to reformulate Eq. (13) into a sum over all elements visible to the module. By introducing the ground view factor VF$_g$ the geometrical calculations are not depending on the timestamp anymore, which is very beneficial for the speed of numerical evaluation:

(15)$$I(s,t)=\frac{A}{2}\int_{\alpha_{g1}(s)}^{\alpha_{g2}}\gamma\left[x_g\left(s,\alpha\right),t\right]\cos\alpha\,\,\textrm{d}\alpha = \sum_{i=1}^n\gamma\left[x_g\left(s,\alpha_i\right),t\right]\cdot\textrm{VF}_g^i$$

where the VF$_g$ is given by

(16)$$\textrm{VF}_g^i=\int_{\alpha_i}^{\alpha_{i+1}}\cos\alpha\,\,\textrm{d}\alpha=\sin\alpha_{i+1}-\sin\alpha_i$$

The numerical implementation allows to perform a one year simulation with irradiance data in one-hour intervals (4400 data points) in less then 1 s on normal personal computer. The code was published on github and can be accessed online [17].

3. Some examples

In this section we demonstrate the model with a few examples. We assume a PV system which is located in Berlin, Germany (52.5°N, 13.25°E). Based on a rule-of-thumb rule we set the module tilt approximately to the geographical latitude, $\theta _m = 52^\circ$ and determine the distance such that a module does not shade a module behind on 21 December, 12 noon [11,18]. For modules with $\ell = 1.96$ m length this leads to a distance between module rows of $d=7.3$ m. A sketch of this configuration is shown in Fig. 3(a). The solar positions are calculated using the Python package pysolar [19].

Figure 3(b) shows the geometrical distribution functions $\gamma$ for direct and diffuse components on the ground for 20 June 2019, 11:52 am CEST. $\gamma _{\textrm {diff}}$ is minimal below the module where the angle covered by the module is largest; and maximal at $x'$, because here the ground sees (almost) no shadow from module #1.

Figure 4 shows how the direct geometrical distribution function on the ground $\gamma _{\textrm {dir}}$ develops during three days (20 June, 23 September and 20 November 2019). For 20 June, a large fraction of the ground is illuminated throughout the day. After sunrise and before sunset the Sun is below the module plane and direct sunlight hits the module back. During short amounts of time in the morning and evening the Sun lies in the module plane and no direct sunlight hits the module. At these times the shadow on the ground vanishes (around 7 am and 7 pm). In between these times the Sun is above the plane and direct light hits the module front. On 23 September, the length of the shadow is almost the same during the whole day: the module tilt ($52^\circ$) is almost the same as the latitude of Berlin ($52.5^\circ )$ - hence at equinox the ecliptic is practically normal to the module plane. On 20 November, only short stretches of the ground are illuminated.

Fig. 4. Illumination of the ground by direct sunlight on three days: (a) 20 June, (b) 23 September, and (c) 20 November 2019. The color corresponds to the geometrical distribution function $\gamma _{\textrm {dir}}$ and is equivalent to $\cos (\theta _S)$, if it is not zero due to shadowing [see Eq. (9)]. The following parameters were used: $\ell = 1.96$ m, $d=7.30$ m, $h=0.50$ m and $\theta _m=52$°. The lines in (a) correspond to the timestamps also used in Figs. 3(b), 5 and 6. The number of positions $x_g$ on the ground, on which the functions are evaluated, is controlled by the parameter $N_g$, which was 101 in this example.

Download Full Size | PDF

Figure 5 shows an example for the eight geometrical distribution functions $\iota$ corresponding to the irradiance components hitting the PV module on its front and back sides. While the functions originating from the sky (a) are stronger on the front side, the components originating from the ground (b) are stronger on the back side. This can be understood by the opening angles: the opening angle towards the sky is larger on the front side, but the opening angle of the ground is larger at the back. The largest relative variations are observed for the ground component at the module back side originating from direct sunlight: it is strongest for short $s$-values. This can be understood when looking at Fig. 3(b): the ground is illuminated by sunlight directly below the module.

Fig. 5. Example of geometrical distribution functions $\iota$ on the module for light the module receives (a) from the sky and (b) the ground. The following parameters were used: $\ell = 1.96$ m, $d=7.30$ m, $h=0.50$ m, $\theta _m=52$°, and albedo $A=30\%$. The solar position for the direct components was $a_S=31.9$° and $\phi _S=144.1$° (Berlin, 20 June 2019, 11:52 am CEST). For the calculation we used DHI $= 1$ and DNI $= 1$. The number of positions $s$ along the module, on which the functions are evaluated, is controlled by the parameter $N_m$, which is 12 in this example.

Download Full Size | PDF

In Fig. 6 we look at two examples on how the irradiance distribution along the module varies depending on the irradiation conditions. We utilized irradiation data measured at the Hochschule für Technik und Wirtschaft Berlin on 20 June 2019. DNI data was measured with the SHP1 pyrheliometer and DHI data was obtained with the SMP1 pyranometer [20]. Figure 6(a) shows results for 11:52 am CEST. Then, no DNI (0 W/m$^2$) was measured while DHI was 144 W/m$^2$ – hence clouds covered the Sun. As discussed in [14], for many PV modules, such as silicon modules composed of several single solar cells connected in series, the overall module performance is not determined by the mean irradiance but by the lowest irradiance observed by the module. Both the front and total irradiation are lowest at the lower end of the module. The total irradiance at the upper end is 10.7% higher than at the lower end. Figure 6(b) shows the situation 46 minute later, at 12:38 pm. Then, the DNI (883 W/m$^2$) is much higher than the DHI (134 W/m$^2$) – direct sunlight hits the PV module. Here, the irradiation onto the module front is almost independent of the module position, the maximum variations in total irradiance are below 0.5%. The total irradiation, which would be relevant for bifacial modules, is minimal at $s\approx 0.57$ m

Fig. 6. Examples of irradiation onto the module for irradiation conditions in Berlin on 20 June 2019 at (a) 11:52 am and (b) 12:38 pm CEST. The following parameters were used: $\ell = 1.96$ m, $d=7.350$ m, $h=0.50$ m, $\theta _m=52$°, and albedo $A=30\%$. The radiation data was obtained from HTW Berlin weather station [20]. The number of positions $s$ along the module, on which the functions are evaluated, is controlled by the parameter $N_m$, which is 12 in this example. In (a), also an example for the total irradiance mean is shown, which would be calculated by models based on a view factor approach.

Download Full Size | PDF

The varying position of the irradiance minimum shows that it is very relevant to determine the irradiance components as functions along the module position $s$, which is automatically done by the model presented in this manuscript. Models based on view factors often have only mean irradiances as output, which might overestimate the final module performance.

4. Conclusions and outlook

Accurate illumination models are crucial for yield estimations of PV power plants with bifacial solar cells. In this manuscript, we derived a detailed model that takes all relevant irradiation components onto PV modules into account and allows for a quick calculation. The code can be accessed online [17]. The irradiation components are calculated locally resolved along the module, which enables the user to determine the minimal irradiance, which is dependent on the illumination conditions, as illustrated in Fig. 6. This is relevant for accurately estimating the performance for many types of PV modules.

In a next step, we will expand the model such that the irradiance components onto the module are angular resolved. This is especially relevant when connecting the illumination model to optical models of solar cells. Further, it would be desirable to experimentally assess the model described in this manuscript in the future.

5. Appendix. Integration of radiance across a spherical wedge

As shown in Fig. 1, the different diffuse components reach a point $P_m$ on the solar module from certain angular ranges. However, diffuse light does not only reach the module from directions within the $xz$-plane, but from a spherical wedge, as illustrated in Fig. 7.

Fig. 7. The spherical wedge indicating the directions from which diffuse light can reach the module. Here, the wedge angle is $(\alpha _2-\alpha _1)$, where $\alpha _1$ and $\alpha _2$ denote angles with respect to the normal direction $\mathbf {n}$. Depending on the application, $\mathbf {n}$ is either the module normal $\mathbf {n}_m$ or the normal to the ground $\mathbf {n}_g$. When $\textbf {n}=\textbf {n}_g$, the directions $\xi$ and $\zeta$ coincide with $x$ and $z$ in Fig. 1. When $\textbf {n}=\textbf {n}_m$, $\xi$ and $\zeta$ are rotated with respect to $x$ and $z$ by the module tilt angle $\theta _m$. The $xz$- and $\xi \zeta$-planes always coincide.

Download Full Size | PDF

The irradiance can be calculated with

(17)$$I(P) = \int_\Omega L(\alpha,\beta) \cos\theta\,\,\textrm{d}\Omega,$$

where $\Omega$ is the set of solid angles defined by the wedge and $\theta$ is the angle spanned between the directional unit vector $\mathbf {e}(\alpha ,\beta )$ and the normal $\mathbf {n}$. With the spherical coordinate system defined in Fig. 7, we have

(18)$$\,\textrm{d}\Omega = \sin\beta\,\,\textrm{d}\alpha\,\,\textrm{d}\beta.$$

The cosine factor can be expressed as

(19)$$\cos\theta = \mathbf{n}\cdot\mathbf{e}(\alpha,\beta) = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \cdot \begin{pmatrix} \sin\beta\sin\alpha \\ \cos\beta \\ \sin\beta\cos\alpha \end{pmatrix} = \sin\beta\cos\alpha.$$

Hence, we find for the irradiance

(20)$$I(P) = \int_{\alpha = \alpha_1}^{\alpha_2}\int_{\beta = 0}^\pi L(\alpha,\beta)\cos\alpha\sin^2\beta\,\,\textrm{d}\alpha\,\,\textrm{d}\beta.$$

As we assume that the PV field is so big that boundary effects can be neglected. As a consequence, the illumination is invariant under translation along the $y$-axis and the problem can be treated as 2-dimensional. Therefore, the radiance does not depend on $\beta$, $L(\alpha ,\beta ) \equiv L(\alpha )$. Integrating over $\beta$ leads to

(21)$$I(P) = \frac{\pi}{2}\int_{\alpha_1}^{\alpha_2}L(\alpha)\cos\alpha\,\,\textrm{d}\alpha.$$

When the radiance is isotropic, $L(\alpha ) \equiv L$, we have

(22)$$I(P) = \frac{L\pi}{2}\left(\sin\alpha_2-\sin\alpha_1\right).$$

For small angle intervals it is sometimes useful to work with the angles $\alpha = 0.5(\alpha _1 + \alpha _2)$ and $\delta = \alpha _2-\alpha _1$. Then, Eq. (22) becomes

(23)$$I(P) = L\pi\cos\alpha\sin\frac{\delta}{2}.$$

In our model, the diffuse light from the sky is given as diffuse horizontal irradiance DHI. In this work we assume the diffuse skylight to be isotropic, hence the diffuse sky radiance $L_s$ is independent of the direction. With Eq. (22), we find

(24)$$L_s = \frac{1}{\pi}\textrm{DHI},$$

where we used $\alpha _1 = -\frac {\pi }{2}$ and $\alpha _2 = \frac {\pi }{2}$.

Funding

Helmholtz Einstein International Berlin Research School in Data Science (HEIBRiDS); Helmholtz Association via Helmholtz Excellence Network SOLARMATH (ExNet-0042-Phase-2-3).

Acknowledgments

We thank Lev Kreinin and Asher Karsenti from SolAround for fruitful discussions regarding the illumination model for bifacial solar cells. The results were obtained at the Berlin Joint Lab for Optical Simulations for Energy Research (BerOSE) of Helmholtz-Zentrum Berlin für Materialien und Energie, Zuse Institute Berlin and Freie Universität Berlin, and within the Helmholtz Excellence Network SOLARMATH, a strategic collaboration of the DFG Excellence Cluster MATH+ and Helmholtz-Zentrum Berlin für Materialien und Energie.

Disclosures

The authors declare no conflict of interest.

References

1. F. Creutzig, P. Agoston, J. C. Goldschmidt, G. Luderer, G. Nemet, and R. C. Pietzcker, “The underestimated potential of solar energy to mitigate climate change,” Nat. Energy 2(9), 17140 (2017). [CrossRef]

2. N. M. Haegel, H. Atwater, T. Barnes, C. Breyer, A. Burrell, Y.-M. Chiang, S. D. Wolf, B. Dimmler, D. Feldman, S. Glunz, J. C. Goldschmidt, D. Hochschild, R. Inzunza, I. Kaizuka, B. Kroposki, S. Kurtz, S. Leu, R. Margolis, K. Matsubara, A. Metz, W. K. Metzger, M. Morjaria, S. Niki, S. Nowak, I. M. Peters, S. Philipps, T. Reindl, A. Richter, D. Rose, K. Sakurai, R. Schlatmann, M. Shikano, W. Sinke, R. Sinton, B. Stanbery, M. Topic, W. Tumas, Y. Ueda, J. van de Lagemaat, P. Verlinden, M. Vetter, E. Warren, M. Werner, M. Yamaguchi, and A. W. Bett, “Terawatt-scale photovoltaics: Transform global energy,” Science 364(6443), 836–838 (2019). [CrossRef]

3. M. Z. Jacobson, M. A. Delucchi, Z. A. Bauer, S. C. Goodman, W. E. Chapman, M. A. Cameron, C. Bozonnat, L. Chobadi, H. A. Clonts, P. Enevoldsen, J. R. Erwin, S. N. Fobi, O. K. Goldstrom, E. M. Hennessy, J. Liu, J. Lo, C. B. Meyer, S. B. Morris, K. R. Moy, P. L. O’Neill, I. Petkov, S. Redfern, R. Schucker, M. A. Sontag, J. Wang, E. Weiner, and A. S. Yachanin, “100% Clean and Renewable Wind, Water, and Sunlight All-Sector Energy Roadmaps for 139 Countries of the World,” Joule 1(1), 108–121 (2017). [CrossRef]

4. M. Child, C. Kemfert, D. Bogdanov, and C. Breyer, “Flexible electricity generation, grid exchange and storage for the transition to a 100% renewable energy system in Europe,” Renew. Energy 139, 80–101 (2019). [CrossRef]

5. T. Dullweber, C. Kranz, R. Peibst, U. Baumann, H. Hannebauer, A. Fülle, S. Steckemetz, T. Weber, M. Kutzer, M. Müller, G. Fischer, P. Palinginis, and H. Neuhaus, “PERC+: industrial PERC solar cells with rear Al grid enabling bifaciality and reduced Al paste consumption,” Prog. Photovoltaics 24(12), 1487–1498 (2016). [CrossRef]

6. S. Chunduri and M. Schmela, “Bifacial Solar Technology Report 2018 edition,” Tech. rep., TaiYang News (2018).

7. “International Technology Roadmap for Photovoltaic (ITRPV), Tenth edition,” Tech. rep., VDMA (2019).

8. T. S. Liang, M. Pravettoni, C. Deline, J. S. Stein, R. Kopecek, J. P. Singh, W. Luo, Y. Wang, A. G. Aberle, and Y. S. Khoo, “A review of crystalline silicon bifacial photovoltaic performance characterisation and simulation,” Energy Environ. Sci. 12(1), 116–148 (2019). [CrossRef]

9. R. Schmager, M. Langenhorst, J. Lehr, U. Lemmer, B. S. Richards, and U. W. Paetzold, “Methodology of energy yield modelling of perovskite-based multi-junction photovoltaics,” Opt. Express 27(8), A507–A523 (2019). [CrossRef]

10. U. A. Yusufoglu, T. M. Pletzer, L. J. Koduvelikulathu, C. Comparotto, R. Kopecek, and H. Kurz, “Analysis of the Annual Performance of Bifacial Modules and Optimization Methods,” IEEE J. Photovolt. 5(1), 320–328 (2015). [CrossRef]

11. L. Kreinin, A. Karsenty, D. Grobgeld, and N. Eisenberg, “PV systems based on bifacial modules: Performance simulation vs. design factors,” in Proceedings of the 43rd IEEE Photovoltaic Specialists Conference (PVSC), (IEEE, 2016), pp. 2688–2691.

12. B. Marion, S. MacAlpine, C. Deline, A. Asgharzadeh, F. Toor, D. Riley, J. Stein, and C. Hansen, “A practical irradiance model for bifacial PV modules,” in Proceedings of the 44th IEEE Photovoltaic Specialist Conference (PVSC), (IEEE, 2017), pp. 1537–1542.

13. A. Calcabrini, H. Ziar, O. Isabella, and M. Zeman, “A simplified skyline-based method for estimating the annual solar energy potential in urban environments,” Nat. Energy 4(3), 206–215 (2019). [CrossRef]

14. P. Tillmann, K. Jäger, and C. Becker, “Minimising levelised cost of electricity of bifacial solar panel arrays using bayesian optimisation,” Sustain. Energy Fuels 4(1), 254–264 (2020). [CrossRef]

15. R. Perez, R. Stewart, R. Seals, and T. Guertin, “The development and verification of the Perez diffuse radiation model,” Tech. rep., Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA (United States) (1988).

16. J. E. Harvey, C. L. Vernold, A. Krywonos, and P. L. Thompson, “Diffracted radiance: A fundamental quantity in nonparaxial scalar diffraction theory,” Appl. Opt. 38(31), 6469–6481 (1999). [CrossRef]

17. P. Tillmann and K. Jäger, Bifacial illumination model, (2019). DOI: 10.5281/zenodo.3543570.

18. P. Grana, “The new rules for latitude and solar system design,” (2018). https://www.solarpowerworldonline.com/2018/08/new-rules-for-latitude-and-solar-system-design/.

19. B. Stafford, “Pysolar,” (2018). DOI: 10.5281/zenodo.1461066.

20. V. Quaschning, T. Tjaden, J. Weniger, and J. Bergner, “HTW weather station,” (2019). https://wetter.htw-berlin.de/History/DIF_SMP5,DIR_SHP1/2019-06-20/2019-06-21.

direct irradiance from the sky	$I_{dir}^{sky} (s)$
diffuse irradiance from the sky	$I_{diff}^{sky} (s)$
diffuse irradiance from the ground…
…originating from direct sunlight	$I_{dir}^{gr.} (s)$
…originating from diffuse skylight	$I_{diff}^{gr.} (s)$

Module parameters (depicted in Fig. 1)
$ℓ$	module length (m)
$d$	module spacing (m)
$h$	module height above the ground (m)
$θ_{m}$	module tilt angle
$ϕ_{m}$	module azimuth angle
Solar and environmental parameters
DNI	direct normal irradiance (W/m $^{2}$ )^¶
DHI	diffuse horizontal irradiance (W/m $^{2}$ )^¶
$θ_{S}$	zenith angle of the Sun
$a_{S}$	solar altitude, connected to zenith via $a_{S} = 90^{\circ} - θ_{S}$
$ϕ_{S}$	azimuth of the Sun
$A$	albedo of the ground
^¶ This parameter also can be spectral. Then, the unit would be W/(m $^{2}$ nm) and we would use the spectral irradiance $[I_{λ}] =$ W/(m $^{2}$ nm) and the spectral radiance $[L_{λ}] =$ W/(sr $\cdot$ m $^{2}$ nm)
Numerical parameters
$N_{g}$	number of points at which irradiance on the ground is evaluated in the interval $[0, d]$
$N_{m}$	number of points at which irradiance is evaluated on the module
$N_{α}$	number of directions at which light impinging onto the module is discretized

direct irradiance from the sky	$I_{dir}^{sky} (s)$
diffuse irradiance from the sky	$I_{diff}^{sky} (s)$
diffuse irradiance from the ground…
…originating from direct sunlight	$I_{dir}^{gr.} (s)$
…originating from diffuse skylight	$I_{diff}^{gr.} (s)$

Module parameters (depicted in Fig. 1)
$ℓ$	module length (m)
$d$	module spacing (m)
$h$	module height above the ground (m)
$θ_{m}$	module tilt angle
$ϕ_{m}$	module azimuth angle
Solar and environmental parameters
DNI	direct normal irradiance (W/m $^{2}$ )^¶
DHI	diffuse horizontal irradiance (W/m $^{2}$ )^¶
$θ_{S}$	zenith angle of the Sun
$a_{S}$	solar altitude, connected to zenith via $a_{S} = 90^{\circ} - θ_{S}$
$ϕ_{S}$	azimuth of the Sun
$A$	albedo of the ground
^¶ This parameter also can be spectral. Then, the unit would be W/(m $^{2}$ nm) and we would use the spectral irradiance $[I_{λ}] =$ W/(m $^{2}$ nm) and the spectral radiance $[L_{λ}] =$ W/(sr $\cdot$ m $^{2}$ nm)
Numerical parameters
$N_{g}$	number of points at which irradiance on the ground is evaluated in the interval $[0, d]$
$N_{m}$	number of points at which irradiance is evaluated on the module
$N_{α}$	number of directions at which light impinging onto the module is discretized

Detailed illumination model for bifacial solar cells

Abstract

1. Introduction

2. Illumination model

2.1 Irradiance components from the sky

2.2 Irradiance components from the ground

2.3 Implementation details

3. Some examples

4. Conclusions and outlook

5. Appendix. Integration of radiance across a spherical wedge

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (7)

Tables (2)

Equations (29)

Optics Express