Introduction

Thin-film solar devices consist of short energy payback times, low active material consumption, and the possibility of low-cost production at moderate temperatures. Therefore, they will play a crucial role in a renewables-based energy economy [1]. Hence, it is not surprising that this technology is under frequent investigations by simulation, both from an electrical [2] and drift-diffusion [3] point of view.

Internally produced current needs to be guided towards the external contact point in order to extract the converted solar power. On the back side of single-junction solar devices, metallic layers are used as contact, whereas on the front side such a contact material prevents solar irradiation to reach the absorber material and instead gets reflected. Therefore, transparent conducting oxides (TCO) are used providing both a decent transparency combined with a passable conductivity. To support the TCO’s conductivity, the usage of metallization grid patterns is widespread. Due to the improved sheet resistance, grid structures allow to aim for larger cells [4]. However, because of their induced shading effect, grid structures always are a tradeoff between enlarged optical shading and improved electrical conductivity. Therefore, to find an appropriate grid design for a certain TCO layer thickness has always deserved special attention [5]. While the thickness of the TCO layer is a single optimization variable, achieving a high level of complexity in the grid design requires a large number of spatial degrees of freedom. Therefore, advanced methods like topology optimization (TO) are an interesting approach [6]. We show that TO can not only be applied to solar cells [7] or solar pin-up modules [8, 9] but also to monolithically integrated cells.

Although grid patterns and TCO layer thicknesses of real-world solar modules need to be optimized, it is widely spread that such algorithms are applied for standard test conditions (STC) [10] implying an irradiation intensity of \(1000\,\mathrm {W/m}^2\). However, this value is only the yearly peak of direct irradiation, which makes a module optimization for lower irradiations much more reasonable. Such yield-optimizations for the grid design have been performed for classical optimization algorithms [11], which differ depending on the geographical location of the solar module [12]. This work shows topology optimization applied to yield-optimizing algorithms, which leads to savings of up to 50 % in TCO and grid materials while simultaneously increasing the yearly yield by over 1 %, which is a significant improvement in terms of the large leverage effect of solar power generation.

Modeling methodology

For device simulations we combine an electrical Finite Element Method (FEM) [13,14,15] with an optical Transfer-Matrix Method (TMM) [16, 17]. For this, we divide a solar module into multiple, monolithically interconnected cells. Using a quadtree triangulation, we further divide these cells into around \(N=20000\) subdomains, where each acts as a finite element. The electrical and optical behavior of each finite element is briefly explained in the following sections.

Electrical finite element simulation

Each element k consists of a front potential \(\Phi _{\text {front}}^{k}\) and a single-diode equivalent-circuit model with a photo current \(I_{\text {ph}}^{k}\), a reverse saturation current \(I_0^{k}\), a diode ideality factor \(n_{\text {d}}^{k}\), and shunt and series resistances \(R_{\text {s}}^{k}\) and \(R_{\text {sh}}^{k}\). In each element, the net generated current is

$$\begin{aligned} I_{\text {net}}^{k}=&-f_{\text {TMM}}\varvec{\cdot }I_{\text {ph}}^{k}\nonumber \\&+ I_0^{k}\varvec{\cdot } \left( \exp \left( \frac{q_e \varvec{\cdot }\left( \Phi _{\text {front}}^{k} - I_{\text {net}}^{k}R_{\text {s}}^{k}\right) }{n_{\text {d}}^{k}k_\text{ B } T}\right) - 1\right) \nonumber \\&+ \frac{\Phi _{\text {front}}^{k} - I_{\text {net}}^{k}R_{\text {s}}^{k}}{R_{\text {sh}}^{k}}, \end{aligned}$$
(1)

with the Boltzmann constant \(k_\text{ B }\) and the linear shading \(f_{\text {TMM}}\). Equation (1) represents the raw material IV curve of the internal solar cell, which is explained in more detail in [15]. It excludes the electrical losses due to the contact structures. We use a fixed temperature of \(T = 298\,\mathrm {K}\) neglecting local module temperatures, a photo current density of \(41.17\,\mathrm {mA/cm}^2\), a reverse saturation current density of \(3\times 10^{-7}\,\mathrm {mA/cm}^2\), a diode ideality factor of 1.52, an area-normalized series resistance \(1.125\times 10^{-11}\,\Omega \mathrm {m}^2\) , and shunt resistance of \(0.067\,\Omega \mathrm {m}^2\).

