PN Diode ProjectPN结二极管项目.doc

ACUTE: PN Diode Modeling Dragica Vasileska and Gerhard Klimeck EQUILIBRIUM SOLVER: You are provided with a MATLAB script of an equilibrium 1D Poisson equation solver for a pn-diode. Please try to understand and run the code for the following doping densities: (a) NA = 1016 cm-3, ND =1016 cm-3 (b) NA = 1016 cm-3, ND =1018 cm-3 (c) NA = 1018 cm-3, ND =1018 cm-3 For each of these cases plot the potential and electric field profiles, the electron and hole densities and the total charge densities. From these plots extract the width of the depletion region and the peak electric field. Compare your simulated data with the depletion charge approximation analytical results. NON-EQUILIBRIUM SOLVER: Develop a one-dimensional (1D) drift-diffusion simulator for modeling pn-junctions (diodes) under forward and reverse bias conditions. Include both types of carriers in your model (electrons and holes). Use the finite-difference expressions for the electron and hole current continuity equations that utilize the Sharfetter-Gummel discretization scheme. Model: Silicon diode, with permittivity and intrinsic carrier concentration F/m 1005 . 1 10 sc at T=300K. In all your simulations assume that T=300K. Use concentration- 310 cm 105 . 1 i n dependent and field-dependent mobility models and SRH generation-recombination process. Assume ohmic contacts and charge neutrality at both ends to get the appropriate boundary conditions for the potential and the electron and hole concentrations. For the electron and hole mobility use 1500 and 1000 cm2/V-s, respectively. For the SRH generation-recombination, use TAUN0=TAUP0=0.1 us. To simplify your calculations, assume that the trap energy level coincides with the intrinsic level. Doping: Use and as a net doping of the p- and n-regions, Ncm A 1016 3 Ncm D 1017 3 respectively. Numerical methods: Use the LU decomposition method for the solution of the 1D Poisson and the two 1D continuity equations for electrons and holes individually. Use Gummel's decoupled scheme, described in the class, to solve the resultant set of coupled set of algebraic equations. Outputs: Plot the conduction band edge under equilibrium conditions (no current flow) and for VA=0.625 V. Plot the electron and hole densities under equilibrium conditions (no current flow) and for VA=0.625 V. Plot the electric field profile under equilibrium conditions (no current flow) and for VA=0.625 V. Vary the Anode bias from 0 to 0.625 V, in voltage increments that are fraction of the VA thermal voltage , to have stable convergence. Plot the resulting I-V Vk T q TB / characteristics. The current will be in A/unit area, since you are doing 1D modeling. Check the conservation of current when going from the cathode to the anode, which also means conservation of particles in your system. For the calculation of the current density, use the results given in the notes. For =0.625 V, plot the position of the electron and hole quasi-Fermi levels, with respect VA to the equilibrium Fermi level, assumed to be the reference energy level. Final note: When you submit your project report, in addition to the final results, give a brief explanation of the problem you are solving with reference to the listing of your program that you need to turn in with the report. % % % % 1D Poisson Equation Solver for pn Diodes % % % % % Defining the Fundamental and Material Constants % q = 1.602E-19; % C or J/eV kb = 1.38E-23; % J/K eps = 1.05E-12; % This includes the eps = 11.7 for Si F/cm T = 300; % K ni = 1.45E10; % Intrinsic carrier concentration 1/cm3 Vt = kb*T/q; % eV RNc = 2.82E19; % Effective DOS of the conduction Band dEc = Vt*log(RNc/ni); % Define Doping Values % Na = 1E18; % 1/cm3 Nd = 1E18; % 1/cm3 % Calculate relevant parameters for the simulation % Vbi = Vt*log(Na*Nd/(ni*ni); W = sqrt(2*eps*(Na+Nd)*Vbi/(q*Na*Nd) % cm Wn = W*sqrt(Na/(Na+Nd) % cm Wp = W*sqrt(Nd/(Na+Nd) % cm Wone = sqrt(2*eps*Vbi/(q*Na) % cm E_p = q*Nd*Wn/eps % V/cm Ldn = sqrt(eps*Vt/(q*Nd); Ldp = sqrt(eps*Vt/(q*Na); Ldi = sqrt(eps*Vt/(q*ni); % Calculate relevant parameters in an input file % % Write to a file save input_params.txt Na Nd Vbi W Wn Wp E_p Ldn Ldp %Material_Constants %Define some material constants % Setting the size of the simulation domain based % on the analytical results for the width of the depletion regions % for a simple pn-diode % x_max = 0; if(x_max Ldp) dx=Ldp; end dx = dx/20; % Calculate the required number of grid points and renormalize dx % n_max = x_max/dx; n_max = round(n_max); dx = dx/Ldi; % Renormalize lengths with Ldi % Set up the doping C(x) = Nd(x) - Na(x) that is normalized with ni % for i = 1:n_max if(i n_max/2) dop(i) = Nd/ni; end end % Initialize the potential based on the requirement of charge % neutrality throughout the whole structure for i = 1: n_max zz = 0.5*dop(i); if(zz 0) xx = zz*(1 + sqrt(1+1/(zz*zz); elseif(zz delta_max) delta_max=xx; end %sprintf('delta_max = %d',delta_max) %'k_iter = %d',k_iter,' end %delta_max=max(abs(delta); % Test convergence and recalculate forcing function and % central coefficient b if necessary if(delta_max delta_acc) flag_conv = 1; else for i = 2: n_max-1 b(i) = -(2/dx2 + exp(fi(i) + exp(-fi(i); f(i) = exp(fi(i) - exp(-fi(i) - dop(i) - fi(i)*(exp(fi(i) + exp(-fi(i); end end end % Write the results of the simulation in files % xx1(1) = dx*1e4; for i = 2:n_max-1 Ec(i) = dEc - Vt*fi(i); %Values from the second Node% ro(i) = -ni*(exp(fi(i) - exp(-fi(i) - dop(i); el_field1(i) = -(fi(i+1) - fi(i)*Vt/(dx*Ldi); el_field2(i) = -(fi(i+1) - fi(i-1)*Vt/(2*dx*Ldi); n(i) = exp(fi(i); p(i) = exp(-fi(i); xx1(i) = xx1(i-1) + dx*Ldi*1e4; end Ec(1) = Ec(2); Ec(n_max) = Ec(n_max-1); xx1(n_max) = xx1(n_max-1) + dx*Ldi*1e4; el_field1(1) = el_field1(2); el_field2(1) = el_field2(2); el_field1(n_max) = el_field1(n_max-1); el_field2(n_max) = el_field2(n_max-1); nf = n*ni; pf = p*ni; ro(1) = ro(2); ro(n_max) = ro(n_max-1); figure(1) plot(xx1, fi) xlabel('x um'); ylabel('Potential eV'); title('Potential vs Position'); figure(2) plot(xx1, el_field1) hold on; plot(xx1, el_field2) xlabel('x um'); ylabel('Electric Field V/cm'); title('Field Profile vs Position'); figure(3) plot(xx1, nf) hold on; plot(xx1, pf) xlabel('x um'); ylabel('Electron title('Electron legend('n','p'); %axis(0 6.75 0 10.2e17) figure(4) plot(xx1, ro) xlabel('x um'); ylabel('Total Charge Density 1/cm3'); title('Total Charge Density vs Position'); %axis(0.5 5 -3e17 8e17) %figure(5) %plot(xx1, n) %hold all %plot(xx1, p) save cond_band xx1 Ec; save tot_charge xx1 ro; save el_field xx1 el_field1 el_field2; save np_data xx1 nf pf; save pot_1 fi; clear all;