possibility to calculate recombination currents

examples/Ex103_PSC_IVMeasurement.jl

@@ -283,7 +283,9 @@ function main(;n = 3, Plotter = PyPlot, plotting = false, verbose = false, test
     tvalues    = range(0, stop = tend, length = number_tsteps)
 
     ## for saving I-V data
-    IV         = zeros(0) # for IV values
+    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 = 2:number_tsteps
 
@@ -302,9 +304,24 @@ function main(;n = 3, Plotter = PyPlot, plotting = false, verbose = false, test
         ## from last timestep
         solution = solve(ctsys; inival = inival, control = control, tstep = Δt)
         ## get I-V data
-        current  = get_current_val(ctsys, solution, inival, Δt)
+        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 = 1:numberOfRegions
+            IntSRHnSum = IntSRHnSum - IntSRH[iphin, ii]
+            IntRadnSum = IntRadnSum - IntRad[iphin, ii]
+
+            IntSRHpSum = IntSRHpSum + IntSRH[iphip, ii]
+            IntRadpSum = IntRadpSum + IntRad[iphip, ii]
+        end
 
         push!(IV, current)
+        push!(ISRHn, IntSRHnSum); push!(ISRHp, IntSRHpSum)
+        push!(IRadn, IntRadnSum); push!(IRadp, IntRadpSum)
 
         inival = solution
 
@@ -335,6 +352,17 @@ function main(;n = 3, Plotter = PyPlot, plotting = false, verbose = false, test
         Plotter.ylabel("bias [V]")
         Plotter.figure()
         plot_IV(Plotter, biasValues[2:end], IV, "bias \$\\Delta u\$ = $(endVoltage)")
+        ###############
+        Plotter.figure()
+        semilogy(biasValues[2:end], ISRHn.*(cm^2).*1.0e3, linewidth = 5, color = "darkblue",  label ="SRH recombination")
+        semilogy(biasValues[2:end], ISRHp.*(cm^2).*1.0e3, linewidth = 5, color = "lightblue", linestyle = ":")
+        semilogy(biasValues[2:end], IRadn.*(cm^2).*1.0e3, linewidth = 5, color = "darkgreen", label ="Radiative recombination")
+        semilogy(biasValues[2:end], IRadp.*(cm^2).*1.0e3, linewidth = 5, color = "lightgreen", linestyle = ":")
+
+        PyPlot.grid()
+        PyPlot.legend()
+        PyPlot.xlabel("bias [V]")
+        PyPlot.ylabel("current density [mAcm\$^{-2} \$]")
     end
 
     testval = sum(filter(!isnan, solution))/length(solution) # when using sparse storage, we get NaN values in solution

src/ChargeTransport.jl

@@ -66,6 +66,7 @@ export breaction!, bstorage!, reaction!, storage!, flux!
 export zeroVoltage
 export BeerLambert
 export storage!
+export reaction!, SRHRecombination!, RadiativeRecombination!, SRRecombination!, Photogeneration!
 ##################################################################
 
 include("ct_system.jl")

src/ct_physics.jl

@@ -1419,4 +1419,120 @@ function chargeCarrierFlux!(f, u, edge, data, icc, ::Type{GeneralizedSG})
 end
 
 # ##########################################################
-# ##########################################################
+# recombination kernels for calculating recombination currents
+# ##########################################################
+
+
+function SRRecombination!(f, u, bnode, data)
+
+    params          = data.params
+
+    # indices (∈ IN) of electron and hole quasi Fermi potentials specified by user (passed through recombination)
+    iphin           = data.bulkRecombination.iphin # integer index of φ_n
+    iphip           = data.bulkRecombination.iphip # integer index of φ_p
+
+    n               = get_density!(u, bnode, data, iphin)
+    p               = get_density!(u, bnode, data, iphip)
+
+    exponentialTerm = exp((q * u[iphin] - q  * u[iphip] ) / (kB * params.temperature))
+    excessDensTerm  = n * p * (1.0 - exponentialTerm)
+
+    if params.recombinationSRHvelocity[iphip, bnode.region] ≈ 0.0
+        vp =  1.0e30
+    else
+        vp = 1.0/params.recombinationSRHvelocity[iphip, bnode.region]
+    end
+
+    if params.recombinationSRHvelocity[iphin, bnode.region] ≈ 0.0
+        vn = 1.0e30
+    else
+        vn = 1.0/params.recombinationSRHvelocity[iphin, bnode.region]
+    end
+
+    kernelSRH = 1.0 / ( vp * (n + params.bRecombinationSRHTrapDensity[iphin, bnode.region]) + vn * (p + params.bRecombinationSRHTrapDensity[iphip, bnode.region] ) )
+
+    for icc ∈ data.electricCarrierList
+        icc = data.chargeCarrierList[icc]
+        f[icc] =  q * params.chargeNumbers[icc] * kernelSRH *  excessDensTerm
+    end
+
+
+end
+
+function SRHRecombination!(f, u, node, data)
+
+  params = data.params
+  ireg   = node.region
+
+  # indices (∈ IN) of electron and hole quasi Fermi potentials used by user (passed through recombination)
+  iphin  = data.bulkRecombination.iphin
+  iphip  = data.bulkRecombination.iphip
+
+  # based on user index and regularity of solution quantities or integers are used and depicted here
+  iphin  = data.chargeCarrierList[iphin]
+  iphip  = data.chargeCarrierList[iphip]
+
+  n      = get_density!(u, node, data, iphin)
+  p      = get_density!(u, node, data, iphip)
+
+  taun   = params.recombinationSRHLifetime[iphin, ireg]
+  n0     = params.recombinationSRHTrapDensity[iphin, ireg]
+  taup   = params.recombinationSRHLifetime[iphip, ireg]
+  p0     = params.recombinationSRHTrapDensity[iphip, ireg]
+
+  exponentialTerm = exp((q * u[iphin] - q * u[iphip]) / (kB * data.params.temperature))
+  excessDensTerm  = n * p * (1.0 - exponentialTerm)
+
+  kernelSRH   = params.prefactor_SRH / ( taup * (n + n0) + taun * (p + p0) )
+  ###########################################################
+  ####       right-hand side of continuity equations     ####
+  ####       for φ_n and φ_p (bipolar reaction)          ####
+  ###########################################################
+  f[iphin] = q * params.chargeNumbers[iphin] * kernelSRH * excessDensTerm
+  f[iphip] = q * params.chargeNumbers[iphip] * kernelSRH * excessDensTerm
+
+end
+
+
+function RadiativeRecombination!(f, u, node, data)
+
+  params = data.params
+  ireg   = node.region
+
+  # indices (∈ IN) of electron and hole quasi Fermi potentials used by user (passed through recombination)
+  iphin  = data.bulkRecombination.iphin
+  iphip  = data.bulkRecombination.iphip
+
+  # based on user index and regularity of solution quantities or integers are used and depicted here
+  iphin  = data.chargeCarrierList[iphin]
+  iphip  = data.chargeCarrierList[iphip]
+
+  n      = get_density!(u, node, data, iphin)
+  p      = get_density!(u, node, data, iphip)
+
+  exponentialTerm = exp((q * u[iphin] - q * u[iphip]) / (kB * data.params.temperature))
+  excessDensTerm  = n * p * (1.0 - exponentialTerm)
+
+  # calculate recombination kernel. If user adjusted Auger, radiative or SRH recombination,
+  # they are set to 0. Hence, adding them here, has no influence since we simply add by 0.0.
+  kernelRad   = params.recombinationRadiative[ireg]
+  ###########################################################
+  ####       right-hand side of continuity equations     ####
+  ####       for φ_n and φ_p (bipolar reaction)          ####
+  ###########################################################
+  f[iphin] = q * params.chargeNumbers[iphin] * kernelRad * excessDensTerm
+  f[iphip] = q * params.chargeNumbers[iphip] * kernelRad * excessDensTerm
+
+end
+
+function Photogeneration!(f, u, node, data)
+
+    generationTerm = generation(data, node, data.generationModel)
+
+    for icc ∈ data.electricCarrierList
+        icc    = data.chargeCarrierList[icc] # based on user index and regularity of solution quantities or integers are used and depicted here
+        f[icc] = q * data.params.chargeNumbers[icc] * generationTerm
+    end
+
+
+end