CIGS: stationary with Schottky contacts.

(source code)

Simulating stationary charge transport for CIGS with mixed Schottky/Ohmic contact conditions. Assume that SRH recombination only happens within a small regime.

module Ex108_CIGS

using ChargeTransport
using ExtendableGrids
using GridVisualize
using LaTeXStrings

# function to initialize the grid for a possible extension to other p-i-n devices.
function initialize_pin_grid(refinementfactor, h_ndoping, h_pdoping_left, h_pdoping_trap, h_pdoing_right)
    coord_ndoping = collect(range(0.0, stop = h_ndoping, length = 2 * refinementfactor))
    coord_pdoping_left = collect(range(h_ndoping, stop = (h_ndoping + h_pdoping_left), length = 3 * refinementfactor))
    coord_pdoping_plus = collect(
        range(
            (h_ndoping + h_pdoping_left),
            stop = (h_ndoping + h_pdoping_left + h_pdoping_trap),
            length = refinementfactor
        )
    )
    coord_pdoping_right = collect(
        range(
            (h_ndoping + h_pdoping_left + h_pdoping_trap),
            stop = (h_ndoping + h_pdoping_left + h_pdoping_trap + h_pdoing_right),
            length = 3 * refinementfactor
        )
    )
    coord = glue(coord_ndoping, coord_pdoping_left)
    coord = glue(coord, coord_pdoping_plus)
    coord = glue(coord, coord_pdoping_right)

    return coord
end

supported Plotters are GLMakie and PythonPlot you can set verbose also to true to display some solver information

function main(; n = 3, Plotter = nothing, verbose = false, test = false)

    ################################################################################
    if test == false
        println("Set up grid and regions")
    end
    ################################################################################

    @local_unitfactors μm cm s ns V K ps Hz W m

    constants = ChargeTransport.constants
    (; q, k_B, ε_0, Planck_constant, m_e) = constants

    eV = q * V

    # region numbers
    regionDonor = 1                           # n doped region
    regionAcceptorLeft = 2                           # p doped region
    regionAcceptorTrap = 3                           # p doped region with trap
    regionAcceptorRight = 4                           # p doped region
    regions = [regionDonor, regionAcceptorLeft, regionAcceptorTrap, regionAcceptorRight]
    numberOfRegions = length(regions)

    # boundary region numbers
    bregionDonor = 1
    bregionAcceptor = 2
    bregionDALeft = 3
    bregionALeftATrap = 4
    bregionATrapARight = 5

    # grid
    refinementfactor = 2^(n - 1)
    h_ndoping = 0.5 * μm
    h_pdoping_left = 1.0 * μm
    h_pdoping_trap = 0.1 * μm
    h_pdoing_right = 1.0 * μm
    w_device = 0.5 * μm  # width of device
    z_device = 1.0e-4 * cm  # depth of device
    h_total = h_ndoping + h_pdoping_left + h_pdoping_trap + h_pdoing_right
    coord = initialize_pin_grid(
        refinementfactor,
        h_ndoping,
        h_pdoping_left,
        h_pdoping_trap,
        h_pdoing_right
    )

    grid = simplexgrid(coord)

    # set different regions in grid
    cellmask!(grid, [0.0 * μm], [h_ndoping], regionDonor) # n doped
    cellmask!(grid, [h_ndoping], [h_ndoping + h_pdoping_left], regionAcceptorLeft) # p doped
    cellmask!(grid, [h_ndoping + h_pdoping_left], [h_ndoping + h_pdoping_left + h_pdoping_trap], regionAcceptorTrap) # p doped with traps
    cellmask!(grid, [h_ndoping + h_pdoping_left + h_pdoping_trap], [h_total], regionAcceptorRight) # p doped

    bfacemask!(grid, [h_ndoping], [h_ndoping], bregionDALeft, tol = 1.0e-18)
    bfacemask!(grid, [h_ndoping + h_pdoping_left], [h_ndoping + h_pdoping_left], bregionALeftATrap, tol = 1.0e-18)
    bfacemask!(grid, [h_ndoping + h_pdoping_left + h_pdoping_trap], [h_ndoping + h_pdoping_left + h_pdoping_trap], bregionATrapARight, tol = 1.0e-18)

TODO MO: Hier habe ich die Legende weggelassen, sonst werden die Bilder nicht richtig angezeigt bei PythonPlot

    if Plotter !== nothing
        vis = GridVisualizer(; Plotter, layout = (3, 4), size = (1550, 800))
        gridplot!(vis[1, 1], grid; Plotter, legend = :none, title = "Grid", xlabel = L"\text{space [m]}", show = true)
    end

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Define physical parameters and model")
    end
    ################################################################################

    iphin = 1 # index electron quasi Fermi potential
    iphip = 2 # index hole quasi Fermi potential
    numberOfCarriers = 2 # electrons and holes

    # physical data
    T = 300.0 * K

    # band edge energies
    Ec_ZnO = 3.4 * eV
    Ev_ZnO = 0.0 * eV

    Ec_CIGS = 3.4 * eV
    Ev_CIGS = 2.3 * eV

    EC = [Ec_ZnO, Ec_CIGS, Ec_CIGS, Ec_CIGS]
    EV = [Ev_ZnO, Ev_CIGS, Ev_CIGS, Ev_CIGS]