All elements k are connected to their corresponding neighbors n via the total resistance \(R_{\text {front}}^{k,n} + R_{\text {front}}^{n,k}\) from the k-th to the n-th element. This leads to the following coupled system of equations

$$\begin{aligned} 0&= \sum _{n\,\in \,\mathcal {N}(k)} \underbrace{\frac{\Phi _{\text {front}}^{k} - \Phi _{\text {front}}^{n}}{R_{\text {front}}^{k,n} + R_{\text {front}}^{n,k}}}_{\text {front currents } I_{\text {front}}^{k,n}} + \, I_{\text {net}}^{k}\left( \Phi _{\text {front}}^{k}\right) , \end{aligned}$$
(2)

where \(\mathcal N(k)\) is the set of all M finite elements. In order to correctly simulate entire modules, we used periodic boundary conditions. The connecting resistances are determined from a parallel circuit from the TCO layer and an optional grid layer for each element. Measurements revealed a thickness-dependent TCO resistivity of \(8.556\times 10^{-6}\,\Omega \mathrm {m} + 7.124\times 10^{-5}\,\Omega \mathrm {m} \cdot e^{-8.631\times 10^{6} \cdot d}\) with the thickness d and a constant resistivity for the grid of \(2.7\times 10^{-8}\,\Omega \mathrm {m}\) [17].

Using the total generated current \(I_{\mathrm {out}} = \sum _{k} I_{\text {net}}^{k}\), the power conversion efficiency (PCE) can be computed via

$$\begin{aligned} \text {PCE} = \frac{P_{\text {out}}^{\text {MPP}}}{P_{\text {in}}} = \frac{V_{\text {op}}^{\text {MPP}} \varvec{\cdot } I_{\text {out}}^{\text {MPP}}}{P_{\text {in}}}, \end{aligned}$$
(3)

where \(V_{\text {op}}^{\text {MPP}}\) is the operating voltage at the maximum power point (MPP), \(P_{\text {out}}^{\text {MPP}}\) is the produced power, and \(P_{\text {in}}\) is the power of all incident photons. The resulting IV curves are fitted by a numerically robust method [18].

Optical transfer-matrix method

Our optical approach is a transfer-matrix method [19] modified with a scalar scattering theory [16] in order to be able to account for rough interfaces up to its limit of validity. The complex electric fields \({\hat{E}}_{0}\) and \({\hat{E}}_{M}\) of the first and last layer are related via

$$\begin{aligned} {\hat{E}}_0 \left( \lambda \right)&= \underbrace{ \mathbf {D}_{0}^{-1} \cdot \left[ \prod _{i=0}^{M} \mathbf {D}_i \cdot \mathbf { P}_i \cdot \mathbf {D}_{i}^{-1}\right] \cdot \mathbf { D}_{M+1} }_{=\mathbf { M}} \cdot {\hat{E}}_M \left( \lambda \right) , \end{aligned}$$
(4)

where the total transfer-matrix \(\mathbf { M}\) is calculated as a product of all propagation matrices \(\mathbf { P}_i\) and all diffraction matrices \(\mathbf { D}_i\).

Calculating yearly yield

In order to calculate a module’s yearly yield from the produced temporary power under certain conditions, yearly irradiation data are split into bins of equal irradiation power. For each bin a FEM simulation is performed and the calculated power is multiplied with the amount of hours that the module is in this bin within a year. Then the yield of each bin is summed up to get the final expected mean yield within one year. We used bins of \(10\,\mathrm {W/m}^2\) and used typical meteorological year (TMY) data averaged from 2005 to 2020 from the PVGIS-ERA5 dataset [20].

Topology optimization of metallization grid pattern

