Calculation of absorptivity/emissivity of a silicon heterojunction cell with pyramid textures

This example calculates the absorption (and thus emissivity) in all the individual layers of a silicon heterojunction (HIT) solar cell. This paper gives more details.

Start with relevant imports:

import numpy as np
import os

from solcore.structure import Layer
from solcore import material
from solcore.light_source import LightSource
from solcore.constants import q

from rayflare.textures import regular_pyramids
from rayflare.structure import Interface, BulkLayer, Structure
from rayflare.matrix_formalism import calculate_RAT
from rayflare.matrix_formalism import process_structure
from rayflare.options import default_options

import matplotlib.pyplot as plt
import seaborn as sns

from cycler import cycler

We then generate the wavelengths (equally spaced between 300 nm and 16 microns on a log scale) at which to calculate, and set some options. Note that for this example, we are using relatively few rays and angular bins so the example runs in a reasonable amount of time, but this means the results will have a lot of noise due to the stochastic nature of the ray-tracing.

wavelengths = np.linspace(np.log(300), np.log(16*1000), 104)
wavelengths = np.round(np.floor(np.exp(wavelengths))*1e-9, 12)

options = default_options()
options.wavelengths = wavelengths
options.project_name = 'HIT_notebook'
options.n_rays = 5000
options.n_theta_bins = 20
options.nx = 5
options.ny = 5

Define the materials and the front and back material stacks. These are the layers which are deposited on the surface texture of the Si on each side. Note that these are custom materials (run the non-notebook version of this example in the examples folder to add them to your database).

Si = material('Si_UVtoMIR')()
Air = material('Air')()
ITO_front = material('ITO_front')()
ITO_back = material('ITO_back')()
Ag = material('Ag_Jiang')()
aSi_i = material('aSi_i')()
aSi_p = material('aSi_p')()
aSi_n = material('aSi_n')()


# stack based on doi:10.1038/s41563-018-0115-4
front_materials = [Layer(80e-9, ITO_front), Layer(6.5e-9, aSi_p), Layer(6.5e-9, aSi_i)]
back_materials = [Layer(6.5e-9, aSi_i), Layer(6.5e-9, aSi_p), Layer(240e-9, ITO_back)]

Define the surface textures for the front and back surface. These are pyramids in both cases, but because the pyramids are achieved using a chemical etch and are thus ‘pointing out’ of the surface in both cases. RayFlare defines ‘upright’ relative to the direction from which light is coming in, so the pyramids on the front surface are upright while those on the back surface are inverted (upright=False). The surface textures are used in combination with the layer stacks defined above to define the Interfaces of the cell, which are to be treated with the ‘RT_TMM’ method (ray-tracing + transfer-matrix method to correctly treat the surface layers). Then, finally, the bulk of the cell (Si) and thus the overall structure (front surface, bulk, back surface) is defined, also specifying the incidence medium above the cell (air) and the semi-infinite transmission medium behind the cell (silver).

surf = regular_pyramids(elevation_angle=55, upright=True)
surf_back = regular_pyramids(elevation_angle=55, upright=False)

front_surf = Interface('RT_TMM', texture=surf, layers=front_materials, name='HIT_front',
                       coherent=True)
back_surf = Interface('RT_TMM', texture=surf_back, layers=back_materials, name='HIT_back',
                      coherent=True)


bulk_Si = BulkLayer(170e-6, Si, name = 'Si_bulk') # bulk thickness in m

SC = Structure([front_surf, bulk_Si, back_surf], incidence=Air, transmission=Ag)

Now we are finished setting up, we run the calculation. process_structure calculates all the relevant matrices for the front and back surfaces (this can take a while). Then calculate_RAT does the matrix multiplication (this is fast). We then extract relevant results and use them to calculate the predicted photogenerated current in the cell.

process_structure(SC, options)

results = calculate_RAT(SC, options)

RAT = results[0]
results_per_pass = results[1]


R_per_pass = np.sum(results_per_pass['r'][0], 2)
R_0 = R_per_pass[0]
R_escape = np.sum(R_per_pass[1:, :], 0)

# only select absorbing layers, sum over passes
results_per_layer_front = np.sum(results_per_pass['a'][0], 0)

results_per_layer_back = np.sum(results_per_pass['a'][1], 0)


allres = np.flip(np.stack((results_per_layer_front[:,0],
                        results_per_layer_front[:,1]+results_per_layer_front[:,2],
                        RAT['A_bulk'][0, :],
                        results_per_layer_back[:,0] + results_per_layer_back[:,1],
                        results_per_layer_back[:,2], RAT['T'][0, :])), 0)

# calculated photogenerated current (Jsc with 100% EQE)

spectr_flux = LightSource(source_type='standard', version='AM1.5g',
                          x=wavelengths, output_units='photon_flux_per_m',
                          concentration=1).spectrum(wavelengths)[1]

Jph_Si = q * np.trapz(RAT['A_bulk'][0] * spectr_flux, wavelengths)/10 # mA/cm2

Making lookuptable for element 0 in structure
Existing lookup table found
Making lookuptable for element 2 in structure
Existing lookup table found
Ray tracing with TMM lookup table for element 0 in structure
Existing angular redistribution matrices found
Existing angular redistribution matrices found
Ray tracing with TMM lookup table for element 2 in structure
Existing angular redistribution matrices found
After iteration 1 : maximum power fraction remaining = 0.6097034485957743
After iteration 2 : maximum power fraction remaining = 0.4231168493584152
After iteration 3 : maximum power fraction remaining = 0.29274584120143654
After iteration 4 : maximum power fraction remaining = 0.20242354720194347
After iteration 5 : maximum power fraction remaining = 0.13972317102933673
After iteration 6 : maximum power fraction remaining = 0.09641891122644888
After iteration 7 : maximum power fraction remaining = 0.06653665716366666
After iteration 8 : maximum power fraction remaining = 0.04591593147741026
After iteration 9 : maximum power fraction remaining = 0.03168609279992504
After iteration 10 : maximum power fraction remaining = 0.021866277368570355
After iteration 11 : maximum power fraction remaining = 0.015089717420193792
After iteration 12 : maximum power fraction remaining = 0.01041327533750268
After iteration 13 : maximum power fraction remaining = 0.00718610544936348

Finally, plot the results:

pal = sns.cubehelix_palette(allres.shape[0], start=.5, rot=-.9)
pal.reverse()
cols = cycler('color', pal)

params = {'legend.fontsize': 'small',
          'axes.labelsize': 'small',
          'axes.titlesize': 'small',
          'xtick.labelsize': 'small',
          'ytick.labelsize': 'small',
          'axes.prop_cycle': cols}

plt.rcParams.update(params)

# plot total R, A, T
fig = plt.figure(figsize=(5,4))
ax = plt.subplot(111)
ax.semilogx(options['wavelengths']*1e6, R_escape + R_0, '--k', label=r'$R_{total}$')
ax.semilogx(options['wavelengths']*1e6, R_0, '-.k', label=r'$R_0$')
ax.stackplot(options['wavelengths']*1e6, allres,
             labels=['Ag', 'Back ITO', 'a-Si (back)', 'Bulk Si', 'a-Si (front)', 'Front ITO'])
ax.set_xlabel(r'Wavelength ($\mu$m)')
ax.set_ylabel('Absorption/Emissivity')
ax.set_xlim(min(options['wavelengths']*1e6), max(options['wavelengths']*1e6))
ax.set_ylim(0, 1)
plt.legend()
plt.show()