161: Bipolar drift-diffusion with different definition of current
The problem consists of a Poisson equation for the electrostatic potential $\psi$:
\[-\nabla \lambda^2 \nabla \psi = (p-n + C)\]
and drift-diffusion equations of the transport of electrons and holes:
\[z_n \partial_t n + \nabla\cdot j_n = z_n R(n, p),\\ z_p \partial_t p + \nabla\cdot j_p = z_p R(n, p).\]
In particular, this example explores different definitions of the carrier currents and also includes the displacement flux.
module Example161_BipolarDriftDiffusionCurrentusing VoronoiFVMusing ExtendableGridsusing GridVisualizefunction main(; n = 20, # number of nodes Plotter = nothing, verbose = false, unknown_storage = :sparse, assembly = :edgewise ) ################################################################################ #### grid ################################################################################ h1 = 0.5; h2 = 4.0; h3 = 0.5 h_total = h1 + h2 + h3region numbers
region1 = 1 region2 = 2 region3 = 3 regions = [region1, region2, region3]boundary region numbers
bregion1 = 1 bregion2 = 2558 nodes; spacing ≈ 0.0126;
coord1 = collect(range(0.0; stop = h1, length = n)) coord2 = collect(range(h1; stop = h1 + h2, length = 4 * n)) coord3 = collect(range(h1 + h2; stop = h_total, length = 2 * n)) coord = glue(coord1, coord2) coord = glue(coord, coord3) grid = simplexgrid(coord) cellmask!(grid, [0.0], [h1], region1) cellmask!(grid, [h1], [h1 + h2], region2) cellmask!(grid, [h1 + h2], [h_total], region3)specify outer regions
bfacemask!(grid, [0.0], [0.0], bregion1) bfacemask!(grid, [h_total], [h_total], bregion2) ################################################################################ ######### system ################################################################################ tPrecond = 5.0 tExt = 75.0 tRamp = 1.0e-5 tEnd = tPrecond + tRamp + tExt Vprecond = 1.0 VExt = -3.0 # Define scan protocol function function scanProtocol(t) if 0.0 <= t && t <= tPrecond biasVal = Vprecond elseif tPrecond <= t && t <= tPrecond + tRamp # ramping down biasVal = Vprecond - (Vprecond - VExt) / tRamp * (t - tPrecond) elseif tPrecond + tRamp <= t && t <= tEnd biasVal = VExt else biasVal = 0.0 end return biasVal end Cn = 10 # n-doped Cp = 10 # p-doped Ca = 0.0 # intrinsic En = 1.0 # conduction band edge Ep = 0.0 # valence band edge μn = 1.0e1 # electron mobility μp = 1.0e1 # hole mobility lambda = 1.0e-1 # Debye length DirichletVal = 0.0 sys = VoronoiFVM.System(grid; unknown_storage, assembly) iphin = 1; zn = -1; enable_species!(sys, iphin, regions) iphip = 2; zp = 1; enable_species!(sys, iphip, regions) ipsi = 3; enable_species!(sys, ipsi, regions) function reaction!(f, u, node, data) nn = exp(zn * (u[iphin] - u[ipsi] + En)) np = exp(zp * (u[iphip] - u[ipsi] + Ep)) if node.region == region1 C = Cn elseif node.region == region2 C = Ca elseif node.region == region3 C = -Cp end f[ipsi] = -(C + zn * nn + zp * np) ######################## r0 = 1.0 recomb = (r0 + 1 / (nn + np)) * (nn * np * (1.0 - exp(u[iphin] - u[iphip]))) f[iphin] = zn * recomb f[iphip] = zp * recomb return nothing end function flux!(f, u, node, data) f[ipsi] = - lambda^2 * (u[ipsi, 2] - u[ipsi, 1]) ######################## bp, bm = fbernoulli_pm(-(u[ipsi, 2] - u[ipsi, 1])) nn1 = exp(zn * (u[iphin, 1] - u[ipsi, 1] + En)) np1 = exp(zp * (u[iphip, 1] - u[ipsi, 1] + Ep)) nn2 = exp(zn * (u[iphin, 2] - u[ipsi, 2] + En)) np2 = exp(zp * (u[iphip, 2] - u[ipsi, 2] + Ep)) f[iphin] = - zn * μn * (bm * nn2 - bp * nn1) f[iphip] = - zp * μp * (bp * np2 - bm * np1) return nothing end function breaction!(f, u, bnode, data) boundary_dirichlet!(f, u, bnode, iphin, bregion1, 0.0) boundary_dirichlet!(f, u, bnode, iphip, bregion1, 0.0) boundary_dirichlet!(f, u, bnode, ipsi, bregion1, 0.5 * (En + Ep) + asinh(Cn / (2 * sqrt(exp(- (En - Ep)))))) boundary_dirichlet!(f, u, bnode, iphin, bregion2, DirichletVal + scanProtocol(bnode.time)) boundary_dirichlet!(f, u, bnode, iphip, bregion2, DirichletVal + scanProtocol(bnode.time)) boundary_dirichlet!(f, u, bnode, ipsi, bregion2, 0.5 * (En + Ep) + asinh(-Cp / (2 * sqrt(exp(- (En - Ep))))) + DirichletVal + scanProtocol(bnode.time)) return nothing end function storage!(f, u, node, data) nn = exp(zn * (u[iphin] - u[ipsi] + En)) np = exp(zp * (u[iphip] - u[ipsi] + Ep)) f[iphin] = zn * nn f[iphip] = zp * np return nothing end physics!( sys, VoronoiFVM.Physics(; flux = flux!, reaction = reaction!, breaction = breaction!, storage = storage! ) ) ################################################################################ ######### time loop ################################################################################ control = SolverControl() control.Δu_opt = Inf control.max_round = 3 control.damp_initial = 0.9 control.damp_growth = 1.61 # >= 1 control.verbose = verbose control.abstol = 1.0e-5 control.reltol = 1.0e-5 control.tol_round = 1.0e-5 # Create a solution array inival2 = unknowns(sys) inival2[iphin, :] = 0.0 .+ DirichletVal ./ h_total .* coord inival2[iphip, :] = 0.0 .+ DirichletVal ./ h_total .* coord inival2[ipsi, :] = asinh(Cn / 2) .+ ((asinh(-Cp / 2) + DirichletVal) - asinh(Cn / 2)) ./ h_total .* coord solPrecond = solve(sys, inival = inival2, times = (0.0, tPrecond), control = control) control.Δt_min = 1.0e-8 control.Δt = 1.0e-8 solRamp = solve(sys, inival = solPrecond.u[end], times = (tPrecond, tPrecond + tRamp), control = control) ##### control.Δt_min = 1.0e-10 control.Δt = 1.0e-10 control.Δt_grow = 1.7 solExt = solve(sys, inival = solRamp.u[end], times = (tPrecond + tRamp, tEnd), control = control) ################################################################################ # Current calculation ################################################################################for saving Current data
I1 = zeros(0) In1 = zeros(0); Ip1 = zeros(0) I2 = zeros(0) In2 = zeros(0); Ip2 = zeros(0); Iψ2 = zeros(0) tvalues = solExt.t number_tsteps = length(tvalues) factory = TestFunctionFactory(sys) tf = testfunction(factory, [bregion1], [bregion2]) for istep in 2:number_tsteps Δt = tvalues[istep] - tvalues[istep - 1] # Time step size inival = solExt.u[istep - 1] solution = solExt.u[istep] # Variant 1: Define current as sum of charge carrier currents without electric displacement II = integrate(sys, tf, solution, inival, Δt) push!(In1, II[iphin]); push!(Ip1, II[iphip]) push!(I1, In1[istep - 1] + Ip1[istep - 1]) # Variant 2: Define all current contributions as the edge integrals and add the dielectric displacement IIEdge = VoronoiFVM.integrate_∇TxFlux(sys, tf, solution) IIEdgeOld = VoronoiFVM.integrate_∇TxFlux(sys, tf, inival) push!(In2, IIEdge[iphin]) push!(Ip2, IIEdge[iphip]) push!(Iψ2, (IIEdge[ipsi] - IIEdgeOld[ipsi]) / Δt) push!(I2, In2[istep - 1] + Ip2[istep - 1] + Iψ2[istep - 1]) end ################################################################################ ######### Plotting ################################################################################ tvalues_shift = tvalues .- (tPrecond + tRamp) if !isnothing(Plotter) vis1 = GridVisualizer(; layout = (3, 1), xlabel = "space", ylabel = "potential", legend = :lt, Plotter = Plotter, fignumber = 1) sol1 = solExt.u[1] sol2 = solExt.u[end] scalarplot!(vis1[1, 1], coord, sol1[iphin, :]; color = :darkgreen, label = "phin (start)", linewidth = 5, clear = false) scalarplot!(vis1[1, 1], coord, sol2[iphin, :]; color = :lightgreen, label = "phin (end)", linewidth = 5, clear = false) #### scalarplot!(vis1[2, 1], coord, sol1[iphip, :]; color = :darkred, label = "phip (start)", linewidth = 5, clear = false) scalarplot!(vis1[2, 1], coord, sol2[iphip, :]; color = :red, label = "phip (end)", linewidth = 5, clear = false) #### scalarplot!(vis1[3, 1], coord, sol1[ipsi, :]; color = :darkblue, label = "psi (start)", linewidth = 5, clear = false) scalarplot!(vis1[3, 1], coord, sol2[ipsi, :]; color = :lightblue, label = "psi (end)", linewidth = 5, clear = false) ################### vis2 = GridVisualizer(; layout = (3, 1), xlabel = "time", ylabel = "current", legend = :lt, xscale = :log, yscale = :log, Plotter = Plotter, fignumber = 2) scalarplot!(vis2[1, 1], tvalues_shift[2:end], abs.(In1); color = :darkgreen, label = "Jn", linewidth = 5, clear = false) scalarplot!(vis2[1, 1], tvalues_shift[2:end], abs.(Ip1); color = :darkred, label = "Jp", clear = false) scalarplot!(vis2[1, 1], tvalues_shift[2:end], abs.(I1); color = :black, linestyle = :dash, label = "Jtot", clear = false) #### scalarplot!(vis2[2, 1], tvalues_shift[2:end], abs.(In2); color = :darkgreen, label = "Jn", linewidth = 5, clear = false) scalarplot!(vis2[2, 1], tvalues_shift[2:end], abs.(Ip2); color = :darkred, label = "Jp", clear = false) scalarplot!(vis2[2, 1], tvalues_shift[2:end], abs.(Iψ2); color = :darkblue, label = "Jdisp", clear = false) scalarplot!(vis2[2, 1], tvalues_shift[2:end], abs.(I2); color = :black, linestyle = :dash, label = "Jtot", clear = false) ##### scalarplot!(vis2[3, 1], tvalues_shift[2:end], abs.(I1); color = :red, linestyle = :solid, linewidth = 5, label = "Jtot (#1)", clear = false) scalarplot!(vis2[3, 1], tvalues_shift[2:end], abs.(I2); color = :green, linestyle = :dash, label = "Jtot (#2)", clear = false) end return sum(I2)end # mainusing Testfunction runtests() testval = -965.3101329657035 for assembly in [:edgewise, :cellwise] for unknown_storage in [:dense, :sparse] @test main(; assembly, unknown_storage) ≈ testval end end return nothingendend # moduleThis page was generated using Literate.jl.