Based on the FEM calculation, different grid layouts can be constructed within a design domain \(\Omega = [0, 2.6] \times [0, 5] \,\mathrm {mm}^2\). These designs are represented by the density vector \(\underline{x} \in \{0,1\}^N\), whose components of 0 and 1 represent the potential presence of grid on the corresponding finite element. Each of those combinations represents a well-defined grid design. The function

$$\begin{aligned} f: \{0,1\}^{N} \rightarrow \mathbb R, \; \underline{x} \mapsto \mathrm {PCE}, \end{aligned}$$
(5)

maps each grid arrangement \(\underline{x}\) to a certain PCE using the FEM simulation. Topology optimization aims to find the maximum

$$\begin{aligned} \max \limits _{\underline{x} \in \{0,1\}^N} \mathrm {PCE}(\underline{x}), \end{aligned}$$
(6)

of the high-dimensional function f following the approach of [8]. Calculating a PCE is a cost of few seconds while a yield calculation takes several minutes. A topology optimization uses this functions, and therefore the mentioned factor of around 300 would be introduced by an optimization for yield instead of PCE as well. The resulting total runtime for an optimization procedure on a standard PC would increase to multiple weeks, which is far longer than practically usable. Thus, we explicitly aim for the PCE as loss function.

In order to perform a continuous optimization, the Solid Isotropic Material with Penalization (SIMP) approach with exponential functions is used [21,22,23,24]. Depending on the location of the evaluated element, the loss function can be sensitive up to one percent in PCE by reversing the binary grid density. Within the continuous optimization, we observed that the landscape of the loss function is relatively smooth but has many local maxima. For this hilly and therefore highly non-convex optimization problem, the Adaptive Moment Estimation (Adam) method [25] frequently results in the best optimum among several other optimizers like the BFGS algorithm [26,27,28,29] and multiple variations of the gradient descent method [30]. The Adam method requires a calculation of the gradient

$$\begin{aligned} \mathrm {grad}_{\underline{x}}\left( \mathrm {PCE}(\underline{x})\right), \end{aligned}$$
(7)

in each iteration step, using a numerical approximation at one point to avoid a recursive problem. The result of the optimization can be improved by a non-linear localization strategy [31]. Thereby, the finite elements are randomly divided into several batches and the associated components of the density vector are optimized one batch after another. To avoid getting stuck in a local maximum, a Gaussian blur filter is used, averaging the densities of the individual elements with their neighbors. In the optimization, the boundary conditions are taken into account so that the solar cell can already be used as a component of a solar module.

Results and discussion

In this work, we implemented cells with a layer stack of 3 mm soda–lime glass/500 nm molybdenum (Mo)/2200 nm copper indium gallium diselenide (CIGS)/50 nm cadmium sulfide (CdS)/90 nm intrinsic zinc oxide (i-ZnO)/variable thick aluminum-doped zinc oxide (AZO)/2500 nm nickel–aluminum–nickel (Ni/Al/Ni) metallization grid. Typical dimensions for monolithically integrated cells of \(2.6\,\)x\(\,5\,\mathrm {mm}^2\) are used. Other cell widths and geometries have been tested as well, but this configuration resulted in the highest performances. The size of the corresponding modules is not dependent of the cell size and therefore can be chosen arbitrarily. To model monolithcally interconnected cells, we set the boundary conditions to \(\Phi _{\text {front}}^{k} = 0\) at the left edge and \(\Phi _{\text {front}}^{k} = V_{\text {op}}\) at the right edge, where \(V_{\text {op}}\) is the total voltage across the cell.

Fig. 1
figure 1

Designs for monolithically interconnected cells. The topological grid pattern and the thickness of the TCO layer are optimized for an irradiation of \(1000\,\mathrm {W/m}^2\) (a) and \(500\,\mathrm {W/m}^2\) (b). The grid coverage and the TCO layer thickness of design a are both double the value of design b

Full size image