hole trap energy

    Et = 2.8 * eV

    # effective densities of states
    Nc = 4.35195989587969e17 / (cm^3)
    Nv = 9.139615903601645e18 / (cm^3)

    NC = [Nc, Nc, Nc, Nc]
    NV = [Nv, Nv, Nv, Nv]

    # mobilities
    mun_ZnO = 100 * (cm^2) / (V * s)
    mup_ZnO = 25 * (cm^2) / (V * s)
    mun_CIGS = 100.0 * (cm^2) / (V * s)
    mup_CIGS = 25 * (cm^2) / (V * s)

    μn = [mun_ZnO, mun_CIGS, mun_CIGS, mun_CIGS]
    μp = [mup_ZnO, mup_CIGS, mup_CIGS, mup_CIGS]

    # relative dielectric permittivity
    εr_ZnO = 9 * 1.0
    εr_CIGS = 13.6 * 1.0

    ε = [εr_ZnO, εr_CIGS, εr_CIGS, εr_CIGS]

    # recombination information parameters
    ni_ZnO = sqrt(Nc * Nv) * exp(-(Ec_ZnO - Ev_ZnO) / (2 * k_B * T))     # intrinsic concentration
    n0_ZnO = Nc * Boltzmann((Et - Ec_ZnO) / (k_B * T))                   # Boltzmann equilibrium concentration
    p0_ZnO = ni_ZnO^2 / n0_ZnO                                           # Boltzmann equilibrium concentration
    ni_CIGS = sqrt(Nc * Nv) * exp(-(Ec_CIGS - Ev_CIGS) / (2 * k_B * T))  # intrinsic concentration
    n0_CIGS = Nc * Boltzmann((Et - Ec_CIGS) / (k_B * T))                 # Boltzmann equilibrium concentration
    p0_CIGS = ni_CIGS^2 / n0_CIGS                                        # Boltzmann equilibrium concentration

    p0 = [p0_ZnO, p0_CIGS, p0_CIGS, p0_CIGS]
    n0 = [n0_ZnO, n0_CIGS, n0_CIGS, n0_CIGS]

