Peter -
The semiconductor drift-diffusion equations are notoriously challenging to
solve. As you may well be aware, many people use alternative formulations such
as pseudo-Fermi-levels and Slotboom variables to mitigate some of the
difficulties that arise from the enormous range of concentrations. The
Scharfetter-Gummel discretization is also widely applied to improve the
convergence of the combined drift and diffusion terms.
The fact that these flux peaks scale as dx**2 is consistent with this being a
discretization error. Newton iterations would probably help drive down the
error.
I'm not sure what was intended by terms like
`(-phi.arithmeticFaceValue).divergence`. Fluxes (and currents) are naturally
described at the face centers between cells. As such, I'd define these like:
>>> J_n = -mu_n * n.harmonicFaceValue * phi.faceGrad + D_n * n.faceGrad
There is, unfortunately, no easy way to plot the magnitude of this flux. In the
attached, I manually use Matplotlib to do it.
- Jon
> On Nov 30, 2018, at 10:25 AM, Peter Abele <ppab...@web.de> wrote:
>
> Hello,
>
> I am using fipy to simulate semiconductor pn junctions, see the attached
> python script "01_p_n.py".
>
> Calculation of the electrical field, the width of the space charge region,
> and the potential in the semiconductor looks fine even with biasing.
> Evaluating the current densities (drift and diffusion) for the electrons and
> holes results in extreme peaks independent of the applied biasing, see
> "Current_Density.png".
> This graph was calculated with no bias voltage.
> In steady state current density should be constant along the semiconductor.
> Outside the space charge region the current density looks ok. I even get the
> exponential current density increase with applied voltage.
>
> Observation:
> It seems that the current density peaks coincide with the second derivative
> of the electron and/or hole concentration.
> When I additionally multiply the second derivative of the electron
> concentration with the electron concentration, I get nearly the same shape as
> for the current density!
> >>> viewer_n = fipy.Viewer(vars=(n.grad[0].grad[0]*n), title="2.
> >>> derivative n * n", legend='lower left')
>
> An other observation is that the peak current densities (j_peak) depends on
> the mesh size (dx), see "eval.eps".
> It seems that j_peak ~ dx**2!
>
> Do I see the limit of the numerical calculation due to the extreme change in
> electron / hole concentration at the border of the space charge region?
> Is the big difference in the cell variables (potential (-10 ... 1) and
> hole/electron concentration (0 ... 1E24)) a problem?
> Any proposal what I can investigate to get a constant current density along
> the semiconductor?
>
>
> Thanks and best regards
> Peter
>
>
> <eval.eps><Current_Density.png><01_p_n.py>_______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
from __future__ import print_function, division
import fipy
import numpy as np
from matplotlib import 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 = 1e24 # 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 = 1e24 # 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 = 1*k*T*mu_p/q # Hole diffusion constant, m^2/s
D_n = 1*k*T*mu_n/q # Electron diffusion constant, m^2/s
eps_r = 11.8 # Relative dielectric constant
n_i_Si = 10E15
# Grenze bei 1E11 !!
tau = 1E-3
eps = eps_0*eps_r
# ****************
# generate mesh
# ****************
Lx = 150E-9
dx = .15E-9
nx = int(Lx/dx)
mesh = fipy.Grid1D(dx=dx, nx=nx)
Lcenter = Lx / 2
n = fipy.CellVariable(mesh=mesh, hasOld=True, name='n')
p = fipy.CellVariable(mesh=mesh, hasOld=True, name='p')
phi = fipy.CellVariable(mesh=mesh, hasOld=True, name='phi')
rho = fipy.CellVariable(mesh=mesh, name='rho')
Nd = fipy.CellVariable(mesh=mesh, name='Nd')
Na = fipy.CellVariable(mesh=mesh, name='Na')
x = mesh.cellCenters[0]
# ******************************
# ******************************
# Doping profile
# ******************************
# ******************************
Na.setValue((+N_ap + np.sqrt(N_ap**2 + 4*n_i_Si**2))/2, where=(x<=Lcenter))
p.setValue( (+N_ap + np.sqrt(N_ap**2 + 4*n_i_Si**2))/2, where=(x<=Lcenter))
n.setValue( (-N_ap + np.sqrt(N_ap**2 + 4*n_i_Si**2))/2, where=(x<=Lcenter))
Nd.setValue((+N_dn + np.sqrt(N_dn**2 + 4*n_i_Si**2))/2, where=(x>Lcenter))
n.setValue( (+N_dn + np.sqrt(N_dn**2 + 4*n_i_Si**2))/2, where=(x>Lcenter))
p.setValue( (-N_dn + np.sqrt(N_dn**2 + 4*n_i_Si**2))/2, where=(x>Lcenter))
equ_cont_n = fipy.TransientTerm(coeff=1, var=n) == (
- 0.0 * (fipy.ImplicitSourceTerm(coeff=p,var=n) - n_i_Si**2)/tau
+ fipy.ConvectionTerm(coeff=mu_n * -phi.faceGrad, var=n)
+ fipy.DiffusionTerm(coeff=D_n, var=n)
)
equ_cont_p = fipy.TransientTerm(coeff=1, var=p) == (
- 0.0 * (fipy.ImplicitSourceTerm(coeff=n,var=p) - n_i_Si**2)/tau
- fipy.ConvectionTerm(coeff=mu_p * -phi.faceGrad, var=p)
+ fipy.DiffusionTerm(coeff=D_p, var=p)
)
equ_poisson = fipy.ImplicitDiffusionTerm(coeff=-eps, var=phi) == rho
# ***************************************
# ***************************************
# Ohmic contacts boundary condition
# ***************************************
# ***************************************
# ***************
# Potential
# ***************
#phi.setValue(0)
phi.faceGrad.constrain(0, mesh.facesRight)
phi.faceGrad.constrain(0, mesh.facesLeft)
# *******************************************
#
# electrons and holes at the Ohmic contacts
#
# *******************************************
# **************
# p-contact
# **************
p.constrain((+N_ap + np.sqrt(N_ap**2 + 4*n_i_Si**2))/2, mesh.facesLeft)
n.constrain((-N_ap + np.sqrt(N_ap**2 + 4*n_i_Si**2))/2, mesh.facesLeft)
#n.constrain(1E22, mesh.facesLeft)
# **************
# n-contact
# **************
n.constrain((+N_dn + np.sqrt(N_dn**2 + 4*n_i_Si**2))/2, mesh.facesRight)
p.constrain((-N_dn + np.sqrt(N_dn**2 + 4*n_i_Si**2))/2, mesh.facesRight)
viewer_np = fipy.Viewer(vars=(n,p), log=True)
viewer_rho = fipy.Viewer(vars=rho, limits={'ymin': -180000, 'ymax': 180000})
viewer_phi = fipy.Viewer(vars=phi)
viewer_EField = fipy.Viewer(vars=(-phi.arithmeticFaceValue).divergence, limits={'ymin': -50E6, 'ymax': 0}, title="E-Field")
J_n_drift = -mu_n * n.harmonicFaceValue * phi.faceGrad
J_n_drift.name = r"$J_{n,drift}$"
J_n_diff = D_n * n.faceGrad
J_n_diff.name = r"$J_{n,diff}$"
J_n = -mu_n * n.harmonicFaceValue * phi.faceGrad + D_n * n.faceGrad
J_n.name = r"$J_n$"
J_p = mu_p * p.harmonicFaceValue * phi.faceGrad + D_p * p.faceGrad
J_p.name = r"$J_p$"
J_tot = J_n + J_p
J_tot.name = r"$J_{tot}$"
# viewer_j_n = fipy.Viewer(vars=J_n[0], legend='lower left')
#
# viewer_j_n_drift = fipy.Viewer(vars=J_n_drift[0], legend='lower left')
#
# viewer_j_n_diff = fipy.Viewer(vars=J_n_diff[0], legend='lower left')
#
# viewer_j_tot = fipy.Viewer(vars=J_tot[0], legend='lower left')
fig = plt.figure(10)
ax_j = fig.add_subplot(111)
ax_j.plot(mesh.faceCenters[0].value, J_n[0].value, label="$J_n$")
ax_j.plot(mesh.faceCenters[0].value, J_p[0].value, label="$J_p$")
ax_j.plot(mesh.faceCenters[0].value, J_tot[0].value, label="$J_{tot}$")
ax_j.legend(loc="upper right")
dt = 0.9 * dx**2 / (2 * D_n)
dt = 1E-15
steps =1000
sweeps = 8
raw_input("Exit?")
bias_points = [0.0]
for bias in bias_points:
# ***************
# Potential
# ***************
V_bi_p = k*T/q*np.arcsinh(-N_ap/(2*n_i_Si))
V_bi_n = k*T/q*np.arcsinh(+N_dn/(2*n_i_Si)) + bias
phi.constrain(V_bi_p, mesh.facesLeft)
phi.constrain(V_bi_n, mesh.facesRight)
for step in range(steps):
rho.setValue(q * (p - n + Nd - Na))
n.updateOld()
p.updateOld()
phi.updateOld()
#rho.updateOld()
if (step%10==0):
viewer_np.plot()
viewer_rho.plot()
viewer_phi.plot()
viewer_EField.plot()
ax_j.cla()
ax_j.plot(mesh.faceCenters[0].value, J_n[0].value, label="$J_n$")
ax_j.plot(mesh.faceCenters[0].value, J_p[0].value, label="$J_p$")
ax_j.plot(mesh.faceCenters[0].value, J_tot[0].value, label="$J_{tot}$")
ax_j.legend(loc="upper right")
# viewer_j_n.plot()
# viewer_j_n_drift.plot()
# viewer_j_n_diff.plot()
# viewer_j_tot.plot()
print("Step: "+ str(step))
res0 = res1 =res2 = 1E100
for sweep in range(sweeps):
res0 = equ_poisson.sweep(var=phi, dt=dt)
rho.setValue(q * (p - n + Nd - Na))
res1 = equ_cont_n.sweep(var=n, dt=dt)
rho.setValue(q * (p - n + Nd - Na))
res2 = equ_cont_p.sweep(var=p, dt=dt)
rho.setValue(q * (p - n + Nd - Na))
print("Res: " + str(res0) + " " + str(res1) + " " + str(res2))
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
