Solve the scaled stationary bipolar DD equation system using Gummel algorithm.
[n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_gummel_map (x, D, Na, Nd,
pin, nin, Vin, Fnin,
Fpin, l2, er, u0n,
uminn, vsatn, betan,
Nrefn, u0p, uminp, vsatp,
betap, Nrefp, theta, tn, tp,
Cn, Cp, an, ap, Ecritnin, Ecritpin,
toll, maxit, ptoll, pmaxit)
input:
x spatial grid
D, Na, Nd doping profile
pin initial guess for hole concentration
nin initial guess for electron concentration
Vin initial guess for electrostatic potential
Fnin initial guess for electron Fermi potential
Fpin initial guess for hole Fermi potential
l2 scaled Debye length squared
er relative electric permittivity
u0n, uminn, vsatn, Nrefn electron mobility model coefficients
u0p, uminp, vsatp, Nrefp hole mobility model coefficients
theta intrinsic carrier density
tn, tp, Cn, Cp,
an, ap,
Ecritnin, Ecritpin generation recombination model parameters
toll tolerance for Gummel iterarion convergence test
maxit maximum number of Gummel iterarions
ptoll convergence test tolerance for the non linear
Poisson solver
pmaxit maximum number of Newton iterarions
output:
n electron concentration
p hole concentration
V electrostatic potential
Fn electron Fermi potential
Fp hole Fermi potential
Jn electron current density
Jp hole current density
it number of Gummel iterations performed
res total potential increment at each step
The following code
% physical constants and parameters
secs1d_physical_constants;
secs1d_silicon_material_properties;
% geometry
L = 10e-6; % [m]
xm = L/2;
Nelements = 1000;
x = linspace (0, L, Nelements+1)';
sinodes = [1:length(x)];
% dielectric constant (silicon)
er = esir * ones (Nelements, 1);
% doping profile [m^{-3}]
Na = 1e23 * (x <= xm);
Nd = 1e23 * (x > xm);
% avoid zero doping
D = Nd - Na;
% initial guess for n, p, V, phin, phip
V_p = -1;
V_n = 0;
Fp = V_p * (x <= xm);
Fn = Fp;
p = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni./abs(D)) .^2)) .* (x <= xm) + ...
ni^2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x > xm);
n = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^ 2)) .* (x > xm) + ...
ni ^ 2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x <= xm);
V = Fn + Vth * log (n / ni);
% scaling factors
xbar = L; % [m]
nbar = norm(D, 'inf'); % [m^{-3}]
Vbar = Vth; % [V]
mubar = max (u0n, u0p); % [m^2 V^{-1} s^{-1}]
tbar = xbar^2 / (mubar * Vbar); % [s]
Rbar = nbar / tbar; % [m^{-3} s^{-1}]
Ebar = Vbar / xbar; % [V m^{-1}]
Jbar = q * mubar * nbar * Ebar; % [A m^{-2}]
CAubar = Rbar / nbar^3; % [m^6 s^{-1}]
abar = 1/xbar; % [m^{-1}]
% scaling procedure
l2 = e0 * Vbar / (q * nbar * xbar^2);
theta = ni / nbar;
xin = x / xbar;
Din = D / nbar;
Nain = Na / nbar;
Ndin = Nd / nbar;
pin = p / nbar;
nin = n / nbar;
Vin = V / Vbar;
Fnin = Vin - log (nin);
Fpin = Vin + log (pin);
tnin = tn / tbar;
tpin = tp / tbar;
u0nin = u0n / mubar;
uminnin = uminn / mubar;
vsatnin = vsatn / (mubar * Ebar);
u0pin = u0p / mubar;
uminpin = uminp / mubar;
vsatpin = vsatp / (mubar * Ebar);
Nrefnin = Nrefn / nbar;
Nrefpin = Nrefp / nbar;
Cnin = Cn / CAubar;
Cpin = Cp / CAubar;
anin = an / abar;
apin = ap / abar;
Ecritnin = Ecritn / Ebar;
Ecritpin = Ecritp / Ebar;
% tolerances for convergence checks
toll = 1e-3;
maxit = 1000;
ptoll = 1e-12;
pmaxit = 1000;
% solve the problem using the full DD model
[nout, pout, Vout, Fnout, Fpout, Jnout, Jpout, it, res] = ...
secs1d_dd_gummel_map (xin, Din, Nain, Ndin, pin, nin, Vin, Fnin, Fpin, ...
l2, er, u0nin, uminnin, vsatnin, betan, Nrefnin, ...
u0pin, uminpin, vsatpin, betap, Nrefpin, theta, ...
tnin, tpin, Cnin, Cpin, anin, apin, ...
Ecritnin, Ecritpin, toll, maxit, ptoll, pmaxit);
% Descaling procedure
n = nout*nbar;
p = pout*nbar;
V = Vout*Vbar;
Fn = V - Vth*log(n/ni);
Fp = V + Vth*log(p/ni);
dV = diff(V);
dx = diff(x);
E = -dV./dx;
% band structure
Efn = -Fn;
Efp = -Fp;
Ec = Vth*log(Nc./n)+Efn;
Ev = -Vth*log(Nv./p)+Efp;
plot (x, Efn, x, Efp, x, Ec, x, Ev)
legend ('Efn', 'Efp', 'Ec', 'Ev')
axis tight
Produces the following figure
| Figure 1 |
|---|
![]() |
Package: secs1d