A Tool for a Fast and Accurate Evaluation of the Energy Production of Bifacial PV Modules

1. Introduction

The global energy landscape is rapidly shifting towards renewable energy sources to alleviate the climate change and reduce the dependence on fossil fuels. Among the renewable sources, photovoltaic (PV) technology stands out for its scalability, decreasing costs, and widespread adoption.

Traditionally, the market has been dominated by monofacial PV modules (hereinafter also referred to as panels), which capture sunlight only on one side (front). However, in recent times bifacial PV modules, engineered to capture sunlight also on the backside (henceforth denoted as rear) have emerged as a valuable alternative, as they are supposed to offer higher energy yields at a price of a minor and affordable upgrade of the fabrication equipment; promising results have been obtained in many applications like agrivoltaics, aquavoltaics, or in high-reflection environments such as deserts, snow-covered areas, and urban architectures with reflective surfaces (e.g., [1,2,3,4,5,6,7,8,9,10,11] and references therein).

The literature is populated by papers dealing with the comparison between monofacial and bifacial panels in terms of cost, efficiency, and production under various operating conditions [12,13,14,15,16,17,18]. However, there is still a lack of unambiguous guidelines defining in which circumstances (in terms of orientation, inclination, reflection from ground, partial shading on the front or rear, and so on) bifacial modules offer clear and distinct advantages with respect to the monofacial counterparts. This challenge can only be addressed resorting to advanced simulation tools suited to predict the performance of bifacial and monofacial modules in any environment, thereby supporting engineers and researchers to make choices oriented to energy yield and return-on-investment maximization.

In this paper, we present a comprehensive simulation tool conceived and designed for bifacial PV modules, a simplified variant of which is dedicated to the monofacial ones. The tool allows for a fast, yet accurate enough, evaluation of the energy performance of the modules, it is user-friendly and totally immune to convergence problems. The information to be provided are the day of the year, the coordinates of the geographical site where the module is installed, the evolution along the day of total and diffuse irradiance on the horizontal plane, as well as of ambient temperature, in that site (available on online databases), orientation, tilt, albedo, and some datasheet parameters. The tool determines the I–V characteristic of the module, and thus the maximum produced power, at each clock time in the day. The computation is automated and very fast, taking only tens of seconds to cover a whole day with a 15-minute discretization.

The tool is entirely developed as a Matlab code, except when there are discrepancies among the cells in terms of irradiances (on the front and/or on the rear), temperature, or other parameters. In such a nonuniform case, Matlab is still used to determine front/rear irradiances and temperature on the cells behaving uniformly, but the circuit simulation program OrCAD PSPICE [19] is invoked to compute the I–V characteristics along the day with a high-granularity cell-level approach. By virtue of this feature, it is possible to predict the impact of partial architectural shading, bird drops, snow, and defects (e.g., cracks).

The reminder of this work is articulated as follows. In Section 2, the details of the tool are given. A clear and extensive overview of the adopted models is provided using a tutorial style, with the aim to inspire and ease the reader’s own implementation. In Section 3, the tool is adopted to explore the daily, monthly, and yearly energy production of bifacial PV modules located in Naples by varying orientation and tilt. The bifacial performance is compared to that of a south-oriented monofacial module tilted by 30°, which is the well-established optimum condition in Naples and is therefore taken as a reference. It is found that west- or east-oriented vertical bifacial modules allow matching or even exceeding (depending on the albedo) the yearly-produced energy of the reference monofacial one; this suggests their adoption for e.g., agrivoltaics or noise barriers along highways. Moreover, a counter-intuitive result is also achieved: it is shown that the reduction of the current photogenerated of one cell on the rear (e.g., due to a mud spot) can lead to bypass activation even if the Sun rays are only hitting the front. Lastly, the detrimental impact of cracked cells is quantified. Conclusions are finally drawn in Section 4.

2. The Tool

The block diagram of the tool is represented in Figure 1. As can be seen, the tool is articulated into four blocks. Blocks #1, #2, and #3 are Matlab codes, while for #4 two options are possible: (1) if all the cells in the module share the same irradiance on the front, the same irradiance on the rear, the same operating temperature, and thus the same photogenerated current on the front and the same on the rear, as well as the same parameters, a Matlab block is used to solve the single-diode model (SDM) described in Section 2.4 at a module level, that is, the voltage drop across the module is obtained by multiplying the voltage drop across the cell V_cell by the number of cells embedded in the panel; (2) if there is a nonuniformity over the cells in the module (either on the front or on the rear), in terms of e.g. irradiances, electrical parameters (dictated by e.g., architectural shading, defects) then the solution is demanded to the circuit simulator PSPICE, in which each cell is individually described by its own subcircuit. PSPICE was chosen since it is equipped with a robust engine suited for the solution of systems of nonlinear algebraic equations, which ensures short CPU times and unlikely occurrence of convergence problems.

The detailed description of each block is provided in the following.

2.1. Block #1

Block #1 is fed with

• The latitude ϕ and the local longitude λ, denoted as λ_local, of the geographical site where the PV module is installed (ϕ is equal to 0° at the equator, and ranges from 0° to 90° to the north, from 0° to -90° to the south; λ_local is equal to 0° at the Prime, or Greenwich, meridian, and ranges from 0° to 180° to east, from 0° to -180° to west).
• The longitude λ_standard of the standard meridian on which the clock time (CKT, also referred to as watch time or standard local time) is based.
• The day of the year n.
• The CKT values during the day, i.e., the daytime discretization.

and calculates

