How do I solve equations using the Scharfetter-Gummel scheme in FiPy?

Question

I'm trying to use FiPy to simulate solar cells but I'm struggling to get reasonable results even for simple test cases.

My test problem is an abrupt 1D p-n homojunction in the dark in equilibrium. The governing system of equations are the semiconductor equations with no additional generation or recombination.

Poisson's equation determines the electric field (φ) in a semiconductor with dielectric constant, ε, given the densities of electrons (n), holes (p), donors (N_D), and acceptors (N_A), where the charge of an electron is q:

∇²φ = q(p − n + N_D − N_A) / ε

Electrons and holes drift and diffuse with current densities, J, depending on their mobilities (μ) and diffusion constants (D):

J_n = qμ_nnE + qD_n∇n

J_p = qμ_ppE − qD_p∇n

The evolution of the charge in the system is accounted for with the electron and hole continutiy equations:

∂n/∂t = (∇·J_n) / q

∂p/∂t = − (∇·J_p) / q

which can be expressed in FiPy canonical form as:

∂n/∂t = μ_n∇·(−n∇φ) + D_n∇²n

∂p/∂t = − (μ_p∇·(−p∇φ) − D_p∇²n)

To attempt to solve the problem in FiPy I first import modules and define the physical parameters.

from __future__ import print_function, division
import fipy
import numpy as np
import matplotlib.pyplot as plt

eps_0 = 8.8542e-12     # Permittivity of free space, F/m
q     = 1.6022e-19     # Charge of an electron, C
k     = 1.3807e-23     # Boltzmann constant, J/K
T     = 300            # Temperature, K
Vth   = (k*T)/q        # Thermal voltage, eV

N_ap  = 1e22           # Acceptor density in p-type layer, m^-3
N_an  = 0              # Acceptor density in n-type layer, m^-3
N_dp  = 0              # Donor density in p-type layer, m^-3
N_dn  = 1e22           # Donor density in n-type layer, m^-3

mu_n  = 1400.0e-4      # Mobilty of electrons, m^2/Vs
mu_p  = 450.0e-4       # Mobilty of holes, m^2/Vs
D_p   = k*T*mu_p/q     # Hole diffusion constant, m^2/s
D_n   = k*T*mu_n/q     # Electron diffusion constant, m^2/s
eps_r = 11.8           # Relative dielectric constant

n_i   = (5.29e19 * (T/300)**2.54 * np.exp(-6726/T))*1e6       
V_bias = 0             

Then create the mesh, solution variables, and doping profile.

nx = 20000
dx = 0.1e-9
mesh = fipy.Grid1D(dx=dx, nx=nx)
Ln = Lp = (nx/2) * dx

phi = fipy.CellVariable(mesh=mesh, hasOld=True, name='phi')
n = fipy.CellVariable(mesh=mesh, hasOld=True, name='n')
p = fipy.CellVariable(mesh=mesh, hasOld=True, name='p')
Na = fipy.CellVariable(mesh=mesh, name='Na')
Nd = fipy.CellVariable(mesh=mesh, name='Nd')

Then I set some initial values on the cell centers and impose Dirichlet boundary conditions on all parameters.

x = mesh.cellCenters[0]

n0 = n_i**2 / N_ap
nL = N_dn
p0 = N_ap
pL = n_i**2 / N_dn
phi_min = -(Vth)*np.log(p0/n_i)
phi_max = (Vth)*np.log(nL/n_i) + V_bias

Na.setValue(N_an, where=(x >= Lp))
Na.setValue(N_ap, where=(x < Lp))
Nd.setValue(N_dn, where=(x >= Lp))
Nd.setValue(N_dp, where=(x < Lp))
n.setValue(N_dn, where=(x > Lp))
n.setValue(n_i**2 / N_ap, where=(x < Lp))
p.setValue(n_i**2 / N_dn, where=(x >= Lp))
p.setValue(N_ap, where=(x < Lp))
phi.setValue((phi_max - phi_min)*x/((Ln + Lp)) + phi_min)

phi.constrain(phi_min, mesh.facesLeft)
phi.constrain(phi_max, mesh.facesRight)
n.constrain(nL, mesh.facesRight)
n.constrain(n_i**2 / p0, mesh.facesLeft)
p.constrain(n_i**2 / nL, mesh.facesRight)
p.constrain(p0, mesh.facesLeft)

I express Poisson's equation as

eps = eps_0*eps_r
rho = q * (p - n + Nd - Na)
rho.name = 'rho'

poisson = fipy.ImplicitDiffusionTerm(coeff=eps, var=phi) ==  -rho

the continuity equations as

cont_eqn_n = (fipy.TransientTerm(var=n) == 
              (fipy.ExponentialConvectionTerm(coeff=-phi.faceGrad*mu_n, var=n) 
               + fipy.ImplicitDiffusionTerm(coeff=D_n, var=n)))

cont_eqn_p = (fipy.TransientTerm(var=p) == 
              - (fipy.ExponentialConvectionTerm(coeff=-phi.faceGrad*mu_p, var=p) 
               - fipy.ImplicitDiffusionTerm(coeff=D_p, var=p)))

and solve by coupling the equations and sweeping:

eqn  = poisson & cont_eqn_n & cont_eqn_p

dt = 1e-12
steps = 50
sweeps = 10    

for step in range(steps):
    phi.updateOld()
    n.updateOld()
    p.updateOld()
    for sweep in range(sweeps):
        eqn.sweep(dt=dt)

I have played around with different values for the mesh size, time step, number of time steps, number of sweeps etc. I see some variation but haven't had any luck finding a set of conditions that give me a realistic solution. I think the problem probably lies in the expressions for the current terms.

Usually when solving these equations the current densities are are approximated using the Scharfetter-Gummel (SG) discretization scheme, rather then the direct discretization. In the SG scheme the electron current density (J) through a cell face is approximated as a function of the values of potential (φ) and charge density (n) defined on the centres of cells K and L either side as

J_n,KL=qμ_nV_T[B(δφ/V_T)n_L − B(−δφ/V_T)n_K)

where q is the charge on an electron, μ_n is the electron mobility, V_T is the thermal voltage, δφ=φ_L−φ_K, and B(x) is the Bernoulli function x/(e^x−1).

It's not obvious to me how to implement the scheme in FiPy. I have seen there is a scharfetterGummelFaceVariable but I can't work out from the documentation whether it's suitable or intended for this problem. Looking at the code it seems to only calculate the Bernoulli function multiplied by a factor e^φ_L. Is it possible to directly use the scharfetterGummelFaceVariable to solve this type of problem? If so, how? If not, is there an alternative approach that will allow me to simulate semiconductor devices using FiPy?