PSC device with ions and different I-V scan protocols (1D).

(source code)

Simulating a three layer PSC device Ti02| MAPI | spiro-OMeTAD with mobile ions where the ion vacancy accumulation is limited by the Fermi-Dirac integral of order -1.

The time-dependent simulations are performed with abrupt interfaces. Two different I-V measurement protocols are included and the corresponding solution vectors and the I-V curve after the scan can be depicted.

module Ex103_PSC_IVMeasurement

using ChargeTransport
using ExtendableGrids
using GridVisualize
using LaTeXStrings

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,
        parameter_set = Params_PSC_TiO2_MAPI_spiro, # choose the parameter set
        otherScanProtocol = false,                  # you can choose between two scan protocols
        vacancyEnergyCalculation = true,            # assume the vacancy energy level is either given or not
    )

    @local_unitfactors μm cm s ns V K ps Hz

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

parameters

    p = parameter_set()

    # contact voltage
    voltageAcceptor = 1.2 * V

    # primary data for I-V scan protocol
    scanrate = 1.0 * V / s
    endVoltage = voltageAcceptor # bias goes until the given voltage at acceptor boundary
    tend = endVoltage / scanrate

    # Define scan protocol function
    function linearScanProtocol(t)
        return if t == Inf
            0.0
        else
            scanrate * t
        end
    end

    # Apply zero voltage on left boundary and a linear scan protocol on right boundary
    contactVoltageFunction = [zeroVoltage, linearScanProtocol]

    # Instead of a linear scan protocol, we can also apply other scan protocols which we
    # define by our own and parse to the model generator via the struct Data
    if otherScanProtocol
        # scan protocol parameter
        frequency = 10.0 * Hz
        amplitude = 0.2 * V
        tend = 1 / frequency

        # Define sinusoidal applied voltage
        function sinusoidalScanProtocol(t)
            return if t == Inf
                0.0
            else
                amplitude * sin(2.0 * pi * frequency * t)
            end
        end

        # Apply zero voltage on left boundary and a linear scan protocol on right boundary
        contactVoltageFunction = [zeroVoltage, sinusoidalScanProtocol]
    end

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

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

    δ = 4 * n        # the larger, the finer the mesh
    t = 0.5 * (cm) / δ # tolerance for geomspace and glue (with factor 10)
    k = 1.5        # the closer to 1, the closer to the boundary geomspace

    coord_n_u = collect(range(0.0, p.h_ndoping / 2, step = p.h_ndoping / (0.8 * δ)))
    coord_n_g = geomspace(
        p.h_ndoping / 2, p.h_ndoping,
        p.h_ndoping / (0.7 * δ), p.h_ndoping / (1.1 * δ),
        tol = t
    )
    coord_i_g1 = geomspace(
        p.h_ndoping, p.h_ndoping + p.h_intrinsic / k,
        p.h_intrinsic / (2.8 * δ), p.h_intrinsic / (2.1 * δ),
        tol = t
    )
    coord_i_g2 = geomspace(
        p.h_ndoping + p.h_intrinsic / k, p.h_ndoping + p.h_intrinsic,
        p.h_intrinsic / (2.1 * δ), p.h_intrinsic / (2.8 * δ),
        tol = t
    )
    coord_p_g = geomspace(
        p.h_ndoping + p.h_intrinsic, p.h_ndoping + p.h_intrinsic + p.h_pdoping / 2,
        p.h_pdoping / (1.6 * δ), p.h_pdoping / (1.6 * δ),
        tol = t
    )
    coord_p_u = collect(range(p.h_ndoping + p.h_intrinsic + p.h_pdoping / 2, p.h_ndoping + p.h_intrinsic + p.h_pdoping, step = p.h_pdoping / (1.3 * δ)))

    coord = glue(coord_n_u, coord_n_g, tol = 10 * t)
    coord = glue(coord, coord_i_g1, tol = 10 * t)
    coord = glue(coord, coord_i_g2, tol = 10 * t)
    coord = glue(coord, coord_p_g, tol = 10 * t)
    coord = glue(coord, coord_p_u, tol = 10 * t)
    grid = ExtendableGrids.simplexgrid(coord)

    # set different regions in grid
    cellmask!(grid, [0.0 * μm], [p.heightLayers[1]], p.regionDonor, tol = 1.0e-18)     # n-doped region   = 1
    cellmask!(grid, [p.heightLayers[1]], [p.heightLayers[2]], p.regionIntrinsic, tol = 1.0e-18) # intrinsic region = 2
    cellmask!(grid, [p.heightLayers[2]], [p.heightLayers[3]], p.regionAcceptor, tol = 1.0e-18)  # p-doped region   = 3

    # bfacemask! for setting different boundary regions
    bfacemask!(grid, [0.0], [0.0], p.bregionDonor, tol = 1.0e-18)     # outer left boundary
    bfacemask!(grid, [p.h_total], [p.h_total], p.bregionAcceptor, tol = 1.0e-18)  # outer right boundary
    bfacemask!(grid, [p.heightLayers[1]], [p.heightLayers[1]], p.bregionJ1, tol = 1.0e-18) # first  inner interface
    bfacemask!(grid, [p.heightLayers[2]], [p.heightLayers[2]], p.bregionJ2, tol = 1.0e-18) # second inner interface

    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 System and fill in information about model")
    end
    ################################################################################

    # Initialize Data instance and fill in predefined data
    # Currently, the way to go is to pass a contact voltage function exactly here.
    data = Data(grid, p.numberOfCarriers, contactVoltageFunction = contactVoltageFunction)

    # Possible choices: Stationary, Transient
    data.modelType = Transient

    # The default for electrons and holes is Boltzmann. Here, we set it to a more general statistics function
    data.F[p.iphin] = FermiDiracOneHalfTeSCA
    data.F[p.iphip] = FermiDiracOneHalfTeSCA

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

    # Possible choices: OhmicContact, SchottkyContact (outer boundary) and InterfaceNone,
    # InterfaceRecombination (inner boundary).
    data.boundaryType[p.bregionAcceptor] = OhmicContact
    data.boundaryType[p.bregionDonor] = OhmicContact

    # With this method, the user enable the ionic carrier parsed to ionicCarrier and gives
    # gives the information on which regions this ionic carrier is defined.
    # In this application ion vacancies only live in active perovskite layer.
    # by default the statistics function is set to FermiDiracMinusOne to limit ion depletion
    enable_ionic_carrier!(data, ionicCarrier = p.iphia, regions = [p.regionIntrinsic])

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

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

    data.params = Params(p)

    if !vacancyEnergyCalculation
        data.params.bandEdgeEnergy[p.iphia, p.regionIntrinsic] = p.Ea[p.regionIntrinsic]
    end

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

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

    if Plotter !== nothing
        ################################################################################
        println("Plot electroneutral potential, band-edge energies and doping")
        ################################################################################
        label_solution, label_density, label_energy, label_BEE = set_plotting_labels(data)

        # add labels for anion vacancy
        label_energy[1, p.iphia] = "\$E_a-q\\psi\$"; label_energy[2, p.iphia] = "\$ - q \\varphi_a\$"; label_BEE[p.iphia] = "\$E_a\$"
        label_density[p.iphia] = "\$ n_a \$";      label_solution[p.iphia] = "\$ \\varphi_a\$"

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

    control = SolverControl()
    control.verbose = verbose
    control.max_round = 5
    control.damp_initial = 0.1
    control.damp_growth = 1.21 # >= 1

    if test == false
        println("*** done\n")
    end
    ################################################################################
    if test == false
        println("Compute solution in thermodynamic equilibrium")
    end
    ################################################################################

    # calculate equilibrium solution and as initial guess
    solution = equilibrium_solve!(ctsys, control = control, vacancyEnergyCalculation = vacancyEnergyCalculation)
    inival = solution

    if Plotter !== nothing
        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("IV Measurement")
    end
    ################################################################################

    control.Δt = 0.5
    control.Δt_grow = 1.0

    if otherScanProtocol
        control.Δt_min = 1.0e-4
        control.Δt = 1.0e-4
        control.Δt_grow = 1.2
    end

    # calculation of solution
    sol = ChargeTransport.solve(ctsys, inival = inival, times = (0.0, tend), control = control)

    tvalues = sol.t
    number_tsteps = length(tvalues)

    # for saving I-V data
    IV = zeros(0)    # for IV values
    ISRHn = zeros(0)
    ISRHp = zeros(0) # for SRH recombination current
    IRadn = zeros(0)
    IRadp = zeros(0) # for radiative recombination current

    for istep in 2:number_tsteps

        Δt = tvalues[istep] - tvalues[istep - 1]  # Time step size
        inival = sol.u[istep - 1]
        solution = sol.u[istep]

        # get I-V data
        current = get_current_val(ctsys, solution, inival, Δt)

        IntSRH = integrate(ctsys, SRHRecombination!, solution)
        IntRad = integrate(ctsys, RadiativeRecombination!, solution)

        IntSRHnSum = 0.0; IntRadnSum = 0.0
        IntSRHpSum = 0.0; IntRadpSum = 0.0

        for ii in 1:p.numberOfRegions
            IntSRHnSum = IntSRHnSum - IntSRH[p.iphin, ii]
            IntRadnSum = IntRadnSum - IntRad[p.iphin, ii]

            IntSRHpSum = IntSRHpSum + IntSRH[p.iphip, ii]
            IntRadpSum = IntRadpSum + IntRad[p.iphip, ii]
        end

        push!(IV, current)
        push!(ISRHn, IntSRHnSum); push!(ISRHp, IntSRHpSum)
        push!(IRadn, IntRadnSum); push!(IRadp, IntRadpSum)

        inival = solution

    end # time loop

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

    # here in res the biasValues and the corresponding current are stored.
    # res = [biasValues IV];

    if Plotter !== nothing
        plot_energies!(vis[2, 1], ctsys, solution, "bias Δu = $(endVoltage)", label_energy)
        plot_densities!(vis[2, 2], ctsys, solution, "bias Δu = $(endVoltage)", label_density)
        plot_solution!(vis[2, 3], ctsys, solution, "bias Δu = $(endVoltage)", label_solution)
    end

    biasValues = contactVoltageFunction[p.bregionAcceptor].(tvalues)

    if Plotter !== nothing
        scalarplot!(
            vis[2, 4],
            tvalues,
            biasValues,
            markershape = :circle,
            xlabel = L"\text{time [s]}",
            ylabel = L"\text{bias [V]}",
            title = "Applied voltage over time"
        )

        ###############
        plot_IV!(vis[3, 1], biasValues[2:end], IV, "bias Δu = $(endVoltage)")

        ###############

