1%% Copyright (C) 2004-2012 Carlo de Falco 2%% 3%% This file is part of 4%% SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator 5%% 6%% SECS1D is free software; you can redistribute it and/or modify 7%% it under the terms of the GNU General Public License as published by 8%% the Free Software Foundation; either version 3 of the License, or 9%% (at your option) any later version. 10%% 11%% SECS1D is distributed in the hope that it will be useful, 12%% but WITHOUT ANY WARRANTY; without even the implied warranty of 13%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14%% GNU General Public License for more details. 15%% 16%% You should have received a copy of the GNU General Public License 17%% along with SECS1D; If not, see <http://www.gnu.org/licenses/>. 18 19%% 20%% Solve the scaled stationary bipolar DD equation system using Gummel algorithm. 21%% 22%% [n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_gummel_map (x, D, Na, Nd, 23%% pin, nin, Vin, Fnin, 24%% Fpin, l2, er, u0n, 25%% uminn, vsatn, betan, 26%% Nrefn, u0p, uminp, vsatp, 27%% betap, Nrefp, theta, tn, tp, 28%% Cn, Cp, an, ap, Ecritnin, Ecritpin, 29%% toll, maxit, ptoll, pmaxit) 30%% 31%% input: 32%% x spatial grid 33%% D, Na, Nd doping profile 34%% pin initial guess for hole concentration 35%% nin initial guess for electron concentration 36%% Vin initial guess for electrostatic potential 37%% Fnin initial guess for electron Fermi potential 38%% Fpin initial guess for hole Fermi potential 39%% l2 scaled Debye length squared 40%% er relative electric permittivity 41%% u0n, uminn, vsatn, Nrefn electron mobility model coefficients 42%% u0p, uminp, vsatp, Nrefp hole mobility model coefficients 43%% theta intrinsic carrier density 44%% tn, tp, Cn, Cp, 45%% an, ap, 46%% Ecritnin, Ecritpin generation recombination model parameters 47%% toll tolerance for Gummel iterarion convergence test 48%% maxit maximum number of Gummel iterarions 49%% ptoll convergence test tolerance for the non linear 50%% Poisson solver 51%% pmaxit maximum number of Newton iterarions 52%% 53%% output: 54%% n electron concentration 55%% p hole concentration 56%% V electrostatic potential 57%% Fn electron Fermi potential 58%% Fp hole Fermi potential 59%% Jn electron current density 60%% Jp hole current density 61%% it number of Gummel iterations performed 62%% res total potential increment at each step 63 64function [n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_gummel_map (x, D, Na, Nd, pin, nin, Vin, 65 Fnin, Fpin, l2, er, u0n, uminn, 66 vsatn,betan,Nrefn, u0p, 67 uminp,vsatp,betap,Nrefp, 68 theta, tn, tp, Cn, Cp, an, ap, 69 Ecritnin, Ecritpin, toll, maxit, 70 ptoll, pmaxit) 71 72 p = pin; 73 n = nin; 74 V = Vin; 75 Fp = Fpin; 76 Fn = Fnin; 77 dx = diff (x); 78 dxm = (dx(1:end-1) + dx(2:end)); 79 80 Nnodes = numel (x); 81 Nelements = Nnodes -1; 82 83 Jn = zeros (Nelements, 1); 84 Jp = zeros (Nelements, 1); 85 86 for it = 1:maxit 87 88 [V(:,2), n(:,2), p(:,2)] = secs1d_nlpoisson_newton (x, [1:Nnodes], V(:,1), n(:, 1), 89 p(:,1), Fn(:,1), Fp(:,1), D, l2, 90 er, ptoll, pmaxit); 91 92 dV = diff (V(:, 2)); 93 E = - dV ./ dx; 94 [Bp, Bm] = bimu_bernoulli (dV); 95 96 [Rn, Rp, Gn, Gp, II] = generation_recombination_model (x, n(:, end), p(:, end), 97 E, Jn, Jp, tn, tp, theta, 98 Cn, Cp, an, ap, Ecritnin, Ecritpin); 99 100 mobility = mobility_model (x, Na, Nd, Nrefn, E, u0n, uminn, vsatn, betan); 101 102 A = bim1a_advection_diffusion (x, mobility, 1, 1, V(:, 2)); 103 M = bim1a_reaction (x, 1, Rn) + bim1a_reaction (x, II, 1); 104 105 R = bim1a_rhs (x, 1, Gn); 106 107 A = A + M; 108 109 n(:,3) = nin; 110 n(2:end-1,3) = A(2:end-1, 2:end-1) \ (R(2:end-1) - A(2:end-1, [1 end]) * nin ([1 end])); 111 Fn(:,2) = V(:,2) - log (n(:, 3)); 112 Jn = mobility .* (n(2:end, 2) .* Bp - n(1:end-1, 2) .* Bm) ./ dx; 113 114 [Rn, Rp, Gn, Gp, II] = generation_recombination_model (x, n(:, end), p(:, end), 115 E, Jn, Jp, tn, tp, theta, 116 Cn, Cp, an, ap, Ecritnin, Ecritpin); 117 118 mobility = mobility_model (x, Na, Nd, Nrefp, E, u0p, uminp, vsatp, betap); 119 120 A = bim1a_advection_diffusion (x, mobility, 1, 1, -V(:, 2)); 121 M = bim1a_reaction (x, 1, Rp) + bim1a_reaction (x, II, 1); 122 R = bim1a_rhs (x, 1, Gp); 123 A = A + M; 124 125 p(:,3) = pin; 126 p(2:end-1,3) = A(2:end-1, 2:end-1) \ (R(2:end-1) - A(2:end-1, [1 end]) * pin ([1 end])); 127 Fp(:,2) = V(:,2) + log (p(:,3)); 128 Jp = -mobility .* (p(2:end, 2) .* Bm - p(1:end-1, 2) .* Bp) ./ dx; 129 130 nrfn = norm (Fn(:,2) - Fn(:,1), inf); 131 nrfp = norm (Fp(:,2) - Fp(:,1), inf); 132 nrv = norm (V(:,2) - V(:,1), inf); 133 res(it) = max ([nrfn; nrfp; nrv]); 134 135 if (res(it) < toll) 136 break 137 endif 138 139 V(:,1) = V(:,end); 140 p(:,1) = p(:,end) ; 141 n(:,1) = n(:,end); 142 Fn(:,1) = Fn(:,end); 143 Fp(:,1) = Fp(:,end); 144 145 endfor 146 147 n = n(:,end); 148 p = p(:,end); 149 V = V(:,end); 150 Fn = Fn(:,end); 151 Fp = Fp(:,end); 152 153endfunction 154 155function u = mobility_model (x, Na, Nd, Nref, E, u0, umin, vsat, beta) 156 157 Neff = Na + Nd; 158 Neff = (Neff(1:end-1) + Neff(2:end)) / 2; 159 160 ubar = umin + (u0 - umin) ./ (1 + (Neff ./ Nref) .^ beta); 161 u = 2 * ubar ./ (1 + sqrt (1 + 4 * (ubar .* abs (E) ./ vsat) .^ 2)); 162 163endfunction 164 165function [Rn, Rp, Gn, Gp, II] = generation_recombination_model (x, n, p, E, Jn, Jp, tn, tp, 166 theta, Cn, Cp, an, ap, Ecritn, 167 Ecritp) 168 169 denomsrh = tn .* (p + theta) + tp .* (n + theta); 170 factauger = Cn .* n + Cp .* p; 171 fact = (1 ./ denomsrh + factauger); 172 173 Rn = p .* fact; 174 Rp = n .* fact; 175 176 Gn = theta .^ 2 .* fact; 177 Gp = Gn; 178 179 II = an * exp(-Ecritn./abs(E)) .* abs (Jn) + ap * exp(-Ecritp./abs(E)) .* abs (Jp); 180 181endfunction 182 183 184%!demo 185%! % physical constants and parameters 186%! secs1d_physical_constants; 187%! secs1d_silicon_material_properties; 188%! 189%! % geometry 190%! L = 10e-6; % [m] 191%! xm = L/2; 192%! 193%! Nelements = 1000; 194%! x = linspace (0, L, Nelements+1)'; 195%! sinodes = [1:length(x)]; 196%! 197%! % dielectric constant (silicon) 198%! er = esir * ones (Nelements, 1); 199%! 200%! % doping profile [m^{-3}] 201%! Na = 1e23 * (x <= xm); 202%! Nd = 1e23 * (x > xm); 203%! 204%! % avoid zero doping 205%! D = Nd - Na; 206%! 207%! % initial guess for n, p, V, phin, phip 208%! V_p = -1; 209%! V_n = 0; 210%! 211%! Fp = V_p * (x <= xm); 212%! Fn = Fp; 213%! 214%! p = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni./abs(D)) .^2)) .* (x <= xm) + ... 215%! ni^2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x > xm); 216%! 217%! n = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^ 2)) .* (x > xm) + ... 218%! ni ^ 2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x <= xm); 219%! 220%! V = Fn + Vth * log (n / ni); 221%! 222%! % scaling factors 223%! xbar = L; % [m] 224%! nbar = norm(D, 'inf'); % [m^{-3}] 225%! Vbar = Vth; % [V] 226%! mubar = max (u0n, u0p); % [m^2 V^{-1} s^{-1}] 227%! tbar = xbar^2 / (mubar * Vbar); % [s] 228%! Rbar = nbar / tbar; % [m^{-3} s^{-1}] 229%! Ebar = Vbar / xbar; % [V m^{-1}] 230%! Jbar = q * mubar * nbar * Ebar; % [A m^{-2}] 231%! CAubar = Rbar / nbar^3; % [m^6 s^{-1}] 232%! abar = 1/xbar; % [m^{-1}] 233%! 234%! % scaling procedure 235%! l2 = e0 * Vbar / (q * nbar * xbar^2); 236%! theta = ni / nbar; 237%! 238%! xin = x / xbar; 239%! Din = D / nbar; 240%! Nain = Na / nbar; 241%! Ndin = Nd / nbar; 242%! pin = p / nbar; 243%! nin = n / nbar; 244%! Vin = V / Vbar; 245%! Fnin = Vin - log (nin); 246%! Fpin = Vin + log (pin); 247%! 248%! tnin = tn / tbar; 249%! tpin = tp / tbar; 250%! 251%! u0nin = u0n / mubar; 252%! uminnin = uminn / mubar; 253%! vsatnin = vsatn / (mubar * Ebar); 254%! 255%! u0pin = u0p / mubar; 256%! uminpin = uminp / mubar; 257%! vsatpin = vsatp / (mubar * Ebar); 258%! 259%! Nrefnin = Nrefn / nbar; 260%! Nrefpin = Nrefp / nbar; 261%! 262%! Cnin = Cn / CAubar; 263%! Cpin = Cp / CAubar; 264%! 265%! anin = an / abar; 266%! apin = ap / abar; 267%! Ecritnin = Ecritn / Ebar; 268%! Ecritpin = Ecritp / Ebar; 269%! 270%! % tolerances for convergence checks 271%! toll = 1e-3; 272%! maxit = 1000; 273%! ptoll = 1e-12; 274%! pmaxit = 1000; 275%! 276%! % solve the problem using the full DD model 277%! [nout, pout, Vout, Fnout, Fpout, Jnout, Jpout, it, res] = ... 278%! secs1d_dd_gummel_map (xin, Din, Nain, Ndin, pin, nin, Vin, Fnin, Fpin, ... 279%! l2, er, u0nin, uminnin, vsatnin, betan, Nrefnin, ... 280%! u0pin, uminpin, vsatpin, betap, Nrefpin, theta, ... 281%! tnin, tpin, Cnin, Cpin, anin, apin, ... 282%! Ecritnin, Ecritpin, toll, maxit, ptoll, pmaxit); 283%! 284%! % Descaling procedure 285%! n = nout*nbar; 286%! p = pout*nbar; 287%! V = Vout*Vbar; 288%! Fn = V - Vth*log(n/ni); 289%! Fp = V + Vth*log(p/ni); 290%! dV = diff(V); 291%! dx = diff(x); 292%! E = -dV./dx; 293%! 294%! % band structure 295%! Efn = -Fn; 296%! Efp = -Fp; 297%! Ec = Vth*log(Nc./n)+Efn; 298%! Ev = -Vth*log(Nv./p)+Efp; 299%! 300%! plot (x, Efn, x, Efp, x, Ec, x, Ev) 301%! legend ('Efn', 'Efp', 'Ec', 'Ev') 302%! axis tight