Figure 1 shows smallest periodic structures of two grid designs, each improved for different conditions. Both designs are the manually adapted combination of different topology optimization runs described above. They are therefore the merged design of several individual optimizations, which are not shown here for reasons of clarity. Design a) was optimized for the irradiation of \(1000\,\mathrm {W/m}^2\) distributed over the AM1.5G spectrum, whereas design b) was optimized for \(500\,\mathrm {W/m}^2\) but the same spectral distribution. Both designs are forced to have a symmetry axis in the horizontal direction. Due to their shape, they are referred as radio-antenna design and champagne-glass design in the following. The radio-antenna design is covered with \(2.89\,\%\) grid, whereas the champagne-glass design only has \(1.35\,\%\) grid. Besides the grid design, we also improved the TCO layer thicknesses via a parameter optimization for \(1000\,\mathrm {W/m}^2\) for the radio-antenna design and for \(500\,\mathrm {W/m}^2\) for the champagne-glass design, resulting in \(180\,\mathrm {nm}\) and \(90\,\mathrm {nm}\), respectively. This factor of 2 in the grid coverage as well as in the TCO layer thickness perfectly represents the irradiation difference of a factor of 2, for which both designs were optimized. Such a reduction in grid and TCO materials significantly decreases the material production costs of thin-film solar modules.

Fig. 2
figure 2

Performance of both designs from Fig. 1 as a function of irradiation intensity. Due to the smaller grid coverage and the thinner TCO layer, the champagne-glass design has a higher short circuit current density, whereas the radio-antenna design has a higher fill factor due to its increased conductivity (plot b). Both properties contribute to the PCE in (plot a), making both designs better for a different irradiation intensity. The black lines represent the possible potential of the internal semiconductor material

Full size image

In Fig. 2, the two designs are compared in their performance for different irradiation intensities. Due to its smaller grid coverage and its thinner TCO layer, the champagne-glass design consists of a higher short circuit current density independent of the given irradiation as shown in Fig. 2b. However, because of the increased conductivity of the radio-antenna design (Fig. 1a) due to its grid pattern and its thicker TCO layer, the fill factor is increased, which is particularly noticeable for high currents and hence for high irradiation intensities. For illuminations above \(870\,\mathrm {W/m}^2\), this increased electrical behavior represented by the fill factor leads to a higher PCE in Fig. 2a. For lower irradiations, however, the generated current is small enough to not require a high conductivity, which is why the higher short circuit current density outperforms the non-linear fill factor. This makes the champagne-glass design the more favorable design for small irradiations despite having a worse PCE at standard test conditions (STC) including the \(1000\,\mathrm {W/m}^2\) irradiation. For both designs, the open circuit voltage is almost at the same value for the same irradiation. For reference, the black lines represent the possible potential of the raw semiconductor material with perfect contacts.

Fig. 3
figure 3

Yearly irradiation characteristics for different locations. Hourly irradiation data have been averaged from 2005 to 2020 from the PVGIS-ERA5 dataset [20] and put into bins (gray line). The multiplication of this graph with the x-axis results in the red line, which shows how much power is generated over a year within the corresponding irradiation bin

Full size image

To test the modules on realistic data, irradiation data for multiple cities with different geometrical latitudes reaching from \(64^\circ \,\)N to \(30^\circ \,\)N have been used. The simulated modules in all locations point straight south and have an optimized angle of attack, which is \(48^\circ\) N for Reykjavik, \(42^\circ\) N for Stuttgart, \(39^\circ\) N for Barcelona, and \(30^\circ\) N for Cairo. The gray lines in Fig. 3 represent a histogram of irradiation data for each location with its given angle of attack. Multiplying this histogram data with the corresponding irradiation intensity results in the red lines. In practice, they show at which irradiation intensities the yearly yield is prevalently generated. The weighted mean of this curve indicates the virtual irradiation value, where most of the yearly power is generated.

Table 1 PCE and yield performance comparison of both designs
Full size table

Both designs have been simulated for each of the four locations. Their PCE and yield for each location are listed in Table 1. Although the PCE of the champagne-glass design is 0.51\(\%_\mathrm {rel}\) worse compared with the radio-antenna design, it outperforms it in yield for every location reaching from 0.36\(\,\%_\mathrm {rel}\) in Cairo up to 1.28\(\,\%_\mathrm {rel}\) in Reykjavik. Since the champagne-glass design favors low-light conditions, the superiority is more significant for locations with a higher amount of low-light hours, i.e., further away from the equator. For locations with extremely high irradiations such as the Lago Salar de Arizaro in north-western Argentina both designs produce the same yearly yield of \(439.0\,\mathrm {kWh}/\mathrm {m}^2\). Even at such locations most power is produced at non-peak irradiations due to low sun and clouds. This gives the champagne-glass design an advantage during these times. Therefore, for every yield optimization independent of the geographical location, a consideration of lower irradiations than \(1000\,\mathrm {W/m}^2\) is useful.