Plot Recombination SRH recombination electrons

        scalarplot!(
            vis[3, 2],
            biasValues[2:end],
            ISRHn .* (cm^2) .* 1.0e3;
            clear = false,
            yscale = :log,
            linewidth = 5,
            color = :darkblue,
            label = "SRH recombination (n)",
            legend = :rb,
            xlabel = L"\text{bias [V]}",
            ylabel = L"current density [$\frac{\text{mA}}{\text{cm}^2}$]",
            title = "Recombination"
        )

SRH recombination holes

        scalarplot!(
            vis[3, 2],
            biasValues[2:end],
            ISRHp .* (cm^2) .* 1.0e3;
            clear = false,
            yscale = :log,
            linewidth = 5,
            color = :lightblue,
            linestyle = :dot,
            label = "SRH recombination (p)",
            legend = :rb
        )

Radiative Recombination electrons

        scalarplot!(
            vis[3, 2],
            biasValues[2:end],
            IRadn .* (cm^2) .* 1.0e3;
            clear = false,
            yscale = :log,
            linewidth = 5,
            color = :darkgreen,
            label = "Radiative recombination (n)",
            legend = :rb
        )

Radiative Recombination holes

        scalarplot!(
            vis[3, 2],
            biasValues[2:end],
            IRadp .* (cm^2) .* 1.0e3;
            clear = false,
            yscale = :log10,
            linewidth = 2,
            color = :lightgreen,
            linestyle = :dot,
            label = "Radiative recombination (p)",
            legend = :rb
        )

        reveal(vis)
    end

    if test == false
        integral = integrated_density(ctsys, sol = solution, icc = p.iphia, ireg = p.regionIntrinsic)

        println("Calculated average vacancy density is: ", integral / data.regionVolumes[p.regionIntrinsic])
        println(" ")
        if vacancyEnergyCalculation
            vacancyEnergy = data.params.bandEdgeEnergy[p.iphia, p.regionIntrinsic] / data.constants.q
            println("Value for vacancy energy is: ", vacancyEnergy, " eV. Save this value for later use.")
            println("We recommend to calculate it on a fine grid.")
            println(" ")
        end
    end

    if otherScanProtocol
        return IV[1]
    else
        return sum(IV)
    end

end #  main

function test()
    testval = 313.58311884281136; testvalOther = 0.004948787599489832
    @show main(test = true, otherScanProtocol = false)
    @show main(test = true, otherScanProtocol = true, vacancyEnergyCalculation = false)
    return main(test = true, otherScanProtocol = false) ≈ testval && abs(main(test = true, otherScanProtocol = true, vacancyEnergyCalculation = false) - testvalOther) / testvalOther < 1.0e-7
end

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

end # module

This page was generated using Literate.jl.