• The solar declination δ, i.e., the angle between the Sun rays (the beam radiation) and the equatorial plane, which is positively defined in the northern hemisphere. Angle δ can be reasonably assumed constant during a day (it varies by at most 0.5°), and dependent on the day of the year n through the empirical Cooper’s relation

$$\delta =23.45\cdot \sin \left[{\displaystyle \frac{360}{365}}\cdot \left(284+n\right)\right]$$

(1)

On the summer solstice, δ assumes the maximum value 23.45° (23°27’); on the winter solstice it is equal to the minimum value -23.45° (-23°-27’); δ=0° on the spring and autumn equinoxes.

• The so-called True Local Time (TLT), which can be obtained from the CKT by applying two corrections. The first correction is dictated by the difference between the local longitude λ_local and the longitude of the standard meridian λ_standard; more specifically, the displacement of 1° between these longitudes corresponds to 4 minutes. The second correction is made to account for the non-constancy of the rotation rate of the Earth around the Sun during the year. This effect can be described by introducing a characteristic time referred to as Equation of Time (Eqt) expressed in minutes, which depends on n through the following relation [20,21,22]:

$$Eqt=K_1\cdot \left(K_2+K_3\cos B-K_4\sin B-K_5\cos 2B-K_6\sin 2B\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B=\left(n-1\right)\cdot {\displaystyle \frac{360}{365}}$$

(2)

with K₁=229.2 minutes, K₂=0.000075, K₃=0.001868, K₄=0.032077, K₅=0.014615, K₆=0.04089 [21,22] (or 0.040849 [20] without appreciable Eqt variation). By applying both corrections, under standard time/conditions the relation between TLT (in hours) and CKT is given by

$$TLT=CKT-{\displaystyle \frac{4}{60}}\cdot \left({\lambda }_{standard}-{\lambda }_{local}\right)+{\displaystyle \frac{Eqt}{60}}$$

(3)

which, under daylight-saving time/conditions, is modified into

$$TLT=CKT-1-{\displaystyle \frac{4}{60}}\cdot \left({\lambda }_{standard}-{\lambda }_{local}\right)+{\displaystyle \frac{Eqt}{60}}$$

(4)

that is, a further hour must be deducted.

• The hour angle ω, i.e., the angular displacement of the Sun with respect to the local meridian compared to the case in which TLT=12 due to the rotation from west to east of the Earth around its axis (also denoted as angle subtended by the Sun [23]). Angle ω is negative in the morning, positive in the afternoon, and given by [20,22,24,25,26,27]

$$\omega =15\cdot \left(TLT-12\right)$$

(5)

As can be seen, a one-hour deviation corresponds to a displacement of 15°.

• The solar altitude, or elevation, α, i.e., the angle between the horizontal plane (also denoted as plane of horizon) and the Sun rays; α is a function of latitude ϕ, solar declination δ, and hour angle ω according to [20,21,23,24,25,26]

$$\sin \alpha =\sin \varphi \cdot \sin \delta +\cos \varphi \cdot \cos \delta \cdot \cos \omega$$

(6)

The outcomes of block #1 for some application examples are reported in Table 1.

2.2. Block #2

Block #2 is fed with

• The latitude ϕ of the site.
• The CKT values.
• The azimuth angle of the module front γ, which defines the module orientation, as it is the angular displacement from south of the projection of the normal to the module front onto the horizontal plane. An azimuth γ=0° means that the front is south-oriented; γ is positive clockwise (to west), reaching 180° next to north, and is negative counter-clockwise (to east), reaching -180° next to north [21]. Hence, γ=90° corresponds to a module with west-oriented front, while γ=-90° identifies an east-oriented front.
• The tilt angle of the module front β, which is the inclination with respect to the horizontal plane.
• The solar declination δ and the hour angle ω computed by block #1.

and calculates

• The azimuth angle of the module rear γ_rear=γ-180°.
• The tilt angle of the module rear as the complement of β, that is, β_rear=180°-β.
• The incidence angles of the Sun rays on the module front θ and rear θ_rear (both being ≥0). θ is the angle between the Sun rays and the normal to the module front, and is evaluated from [21,22,23,27]

$$\begin{array}{c}\cos \theta =\cos \delta \cdot \cos \omega \cdot \left(\sin \varphi \cdot \cos \gamma \cdot \sin \beta +\cos \varphi \cdot \cos \beta \right)+\\ +\sin \delta \cdot \left(\sin \varphi \cdot \cos \beta -\cos \varphi \cdot \cos \gamma \cdot \sin \beta \right)+\\ +\cos \delta \cdot \sin \omega \cdot \sin \gamma \cdot \sin \beta \end{array}$$

(7)

It can be inferred that θ=0° for a horizontal module (β=0°) at the equator (ϕ=0°) at TLT=12 (noon, ω=0°) of the equinoxes (δ=0°); angle θ may exceed 90°, which means that the Sun is behind the surface.

As an example, if it is required to determine the angle of incidence of the Sun rays on the front of a module installed at Madison, WI, USA, with orientation γ=15° and tilt β=45° at CKT=10.7 (10:42 AM) on February 13 (n=44), then block #1 evaluates that δ=-13.95°, Eqt=-14.26 minutes, TLT=10.50 (10:30 AM), ω=-22.50°, α=29.38°, and block #2 calculates θ=35° [21].

It is noteworthy that the θ formulation provided in [20] turned out to be wrong.

θ_rear is the angle between the Sun rays and the normal to the module back, and is computed as 180°-θ or equivalently from

$$\begin{array}{c}\cos {\theta }_{rear}=\cos \delta \cdot \cos \omega \cdot \left(\sin \varphi \cdot \cos {\gamma }_{rear}\cdot \sin {\beta }_{rear}+\cos \varphi \cdot \cos {\beta }_{rear}\right)+\\ +\sin \delta \cdot \left(\sin \varphi \cdot \cos {\beta }_{rear}-\cos \varphi \cdot \cos {\gamma }_{rear}\cdot \sin {\beta }_{rear}\right)+\\ +\cos \delta \cdot \sin \omega \cdot \sin {\gamma }_{rear}\cdot \sin {\beta }_{rear}\end{array}$$

(8)

Figure 2 is intended to provide a clear understanding some key angles introduced so far.

2.3. Block #3

The inputs of block #3 are

• The solar altitude α evaluated by block #1.
• The azimuth angles γ, γ_rear, the tilt angles β, β_rear, and the incidence angles θ, θ_rear.
• The total irradiance G_toth, the diffuse irradiance G_dh hitting the horizontal plane (the beam, or direct, irradiance G_bh is determined as G_toth-G_dh), and the ambient temperature T_amb vs. CKT at the selected geographical site. For the analysis performed in Section 3, these data were taken from the PhotoVoltaic Geographical Information System (PVGIS) website [28]. Here it is stated that they were evaluated for the mean day of the chosen month from satellite data through a sophisticate algorithm accounting for sky obstruction (shading) by local terrain features (hills or mountains) calculated from a digital elevation model.
• The albedo value, namely, the ratio of reflected upward radiation from the ground to the incident downward radiation upon it (typical albedo values are 0.04 for fresh asphalt, 0.1-0.15 for soil ground, 0.25-0.3 for green grass, 0.4 for desert sand, 0.55 for fresh concrete, and 0.8 for freshly fallen snow).
• Some key parameters available in the datasheet, i.e., the temperature T_NOCT, the short circuit current I_scnom, and the percentage temperature coefficient TCI_sc of the short-circuit current I_sc, the definitions of which will be provided in the following.
• The block calculates
• The beam irradiance impinging on the module front as [21,23,29,30,31]

$$G_b=G_{bh}\cdot {\displaystyle \frac{\cos \theta }{\sin \alpha }}$$

(9)

and on the rear

$$G_{b,rear}=G_{bh}\cdot {\displaystyle \frac{\cos {\theta }_{rear}}{\sin \alpha }}=-G_{bh}\cdot {\displaystyle \frac{\cos \theta }{\sin \alpha }}$$

(10)

• The diffuse irradiance on the front coming from the sky, for which there are two options.

If the sky is completely and densely overcast, the same irradiance is coming from any point of the sky, which is thus called isotropic. In this case, the diffuse irradiance is expressed as

$$G_{di}=G_{dh}\cdot F_i$$

(11)

where the dimensionless front-sky view factor F_i (<1) accounts for the reduction of the sky dome due to the tilt angle β; using the cross-string approach [32], it can be determined that

$$F_i={\displaystyle \frac{1+\cos \beta }{2}}={\cos }^2{\displaystyle \frac{\beta }{2}}$$

(12)

Formulation (12), sometimes referred to as Kondrat’yev’s view factor [33], is widely accepted [2,9,12,20,21,23,24,25,29,34,35,36,37,38,39,40,41,42], and coincides with the ratio between the projection on the horizontal plane of the portion of hemisphere seen from the titled panel and the projection of the whole sky dome [43].

Analogously, for the module rear

$$G_{di,rear}=G_{dh}\cdot F_{i,rear}$$

(13)

with the rear-sky view factor [2,40,42]

$$F_{i,rear}={\displaystyle \frac{1+\cos {\beta }_{rear}}{2}}={\cos }^2{\displaystyle \frac{{\beta }_{rear}}{2}}={\displaystyle \frac{1-\cos \beta }{2}}={\sin }^2{\displaystyle \frac{\beta }{2}}$$

(14)

If the sky is clear or at least partially cloudy, the diffuse irradiance depends on the position of the Sun in the sky due to effects like the horizon brightening and circumsolar radiation, and the sky is thus referred to as anisotropic. In this case, the diffuse irradiance can be expressed as

$$G_{da}=G_{dh}\cdot F_a$$

(15)

where the front-sky view factor F_a is given by the sum of two terms F_a1 and F_a2 accounting for the horizon brightening and circumsolar radiation, respectively [21,29,31]

$$F_a=F_{a1}+F_{a2}=F_i\cdot \left(1-A_i\right)\cdot \left(1+f\cdot {\sin }^3{\displaystyle \frac{\beta }{2}}\right)+A_i\cdot {\displaystyle \frac{\cos \theta }{\sin \alpha }}$$

(16)

In (16),

$$f=\sqrt{{\displaystyle \frac{G_{bh}}{G_{bh}+G_{dh}}}}$$

(17)

and A_i is the anisotropic index, given by

$$A_i={\displaystyle \frac{G_{bh}}{G_0}}$$

(18)

where G₀ is the solar irradiance incident on a horizontal plane outside the atmosphere, also referred to as extraterrestrial irradiance on a horizontal surface, whose expression is

$$G_0=G_{sc}\cdot \left[1+0.033\cdot \cos \left({\displaystyle \frac{360}{365}}\cdot n\right)\right]\cdot \left(\cos \varphi \cdot \cos \delta \cdot \cos \omega +\sin \varphi \cdot \sin \delta \right)$$

(19)

G_sc being the solar constant, equal to 1353 W/m².

Differently from the isotropic view factor F_i, the anisotropic one F_a can be >1 if the panel is oriented to the portion of the sky dome where the Sun is located. It is worth noting that F_a reduces to F_i if G_bh=0 W/m² (uniformly cloudy, or isotropic, sky), which implies that f=0 and A_i=0.

For the rear of the module,

$$G_{da,rear}=G_{dh}\cdot F_{a,rear}$$

(20)

where

$$F_{a,rear}=F_{a1,rear}+F_{a2,rear}=F_{i,rear}\cdot \left(1-A_i\right)\cdot \left(1+f\cdot {\sin }^3{\displaystyle \frac{{\beta }_{rear}}{2}}\right)+A_i\cdot {\displaystyle \frac{\cos {\theta }_{rear}}{\sin \alpha }}$$

(21)

• The diffuse irradiance on the front due to the reflection from the ground, typically considered as a Lambertian (isotropic) process, calculated as

$$G_{d,albedo}=G_{toth}\cdot albedo\cdot F_{albedo}$$

(22)

where the front-ground view factor F_albedo is expressed as [2,9,10,12,13,20,21,29,30,31,34,35,36,37,38,40,42,44,45]

$$F_{albedo}={\displaystyle \frac{1-\cos \beta }{2}}={\sin }^2{\displaystyle \frac{\beta }{2}}$$

(23)

Model (23). correctly works in the entire range of practical β values, namely, from β=0° (horizontal panel) to β=90° (vertically deployed panel); F_albedo=0 if β=0° since the front does not see the ground; F_albedo=0.5 if β=90° since the front sees half ground.

By referring to the module rear,

$$G_{d,albedo,rear}=G_{toth}\cdot albedo\cdot F_{albedo,rear}$$

(24)

For the rear-ground factor F_albedo,rear we examined the feasibility of the analogous of (23), namely [39,40,41,42]

$$F_{albedo,rear}={\displaystyle \frac{1-\cos {\beta }_{rear}}{2}}={\sin }^2{\displaystyle \frac{{\beta }_{rear}}{2}}={\displaystyle \frac{1+\cos \beta }{2}}={\cos }^2{\displaystyle \frac{\beta }{2}}$$

(25)

Unfortunately, different from (23), model (25) does not provide meaningful results under all the possible configurations. Let us define as d the vertical distance between the lowest point of the module and the ground. If d=0 and β=0° the panel backside lies on the ground, which would imply G_{d,albedo,rear}=0, but (25) predicts F_albedo,rear=1, which is the consequence of its view factor nature. To tackle this issue, we propose an alternative, still simple enough, model for factor F_albedo,rear to also account for the dependence on d, which is given by

$$F_{albedo,rear}={\displaystyle \frac{1-g_{rear}\left(d\right)\cdot \cos {\beta }_{rear}}{2}}={\displaystyle \frac{1+g_{rear}\left(d\right)\cdot \cos \beta }{2}}$$

(26)

where

$$g_{rear}\left(d\right)=1-2\cdot \exp \left(-{\displaystyle \frac{2d}{H_{panel}}}\right)$$

(27)

H_panel being the module height. The behavior of F_albedo,rear expressed by (26), (27) as a function of the ratio d/H_panel for various tilt angles β is reported in Figure 3.

If d=0, then g_rear=-1 and (26) reduces to

$$F_{albedo,rear}={\displaystyle \frac{1+\cos {\beta }_{rear}}{2}}={\displaystyle \frac{1-\cos \beta }{2}}$$

(28)

Which shows that if β=0°, then F_albedo,rear=0, i.e., there is no reflection contribution since the panel lies on the ground, while if β=90°, then F_albedo,rear=0.5 as the module backside sees half ground.

For β<90°, F_albedo,rear increases with d, as the module rear sees more and more ground. If d is significantly higher than H_panel, then then g_rear=1 and (26) saturates to (25), from which, if β=0°, then F_albedo,rear=1 as in this case the rear of the horizontal suspended panel sees all the ground.

If β=90°, then F_albedo,rear=0.5 independently of d, since it can be assumed that the rear sees half ground regardless of the distance from panel to ground.

This simple model given by (24) with (26), (27) inherently assumes that the ground is completely unshaded. If a large portion of the ground is shaded (e.g., in the case of a large PV field composed of horizontal bifacial modules suspended on a framework), (24) can be replaced by

$$G_{d,albedo,rear}=G_{dh}\cdot albedo\cdot F_{albedo,rear}$$

(29)

Accounting for a limited shadow on an irradiated ground, as the one produced by the panel itself, would markedly increase the complexity of the model; examples are the approaches in [12,13,40] under the assumption of infinite-length modules, and the complex procedure in [46].

• The total irradiances on the front (G) and rear (G_rear). The total irradiance on the front is calculated as

$$G=G_b+G_{di}+G_{d,albedo}$$

(30)

for an isotropic sky, or

$$G=G_b+G_{da}+G_{d,albedo}$$

(31)

for an anisotropic one.

The total irradiance on the rear is evaluated as

$$G_{rear}=G_{b,rear}+G_{di,rear}+G_{d,albedo,rear}$$

(32)

for an isotropic sky, and

$$G_{rear}=G_{b,rear}+G_{da,rear}+G_{d,albedo,rear}$$

(33)

for an anisotropic one.

• The operating temperature T of the cells as a function of CKT. By disregarding the self-heating, which is a reasonable approximation when the cells are producing power, T [°C] is a function of T_amb and G according to the following linear law [6,12,26,47,48]:

$$T=T_{amb}+\left(T_{NOCT}-20°\right)\cdot {\displaystyle \frac{G}{G_{nomT}}}$$

(34)

T_NOCT being the normal (nominal) operating cell temperature (NOCT) measured under open-circuit conditions at T_amb=20°C, nominal irradiance G_nomT=800 W/m² impinging on the front, AM 1.5, and wind speed lower than 1 m/s. As mentioned earlier, T_amb vs. CKT is available from PVGIS, and the total irradiance on the panel front G was previously determined by the tool. Formulation (34) is also referred to as NOCT model.

• One the irradiance G and the temperature T vs. CKT are known, the current photogenerated by the module front is evaluated for each CKT as [23,49]

$$I_{ph}=I_{scnom}\cdot {\displaystyle \frac{G}{G_{nom}}}+\kappa \cdot \left(T-25°\right)$$

(35)

where the datasheet parameter I_scnom is the short-circuit current measured by keeping the frontside under standard test conditions (STCs), i.e., nominal irradiance G_nom=1000 W/m², 1.5 Air Mass, and cell temperature equal to 25°C, and covering the backside or at least markedly limiting its irradiance. It is worth noting that the first term on the RHS is the photogenerated current that would be obtained for an irradiance G at T=25°C. The temperature coefficient κ [A/°C] can be easily determined from the datasheet parameter TCI_sc [%/°C], namely, the percentage temperature coefficient of I_sc measured by the module manufacturer under the same conditions mentioned above; TCI_sc is given by

$$TC{I}_{sc}={\displaystyle \frac{100}{I_{scnom}}}\cdot {\displaystyle \frac{\partial I_{sc}}{\partial T}}={\displaystyle \frac{100}{I_{scnom}}}\cdot \kappa$$

(36)

Equation (35) allows accounting for the positive temperature coefficient of I_ph due to the bandgap shrinking and the resulting increase in the number of photons with enough energy to generate electron-hole pairs, and can be recast as

$$I_{ph}=I_{scnom}\cdot {\displaystyle \frac{G}{G_{nom}}}+\kappa \cdot \left(T-T_0+T_0-25°\right)=I_{scnom}\cdot {\displaystyle \frac{G}{G_{nom}}}+\kappa \cdot \left(\Delta T+T_0-25°\right)$$

(37)

where T₀ is a reference temperature and ΔT=T-T₀; by assuming T₀=27°C [39,50]

$$I_{ph}=I_{scnom}\cdot {\displaystyle \frac{G}{G_{nom}}}+\kappa \cdot \left(\Delta T+2°\right)$$

(38)

The current photogenerated by the module rear vs. CKT is calculated as

$$I_{ph,rear}=ERF\cdot \left[I_{scnom}\cdot {\displaystyle \frac{G_{rear}}{G_{nom}}}+\kappa \cdot \left(T-25°\right)\right]=ERF\cdot \left[I_{scnom}\cdot {\displaystyle \frac{G_{rear}}{G_{nom}}}+\kappa \cdot \left(\Delta T+2°\right)\right]$$

(39)

where ERF (<1) is an efficiency reduction factor needed to account for the fact that the PV cell is asymmetric, being technologically designed to maximize light absorption on the front (The parameter usually provided in the datasheet of the module is the so-called bifaciality factor, namely, the ratio between the power produced by the rear illuminated under STCs while the front is covered, and the power produced by the front illuminated under STCs while the rear is covered), and T is given by (34), that is, it assumed to coincide with the temperature due to the irradiance G impinging on the panel front.

2.4. Block #4

Block #4 is intended to evaluate the I–V characteristic of the module (and consequently the maximum produced power P_max) at each CKT during the day.

The individual cell is described through an extended version of the SDM presented in [39,41,48,50] and inspired by [51]. In particular, the cell current I_cell is expressed as

$$I_{cell}=I_{ph}+I_{ph,rear}-I_D-I_{sh}+I_{av}$$

(40)

I_ph and I_ph,rear being given by (38) and (39). In (40), I_D is the current flowing through the intrinsic resistance-free diode. The expression of I_D is derived from the well-known Shockley’s equation

$$I_D\left(T\right)=I_0\left(T\right)\cdot \left[\exp \left({\displaystyle \frac{V_D}{\eta \cdot V_T}}\right)-1\right]$$

(41)

where V_D is the voltage drop across the diode, given by

$$V_D=V_{cell}+R_s\cdot I_{cell}$$

(42)

V_cell and R_s being the voltage drop over the cell and the series resistance, the dimensionless parameter η is the ideality coefficient, and $V_T={\displaystyle \frac{k}{q}}\cdot \left(T+273\right)$ is the thermal voltage, k and q being the Boltzmann’s constant and the elementary charge. As described in detail in [41,50,52], (41) can be conveniently modified into

$$I_D\left(\Delta T\right)=I_0\left(T_0\right)\cdot \left[\exp \left({\displaystyle \frac{V_D+{\varphi }_0\cdot \Delta T}{\eta V_{T0}+\eta \frac{k}{q}\cdot \Delta T}}\right)-1\right]$$

(43)

which better evidences the diode current dependence on the temperature rise ΔT=T-T₀. The value of the key coefficient ϕ₀ is extracted through the following procedure. As illustrated in Figure 1, block #4 also receives the open-circuit voltage V_ocnom of the module front and the temperature coefficient TCV_oc of the open-circuit voltage V_oc, both measured by keeping the frontside under STCs and the backside covered; TCV_oc is given by

$$TC{V}_{oc}={\displaystyle \frac{100}{V_{ocnom}}}\cdot {\displaystyle \frac{\partial V_{oc}}{\partial T}}$$

(44)

from which ${\displaystyle \frac{\partial V_{oc}}{\partial T}}$ is determined. Then block #4 executes Matlab simulations of the module with covered backside (setting I_ph,rear=0 A) with G=G_nom=1000 W/m² and T=25°C and 30°C and automatically calculates ${\displaystyle \frac{\partial V_{oc}}{\partial T}}$; parameter ϕ₀ is calibrated to obtain a good agreement between the ${\displaystyle \frac{\partial V_{oc}}{\partial T}}$ value resulting from simulations and the one calculated from (44).

I_sh is the current traversing the shunt resistance R_sh, given by

$$I_{sh}={\displaystyle \frac{V_D}{R_{sh}}}$$

(45)

Resistance R_sh is assumed to be temperature-insensitive; nevertheless, R_sh might exhibit a negative temperature coefficient induced by the thermal trapping-detrapping of carriers through the defect states in the space-charge region [53].

I_av is the (positive) avalanche current induced by a reiterate impact-ionization mechanism (avalanche multiplication) in the space-charge region of the cell when subject to high reverse voltages, e.g., in the presence of localized shading; I_av is expressed as

$$I_{av}=-{\displaystyle \frac{V_D}{R_{sh}}}\cdot {\displaystyle \frac{a_{II}}{{\left(1-\frac{V_D}{BV}\right)}^{m_{II}}}}$$

(46)

where BV (<0 V) is the breakdown voltage of the junction, while a_II and m_II are dimensionless fitting parameters. From (46) it can be inferred that, if V_D is positive, then the denominator of (46) is high and I_av is negligible; instead, if V_D is negative and approaches the negative BV, then the positive I_av increases and tends to infinity. The temperature-related mitigation of avalanche multiplication can be enabled through the relation

$$BV\left(T\right)=BV\left(T_0\right)\cdot \exp \left(b_{II}\cdot \Delta T\right)$$

(47)

b_II [°C^-1] (>0) being another fitting parameter.

For the series resistance R_s, a positive temperature coefficient is accounted for according to the following power relation:

$$R_s\left(T\right)=R_s\left(T_0\right)\cdot {\left({\displaystyle \frac{T+273}{T_0+273}}\right)}^{m_R}=R_s\left(T_0\right)\cdot {\left({\displaystyle \frac{T_0+273+\Delta T}{T_0+273}}\right)}^{m_R}$$

(48)

with m_R>0, although in principle any other law (e.g., exponential [54] or linear) can be implemented.

It was mentioned that if all the cells in the module share the same G, G_rear, T, and thus the same I_ph, I_ph,rear, as well as the same parameters, the tool resorts to Matlab to solve the SDM given by (40) with (42), (43), (45), (46), (48) at a module level: once the I_cell–V_cell characteristic is evaluated at a selected CKT, the I–V curve of the module is obtained by considering I=I_cell and V=N∙I_cell, N being the number of cells in the module.

Instead, if any nonuniformity is present, the commercial circuit simulator OrCAD PSPICE is enabled to compute the I–V characteristic by extending the approach proposed in [39,41,48,50,55] for monofacial panels to the bifacial counterparts. Such an approach can be described as follows.

• The module is composed by N series-connected cells, each modeled with a subcircuit implementing the SDM described above, which is fed with G, G_rear, ΔT.
• The module can be partitioned into a chosen number of subpanels, each equipped with a bypass diode.
• The PSPICE temperature of all components embedded in the circuit is forced to the reference value T₀=27°C; the temperature rise ΔT, represented as a voltage, is provided to analog behavioral modeling, or ABM, components (nonlinear current/voltage sources) to modify the temperature-sensitive parameters.

An illustrative sketch of the subcircuit modeling the individual cell is shown in Figure 4.

The ABM components A₁ and A₂ calculate the current photogenerated on the module front using (38) and (39), respectively. The standard diode, being at the reference temperature T₀, conducts the current

$$I_D\left(T_0\right)=I_0\left(T_0\right)\cdot \left[\exp \left({\displaystyle \frac{V_D}{\eta V_{T0}}}\right)-1\right]$$

(49)

V_T0 being equal to ${\displaystyle \frac{k}{q}}\cdot \left(T_0+273\right)$. In order to allow the flow of the temperature-dependent current I_D(ΔT) given by (43) in the diode branch, the ABM B is used, which computes and forces the current I_D(ΔT)-I_D(T₀), where I_D(ΔT) is given by (43).

To account for the avalanche multiplication, the cell current I_cell is computed as the sum of an avalanche-free current (i.e., the current that would flow in the absence of avalanche), and an avalanche-dictated current I_av forced by the ABM C and given by (46).

The power dependence on temperature of the series resistance R_s is implemented by forcing on the branch traversed by I_cell a voltage drop given by the product by I_cell and R_s, where R_s is expressed by (48). This is put into practice by making use of the ABM designated as D in Figure 4, which accepts as inputs the temperature rise ΔT and the current I_cell (transformed into a voltage) and imposes the drop I_cell∙R_s(ΔT).

2.5. Simplified variant of the tool for monofacial modules

A simplified tool variant is also available, which allows simulating monofacial modules. In this variant, block #1 coincides with the one of the tool for bifacial panels, block #2 only determines the incidence angle θ, block #3 evaluates G, T, I_ph, and block #4 solves the SDM only with the source for the current I_ph photogenerated by the front.

3. Results and Discussion

3.1. Optimization of Orientation and Tilt for a Bifacial Module

This Section demonstrates the effectiveness of the proposed tool in optimizing the orientation and tilt of a bifacial module at a chosen geographical site.

The site selected as a case-study is Naples (ϕ=40°50’, λ_local=14°15’). It is widely accepted that for the northern hemisphere the yearly energy production for a monofacial panel is maximized if the panel is oriented to south (γ=0°) and the tilt angle is determined by the empirical rule of thumb β=ϕ-15° to ϕ-10°, which allows the best performance (the minimum average θ value) from the spring equinox to the autumn one, i.e., when there is the maximum energy availability [21,22,26,56,57] (A few papers report different findings, an example being [36], where it is stated that the best choice (in terms of daily global solar radiation) for ϕ falling in the range from 36° to 46° is β=ϕ-26° to ϕ-28° from April to September, while along the whole year is β=ϕ-6° to ϕ-12°.). Hence, for the specific case of Naples, γ and β were chosen equal to 0° and 30° (≈ϕ-10°), respectively. Such a condition, schematically shown in Figure 5a, will be henceforth considered as a reference to assess the performance of bifacial modules in the same technology.

The analysis was performed by considering the monofacial ET Solar silicon module ET-M54050 embedding N=40 cells and subdivided in two 20-cell subpanels, each equipped with a bypass diode [39,41,48,50]. The key datasheet parameters are peak power of 50 Wp, I_scnom=2.81 A, T_NOCT=44.4°C, TCI_sc=0.06%/°C, TCV_oc=-0.397%/°C, while the other (either calculated or optimized) model parameters are κ=1.69 mV/°C, I₀(T₀)=8 nA, η=1.2, ϕ₀=4.4 mV/°C, BV=-15 V, R_sh=100 Ω, R_s(T₀)=12 mΩ, m_R=1.5, a_II=0.1, m_II=1.1. The extension of the cell model to the bifacial case was made by adding in parallel to the current source I_ph another current source forcing a photogenerated current I_ph,rear given by (39), where ERF was chosen equal to 0.9, that is, the efficiency on the rear was assumed to be 90% of that on the front side.

The investigation was carried out for d=0 (i) by initially assuming that the cells share the same irradiance and electrical parameters, (ii) by considering clear-sky anisotropic conditions; (iii) by neglecting the reflection from the ground (albedo=0) for both the monofacial panel and for the front of the bifacial panel. Instead, the reflection from the ground was taken into account for the rear of the bifacial panel, according to (33), G_{d,albedo,rear} being given by (26), (27).

Under the above assumptions, the tool was used to evaluate the I–V characteristic of both the monofacial and the corresponding bifacial panel as a function of CKT for a given day of the year n; consequently, the maximum power produced P_max was determined for each CKT, and the daily produced energy was calculated by integrating P_max over the whole day. The monthly energy was obtained by multiplying the energy produced in a representative day of the selected month (the 15th in our analysis) by the number of days in that month; for instance, the energy produced in April was determined multiplying by 30 the energy produced on April 15. This simplified approach allowed speeding up the analysis without significantly sacrificing the accuracy of the results.

Figure 6 shows the energy production along one year normalized to the peak power (50 Wp) for the reference monofacial module (south oriented, i.e., with γ=0° and tilted by β=30°) and for the corresponding (γ=0°, β=30°) bifacial module, the latter for various albedo values. It can be inferred that the energy increment obtained with a bifacial module under this condition is marginal due to the poor irradiance contribution of the rear.

In Figure 7, both the monofacial and bifacial panels are south-oriented (γ=0°) and vertically deployed (β=90°); the vertical arrangement draws interest as more bifacial modules can be installed in the same field. In this case, using a monofacial module is not convenient since from April to September the angle of incidence is on average rather high, and therefore the energy production is low. Differently from the β=30° condition, here using a bifacial panel allows obtaining a significant gain, as the backside is also hit by beam irradiance.

Figure 8 and Figure 9 report the case of west- (γ=90°) and east-oriented (γ=-90°) vertical (β=90°) modules, respectively. Again, the monofacial module is compared to the corresponding bifacial one for different albedo values. The following considerations are in order:

• By specifically referring to the monofacial module, the energy produced by orienting the frontside to west is slightly better than the one obtained by orienting it to east; this result is reasonable since Naples faces the sea to the west while mountains lie to the east.
• As a main finding, it is observed that in these cases the bifacial module allows achieving a significant improvement. While the monofacial panel receives beam irradiance for only half of the day, bifacial panels benefit from effective beam irradiance (hitting the module sides with low incidence angles) both in the morning (rear for a west-oriented panel, as sketched in Figure 5b, and front for an east-oriented one) and afternoon (the other way around).
• West- and east-oriented vertical bifacial panels produce the same amount of energy.

In Figure 10, the normalized energies produced by the vertical bifacial panels for a south-, west-, and east-oriented front (the albedo being 0.2) are compared to the reference monofacial case. By inspecting the curves, it can be inferred that:

• West- and east-oriented vertical bifacial modules (which produce the same energy) offer equal or even improved performance with respect to the reference monofacial counterpart during the time span from April to September, in which they benefit from a low incidence angle of the Sun rays hitting the front and rear of the panel in the mid-morning and the mid-afternoon [1,10]. In terms of yearly energies (all reported in Table 2), the gain compared to the reference monofacial case amounts to 2.5% and 15% for albedo values equal to 0.2 and 0.5, respectively.
• West- and east-oriented vertical bifacial modules also provide a considerable production improvement with respect to the south-oriented vertical bifacial counterpart, which is estimated to be 13-15%, regardless of the albedo. This is again due to the much better performance of west- and east-oriented panels from April to September, while during wintertime the south-oriented module performs better. An in-depth insight into this behavior can be achieved by showing the normalized maximum power over CKT on July 15 (Figure 11a) and December 15 (Figure 11b) for such cases. In July (and more in general during the whole period from late spring to early autumn), it is confirmed that the orientations of the module front to west and east allow for a very effective exploitation of the beam light impinging on one of the two sides in the mid-morning and on the other in the mid-afternoon, whereas in December (and more in general during wintertime) the orientation to south is better over the mid-day. Along the whole year, the first effect markedly prevails over the second.

Lastly, west- and east-oriented vertical bifacial modules are compared to the south-oriented bifacial one tilted by the optimum β=30°. From Table 2, it can be inferred that for albedo=0.2 the latter performs slightly better, while the vertical ones ensure a superior production for higher albedo values (Such findings align fairly well with those derived in [2], where a simplified analysis only based on the evaluation of annual irradiance [kWh/m²] incident on single (i.e., not included in rows) bifacial PV modules was performed by considering Tel Aviv (ϕ=32°04’, λ_local=34°46’) as a geographical site. West- (γ=90°) and east-oriented (γ=-90°) vertical (β=90°) bifacial modules are compared with a south-oriented (γ=0°) bifacial one tilted by β=20°, optimum angle given by the empirical rule ϕ-15° to ϕ-10°. It was found that (i) the annual incident irradiance is significantly higher for the latter panel, although its rear side suffers from a negligible irradiance contribution; (ii) vertical bifacial PV modules can produce more annual energy than a bifacial one with an optimum tilt angle if the ground albedo is very high.).

As key conclusion of this analysis is that west- and east-oriented vertical bifacial panels can be recommended for some applications in Naples (and more in general in Italy), e.g., for agrivoltaics, as the vertical installation does not significantly reduce the land available for cultivation, or for building noise barriers along highways.

3.2. Nonuniform Irradiance Distribution over the Panel Rear

Let us consider the scenario of a vertical (β=90°) bifacial module with east-oriented front (γ=-90°) located in Naples for agrivoltaics applications (albedo=0.2, d=0, ERF=0.9) on July 15. Under clear-sky conditions, the outcomes of block #3 are:

• G=370.49 W/m², T=38.55°C, I_ph=1.064 A, G_rear=153.33 W/m², I_ph,rear=0.408 A at 11:00 AM;
• G=88.23 W/m², T=31.5°C, I_ph=0.259 A, G_rear=436.45 W/m², I_ph,rear=1.114 A at 3:00 PM.

The normalized power vs. voltage characteristics of the module under uniform irradiance conditions are reported in Figure 12a (CKT=11:00 AM) and 12b (CKT=3:00 PM); the maximum powers normalized to the peak value provided in the datasheet are 0.504 and 0.488, respectively. A likely case of a mud spot partially covering the rear of cell #10, sketched in Figure 13, was then considered, which was assumed to reduce the irradiance down to 10% of that of the clean cells. As can be seen, the PSPICE simulations lead to a counter-intuitive result, that is, although the spot is located on the rear, submodule #1 embedding cell #10 is bypassed also in the morning (CKT=11:00 AM) when the Sun rays hit the module front; such unexpected detrimental behavior should also be considered when comparing west- or east-oriented vertical bifacial panels to the reference monofacial counterpart. On the other hand, as expected, the power production dramatically reduces in the afternoon (CKT=3:00 PM), when the Sun rays impinge on the backside.

3.3. Cracked Cell(s)

The proposed tool is also suited to describe a condition with malfunctioning or defected cells. Let us again consider the scenario of a vertical (β=90°) bifacial module with east-oriented front (γ=-90°) located in Naples for agrivoltaics applications (albedo=0.2, d=0, ERF=0.9) on July 15. Under clear-sky conditions, the outcomes of block #3 at 10:00 AM are: G=463.22 W/m², T=40.29°C, I_ph=1.327 A, G_rear=143.19 W/m², I_ph,rear=0.385 A. Figure 14 illustrates the normalized power against voltage for the panel under normal conditions, i.e., with all non-defected cells. Then the PSPICE-based tool version was used to assess the impact of cracked cells [58], which were emulated by setting the shunt resistance R_sh to 1 Ω instead of 100 Ω. More specifically, the cases of one, two, and three cracked cells belonging to submodule #1 were simulated, and the resulting power–voltage curves were reported in Figure 14 as well. It can be observed that the yield degradation is perceptible only for a significant number of damaged cells.

4. Conclusions

Bifacial PV modules are relentlessly drawing attention by virtue of their feature of collecting sunlight on both sides. In this paper, we have proposed an in-house tool specifically conceived to accurately simulate such modules in a short time and without occurrence of convergence problems. Once the day of the year and the geographical site are selected, the tool is suited to account for all the operating and environmental conditions, as well as for a nonuniformity condition across the cells induced by localized shading, dirt, bird drops, or defects. The tool is user-friendly, poorly resource-demanding, and very fast: the calculation of the I–V characteristics along an entire day with a 15-minute discretization takes only tens of seconds. A simulation campaign has been performed with the aim to optimize orientation and tilt of bifacial modules in Naples, chosen as a case-study. The normalized energies produced by such modules under clear-sky conditions have been compared to the reference one obtained with a monofacial panel in the same technology oriented to south and tilted by 30°, which is assumed to be the optimized configuration in Naples. As a main finding, it has been determined that either west- or east-oriented vertical bifacial modules allow equalizing or, for medium/high albedos, even improving the performance with respect to the reference case. This is mainly due to the low incidence angle of the Sun rays impinging on the front and rear of the bifacial module in the mid-morning and mid-afternoon during the time frame from late spring to early autumn. Thanks to the cell-level granularity ensured by the circuital PSPICE-based version of the tool, the impact of dirt obscuring a cell on the rear has been assessed; it has been surprisingly shown that the submodule embedding that cell is driven into bypass even when the Sun rays are hitting the front. Lastly, the influence of cracked cells has been examined. By virtue of such promising results, we believe that our tool is a good candidate to provide trenchant support for researchers and engineers, as it ensures that no critical factors affecting the bifacial module performance are overlooked.