The results for PCE and yield of cross-combined geometries with one grid design and the TCO thickness of the other design and vice versa lie between those of the two standard combinations, for e.g., in the case of Stuttgart, the champagne-glass design with non-typical \(180\,\mathrm {nm}\) of TCO reveals a yield of \(245.4.\,\mathrm {kWh}/\mathrm {m}^2\), whereas the radio-antenna design with its non-typical \(90\,\mathrm {nm}\) of TCO has a yearly yield of \(245.6.\,\mathrm {kWh}/\mathrm {m}^2\). Both optimizations are therefore an improvement on the original PCE-optimized design. However, the best solution for yearly yield is the combination of both improvements, which is listed in Table 1 as champagne-glass design.

Fig. 4
figure 4

Loss mechanisms for both designs. All losses are normalized to the maximum potential of the semiconductor material (internal solar cell). While the radio-antenna design has an advantage in the electrical behavior, the champagne-glass design outperforms it in the optical characteristics

Full size image

With our finite element simulation, we can give a physical reason, why one design has a higher PCE at \(1000\,\mathrm {W/m}^2\) and the other one at \(500\,\mathrm {W/m}^2\). For this, we consider the loss analysis for both designs under both irradiations, shown in Fig. 4. For better comparison, we normalized all losses to the maximum potential of the semiconductor material (internal solar cell). In this visualization, all optical losses stay the same for different irradiations. However, losses due to grid shading and parasitic absorption within the TCO layer are much more dominant for the radio-antenna design due to its larger grid coverage and thicker TCO layer. The same reasons lead to reduced resistive losses within the front contact and front grid. For a reduced irradiation by a factor of 2, all ohmic losses in both designs consistently reduce by a factor of approximately 2. At the same time, relative losses due to local MPP mismatches (arising from a non-equipotential voltage distribution) increase for lower irradiation, since the maximum power points shift away from the internal materials MPP. Moreover, losses due to reverse currents under the grid become larger for lower irradiation intensities, since even though the photo current is reduced the backwards shunt currents stay the same and hence get more dominant in percentage.

Finally, we can conclude that the radio-antenna design deals better with the tradeoff of resistive and optical losses within the TCO and grid layer under the AM1.5G spectrum with \(1000\,\mathrm {W/m}^2\), which is a laboratory condition to measure the efficiency of solar devices. However, the champagne-glass design has its advantages in the low-light regime due to less grid shading and less parasitic absorption within the TCO layer. This low-light advantage makes its appearance in the yield calculation, since most of the time solar devices are not exposed to laboratory irradiation conditions.

Conclusion

In this work, we investigated the topology optimization of metallization grid patterns and optimized TCO layer thickness of monolithically integrated solar devices considering variable irradiation conditions. For different irradiation intensities, not only different TCO layer thicknesses but also different topological grid layouts are preferred. Since throughout a year, the irradiation is most of the time significantly below the laboratory-standardized \(1000\,\mathrm {W/m}^2\), an optimization of the grid structure and the contact layer for lower irradiations is reasonable. We have shown that despite a PCE loss at \(1000\,\mathrm {W/m}^2\) of 0.51\(\%_\mathrm {rel}\), a yield gain for different locations at various geographical latitudes from 0.36\(\,\%_\mathrm {rel}\) in Cairo up to 1.28\(\,\%_\mathrm {rel}\) in Reykjavik can be achieved. Moreover, using designs which are more efficient for lower irradiations can save up to 50 % of the TCO and grid material in production. These improvements can be achieved with comparably little effort and due to the widespread usage of solar modules have large absolute impacts. Using our electrical finite element simulation and optical transfer-matrix calculation, we are able to identify all losses being responsible for a better yield performance.