set the lifetime value high in all other regions, such that SRH recombination can be neglected there

    SRH_LifeTime = [1.0e100, 1.0e100, 1.0e-3 * ns, 1.0e100]

    Auger = 1.0e-29 * cm^6 / s
    Radiative = 1.0e-10 * cm^3 / s

    # Schottky contact information
    An = 4 * pi * q * m_e * k_B^2 / Planck_constant^3
    Ap = 4 * pi * q * m_e * k_B^2 / Planck_constant^3
    vn = An * T^2 / (q * Nc)
    vp = Ap * T^2 / (q * Nv)
    barrier = 0.7 * eV

    # doping information
    Nd = 1.0e18 / (cm^3)
    Na = 5.5e15 / (cm^3)

    # we will impose this applied voltage on one boundary
    voltageAcceptor = 1.0 * V

    if test == false
        println("*** done\n")
    end
    ################################################################################
    if test == false
        println("Define System and fill in information about model")
    end
    ################################################################################

    # initialize Data instance and fill in data
    data = Data(grid, numberOfCarriers)
    data.modelType = Stationary
    data.F .= FermiDiracOneHalfTeSCA

    data.bulkRecombination = set_bulk_recombination(;
        iphin = iphin, iphip = iphip,
        bulk_recomb_Auger = true,
        bulk_recomb_radiative = true,
        bulk_recomb_SRH = true
    )

    data.boundaryType[bregionAcceptor] = SchottkyContact
    data.boundaryType[bregionDonor] = OhmicContact

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Define Params and fill in physical parameters")
    end
    ################################################################################

    # physical parameters
    params = Params(grid[NumCellRegions], grid[NumBFaceRegions], numberOfCarriers)
    params.temperature = T
    params.chargeNumbers[iphin] = -1
    params.chargeNumbers[iphip] = 1

    for ireg in 1:numberOfRegions           # interior region data

        params.dielectricConstant[ireg] = ε[ireg] * ε_0

        # effective DOS, band-edge energy and mobilities
        params.densityOfStates[iphin, ireg] = NC[ireg]
        params.densityOfStates[iphip, ireg] = NV[ireg]
        params.bandEdgeEnergy[iphin, ireg] = EC[ireg]
        params.bandEdgeEnergy[iphip, ireg] = EV[ireg]
        params.mobility[iphin, ireg] = μn[ireg]
        params.mobility[iphip, ireg] = μp[ireg]

        # recombination parameters
        params.recombinationRadiative[ireg] = Radiative
        params.recombinationSRHLifetime[iphin, ireg] = SRH_LifeTime[ireg]
        params.recombinationSRHLifetime[iphip, ireg] = SRH_LifeTime[ireg]
        params.recombinationSRHTrapDensity[iphin, ireg] = n0[ireg]
        params.recombinationSRHTrapDensity[iphip, ireg] = p0[ireg]
        params.recombinationAuger[iphin, ireg] = Auger
        params.recombinationAuger[iphip, ireg] = Auger

    end

    # doping -- since we do not set any doping for the traps it is automatically zero
    params.doping[iphin, regionDonor] = Nd
    params.doping[iphip, regionAcceptorLeft] = Na
    params.doping[iphip, regionAcceptorTrap] = Na
    params.doping[iphip, regionAcceptorRight] = Na

    # values for the schottky contacts
    params.SchottkyBarrier[bregionAcceptor] = barrier
    params.bVelocity[iphin, bregionAcceptor] = vn
    params.bVelocity[iphip, bregionAcceptor] = vp

    data.params = params
    ctsys = System(grid, data, unknown_storage = :sparse)

    if test == false
        show_params(ctsys)
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Define control parameters for Solver")
    end
    ################################################################################

    control = SolverControl()
    control.verbose = verbose
    control.tol_round = 1.0e-7
    control.damp_initial = 0.5
    control.damp_growth = 1.2
    control.maxiters = 30
    control.max_round = 3

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Compute solution in thermodynamic equilibrium")
    end
    ################################################################################

    # solve thermodynamic equilibrium and update initial guess
    solution = equilibrium_solve!(ctsys, control = control)
    inival = solution

    if Plotter !== nothing
        label_solution, label_density, label_energy = set_plotting_labels(data)

        # ##### set legend for plotting routines #####
        plot_energies!(vis[1, 2], ctsys, solution, "Equilibrium", label_energy)
        plot_densities!(vis[1, 3], ctsys, solution, "Equilibrium", label_density)
        plot_solution!(vis[1, 4], ctsys, solution, "Equilibrium", label_solution)
    end

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Stationary bias loop")
    end
    ################################################################################

    endVoltage = voltageAcceptor       # final bias value
    biasValues = collect(range(0, stop = endVoltage, length = 52))

    IV = zeros(0)
    chargeDensities = zeros(0)

    for i in eachindex(biasValues)

        Δu = biasValues[i] # bias

        # Apply new voltage: set non equilibrium boundary conditions
        set_contact!(ctsys, bregionAcceptor, Δu = Δu)

        if test == false
            println("bias: Δu = $(Δu) V")
        end

        # solve time step problems with timestep Δt
        solution = solve(ctsys, inival = inival, control = control)
        inival = solution

        # save IV data
        current = get_current_val(ctsys, solution)
        push!(IV, w_device * z_device * current)

        # store charge density in donor region (ZnO)
        push!(chargeDensities, charge_density(ctsys, solution)[regionDonor])


    end # bias loop

    # compute static capacitance: check this is correctly computed
    staticCapacitance = diff(chargeDensities) ./ diff(biasValues)

    # plot solution and IV curve
    if Plotter !== nothing

        plot_energies!(vis[2, 1], ctsys, solution, "bias Δu = $(endVoltage) V", label_energy)
        plot_densities!(vis[2, 2], ctsys, solution, "bias Δu = $(endVoltage) V", label_density)
        plot_solution!(vis[2, 3], ctsys, solution, "bias Δu = $(endVoltage) V", label_solution)
        plot_IV!(vis[2, 4], biasValues, IV, "bias Δu = $(biasValues[end]) V", plotGridpoints = true) # total current

        scalarplot!(
            vis[3, 1],
            biasValues[1:length(chargeDensities)],
            chargeDensities;
            color = :blue,
            markershape = :circle,
            markersize = 8,
            title = "Charge density in donor region",
            xlabel = L"\text{bias [V]}",
            ylabel = L"\text{Charge density [C]}"
        )

        scalarplot!(
            vis[3, 2],
            biasValues[1:length(staticCapacitance)],
            staticCapacitance;
            color = :blue,
            markershape = :circle,
            markersize = 8,
            title = "Static capacitance in donor region",
            xlabel = L"\text{bias [V]}",
            ylabel = L"Static capacitance [$\frac{C}{V}$]"
        )

        reveal(vis)
    end

    if test == false
        println("*** done\n")
    end

    testval = sum(filter(!isnan, solution)) / length(solution) # when using sparse storage, we get NaN values in solution
    return testval

end #  main

function test()
    testval = 1.3561479172035813

    return main(test = true) ≈ testval
end

if test == false
    println("This message should show when this module has successfully recompiled.")
end


end # module

This page was generated using Literate.jl.