Расчёт конструкции интегрального конденсатора.
Цель работы: написание программы на языке MATLAB, способной рассчитать площадь и толщину диэлектрика пленочного конденсатора.
Основные теоретические положения:
Емкость конденсатора будет определяться посредством соотношения:
Относительная погрешность емкости вычисляется как:
Упростим задачу, рассмотрением квадратного конденсатора. В этом случае:
Относительная погрешность:
Вычислить относительную погрешность температурной составляющей диэлектрической проницаемости можно зная ТКE материала:
Выразим относительную погрешность толщины диэлектрика:
В конечном итоге имеем:
Тогда:
В конечном итоге выберем d как:
Идеальных конденсаторов не существует – у каждого конденсатора есть резистивные потери, которые можно учесть, введя в рассмотрение эквивалентную схему:
Рис. 2. – Эквивалентная схема реального конденсатора.
Одним из важнейших параметров конденсатора является его добротность, которая определяется как:
В этой формуле пассивная мощность:
Активная мощность:
Тогда получим формулу:
После преобразований:
Таким
образом введя параметр
можно характеризовать добротность
конденсатора.
Текст программы и результаты расчетов:
clear all
close all
clc
Up=1; %напряжение питания
Kz=2; %коэффициент запаса
Epr=50e6; %пробивная напряженность
delta_d=0.9e-9; %разброс толщины (задается технологами)
delta_a=9e-8; %разброс длины
delta_C=0.03; %необходимый допуск
C=5e-12; %необходимое значение емкости
delta_eps_c=3e-4; %разброс параметров диэлектрической проницаемости (задается технологами)
TKE=2e-5; %температурный коэффициент (двуокись кремния)
d_min=90e-9; %чтобы не "свистел” ток
Lambda=5e-6; %проектная норма
eps_0=8.85e-12;
eps=4;
Tmin=-40;
Tnom=27;
Tmax=85;
Tmin=Tmin+273;
Tmax=Tmax+273;
Tnom=Tnom+273;
d_u=Up*Kz/Epr;
a=sqrt((d_u*C)/(eps*eps_0));
if a/Lambda-fix(a/Lambda)>0 % остаток от деления
ares=(fix(a/Lambda)+1)*Lambda; %fix - отбрасывает дробную часть, округляем в большую сторону
else
ares=a;
end
ares=ares*1e6 % значение в мкм
d_delta=sqrt(delta_d^2/(delta_C^2-2*(delta_a/a)^2-delta_eps_c^2-(TKE*max([Tmax-Tnom Tnom-Tmin]))^2));
d=max([d_min d_u d_delta]);
d=d*1e9 % значение в нм
f=1e3:500:1e8;
R=1;
l=1e-9;
r=1e9;
ts=C*r;
%построение графика зависимости добротности от частоты
Q=(f.*(2*pi*ts))./(1+(f.^2).*(4*pi^2*2*C*R*ts));
loglog(f,Q, 'LineWidth',2)
xlabel ('f, Hz', 'FontSize',19)
ylabel ('Q', 'FontSize',19)
grid on
Рис. 3. – результат расчетов.
(сторона
обкладки)
(толщина
диэлектрика)
Рис. 4. - зависимость добротности от частоты.
Численное моделирование стационарных процессов в диоде.
Цель работы: написание программы на языке MATLAB, способной рассчитать ВольтАмперную характеристику полупроводникового диода с заданными технологическими параметрами (концентрация, ширина и пр.).
Основные теоретические положения:
В качестве диода можно использовать любой из двух р-n-переходов, расположенных в изолирующем кармане: эмиттерный или коллекторный. Можно также использовать их комбинации. Поэтому, по существу, интегральный диод представляет собой диодное включение интегрального транзистора.
Пять возможных вариантов диодного включения транзистора показаны на рис.1. В табл. 1 приведены типичные параметры этих вариантов. Для них приняты следующие обозначения: до черточки стоит обозначение анода, после черточки - катода; если два слоя соединены, их обозначения пишутся слитно. Из табл. 1 видно, что варианты различаются как по статическим, так и по динамическим параметрам.
Рис. 5. - Интегральные диоды (диодные включения транзистора).
Пробивные напряжения Uпр зависят от используемого перехода: они меньше у тех вариантов, в которых используется эмиттерный переход.
Обратные токи Iобр (без учета токов утечки) — это токи термогенерации в переходах. Они зависят от объема перехода и, следовательно, меньше у тех вариантов, у которых используется только эмиттерный переход, имеющий наименьшую площадь.
Таблица 7 - Типичные параметры интегральных диодов.
Время восстановления обратного тока tв (т.е. время переключения диода из открытого в закрытое состояние) минимально у варианта БК-Э; у этого варианта заряд накапливается только в базовом слое (так как коллекторный переход закорочен). У других вариантов заряд накапливается не только в базе J но и в коллекторе, так что для рассасывания заряда требуется большее время.
Фундаментальная система уравнений:
Процессы переноса и накопления зарядов в полупроводниковых диодах можно описать фундаментальной системой уравнений в диффузионно-дрейфовом приближении (подвижность непосредственно зависит от напряженности поля), в которую входят уравнения непрерывности для электронов и дырок (1), (2), уравнения плотностей электронной и дырочной составляющих электрического тока (3), (4), и уравнение Пуассона (5).
где
-
коэффициенты диффузии электронов и
дырок,
– подвижности электронов и дырок.
В уравнениях плотности тока (3) и (4) первые слагаемые выражают дрейфовую составляющую, определяемую напряженностью электрического поля, а вторые – диффузионную составляющую, определяемую градиентами концентраций электронов и дырок.
Учитывая
и соотношения Эйнштейна
(
),
уравнения (3) и (4) можно представить:
Подставляя (6), (7) в уравнения непрерывности (1), (2), выражая объемную плотность зарядов через концентрации подвижных носителей и ионизированных примесей и считая диэлектрическую проницаемость константой, фундаментальная система уравнений в векторной форме будет выглядеть следующим образом:
где
– концентрации ионов доноров и акцепторов.
В операторной форме:
Таблица 8 – нормировочные коэффициенты.
Нормируемая физическая величина |
Нормировочный коэффициент |
Время |
|
Потенциал |
|
Концентрация |
|
Длина |
|
Коэффициент диффузии |
|
Подвижность |
|
– собственная
концентрация электронов проводимости
или дырок в собственном полупроводнике
Разделив все переменные и параметры на нормировочные коэффициенты, фундаментальная система уравнений в частных производных в нормированном виде примет следующий вид:
Искомыми функциями при решении системы (14) – (16) являются зависимости от координат и времени потенциала и концентраций электронов и дырок для определенных моделей подвижностей и генерации-рекомбинации при известных распределениях по координатам концентраций доноров и акцепторов.
Множество искомых функций системы уравнений – базис переменных. Можно приводить систему уравнений к другим базисам. В данном случае переход к другому базису обусловлен увеличением числа итераций в процессе численного решения.
Распространенными базисами являются:
-
– концентрации электронов и дырок,
электрический потенциал;
-
– квазиуровни Ферми для электронов и
дырок, электрический потенциал;
-
– экспоненты квазиуровней Ферми для
электронов и дырок, электрический
потенциал.
Переход к другому базису переменных осуществляется с использованием выражений, связывающих переменные различных базисов:
Для поиска единственного решения данной задачи необходимо задать граничные и начальные условия, используя следующие упрощения: материал однородный, не учитываются краевые эффекты и изменения концентрации по координате в латеральных направлениях (ось y, z), процесс стационарный.
Используя уравнение электронейтральности и закон действующих масс, условие на границе с p-областью будет представлено в следующем виде:
Отсюда,
Аналогично для p-области:
Выражение
(24) описывает граничное условие в
начальной точке координат n-области,
где
,
выражения (25) и (27) соответствуют граничным
условиям на границе n
и p-областей,
(26) - граничное условие в конечной точке
координат p-области,
где
.
Если
принять потенциал середины запрещенной
зоны равным отрицательному
электростатическому потенциалу
,
то граничные условия для потенциала:
где
– потенциал квазиуровня Ферми для
электронов и дырок,
– энергетический потенциал,
– тепловой потенциал
В терминах базиса Слотбума:
где
,
– экспонента квазиуровней Ферми для
электронов и дырок
Граничные условия:
Тогда,
,
где U – внешнее напряжение.
Следовательно,
Таким образом,
Поиск решения системы уравнений можно осуществить с помощью метода конечных разностей, который предполагает дискретизацию дифференциальных уравнений на прямоугольных координатных сетках, то есть на сетках, элементарные ячейки которых представляют собой прямоугольники для двух измерений или параллелепипеды для трех измерений.
Таким образом, уравнение Пуассона по методу конечных разностей имеет вид:
Уравнения непрерывности:
Граничные условия:
Текст программы и результаты метода Ньютона:
clear all
close all
clc
warning off
constants; % инициализация необходимых физических констант
parameters; % инициализация шага на сетке, значения относительной погрешности к которой стремимся,
structure; % инициализация структуры, длины p и n областей, концентрации доноров и акцепторов в данных областях,
% шага для подгонки, напряжения приложенного к диоду, расчет концентраций
% собственых носителей заряда, а также параметров необходимых для расчета
% подвижности электронов и дырок
grid_functions; % задание сетки координат, концентраций и подвижностей с разным шагом для каждой области.
% загон конечной длины получившейся сетки в переменную I, а также
% инициализация матрицы разностей двух координат (x(2) - x(1)) и т.д.
normirovka; % нормировка координат, матрицы разностей, напряжения, концентрации, подвижностей
prog_param; % значения для толщин линий и размера шрифта, название зашитое в переменной
% для будущих файлов
fi0(1)=log(N(1)/2+sqrt((N(1)/2)^2+1));% граничные условия для потенциала. В начале
fi0(I)=-log(-N(I)/2+sqrt((N(I)/2)^2+1));% В конце
for i=2:I-1 % проход по всей сетке кроме граничных случаев, расчет электростатических
% потенциалов (далее электропотенциал) в ней
fi0(i)=fi0(1)+(fi0(I)-fi0(1))/(x(I)-x(1))*(x(i)-x(1));
end
fi0=fi0'; % транспонирование
E=zeros(I,1); % инициализация матрицы единиц для напряженности электрического
% поля с I строк (кол-во точек в сетке) и 1 столбцом
for j=1:length(U) % проход по всем заданным напряжениям
if j>1 % условие для задания начальных значений, если напряжение является не 1 по списку
for i=1:I
fi0(i)=fi0(i)+(fi0(i)-fi0(1))/(fi0(I)-fi0(1))*(U(j)-U(j-1)); % перерасчет электропотенциалов
end
Fn0=ones(I,1); % сброс прошлой матрицы для экспоненты квазиуровня ферми для электронов + задание начального значения
Fn0(I)=exp(-U(j)); % задание конечного значения для экспоненты квазиуровня ферми для электронов
Fp0=ones(I,1);
Fp0(I)=exp(U(j)); % повтор описанного выше, но уже для дырок
f0=[fi0; Fn0; Fp0]; % запись всех требующихся для нас параметров в одну матрицу
end
ct=0; % инициализация счетчика
err=100*Delta; % задание ошибки такой, чтоб начал работать цикл
while err>Delta
ct=ct+1; % увеличение счетчика на единицу
if U(j)==0 % проверка для определения 0 ли напряжение по списку, и если 0 то
Fn0=ones(I,1); % инициализация матрицы для экспоненты квазиуровня ферми для электронов + задание начального значения
Fp0=ones(I,1); % тоже самое но для дырок
A=zeros(I,I); % инициализация матрицы Якоби
B=zeros(I,1); % инициализация матрицы значений потенциала
A(1,1)=1; % задание начального значения матрицы Якоби
B(1)=fi0(1)-log(N(1)/2+sqrt((N(1)/2)^2+1)); % задание начального значения матрицы значений потенциала
A(I,I)=1; % задание конечного значения матрицы Якоби
B(I)=fi0(I)+log(-N(I)/2+sqrt((N(I)/2)^2+1))-U(j); % задание конечного значения матрицы значений потенциала
for i=2:I-1 % проход по всем элементам матриц с расчетом каждой из них
A(i,i)=2/(dx(i)+dx(i-1))*(-1/dx(i)-1/dx(i-1))-Fn0(i)*exp(fi0(i))-Fp0(i)*exp(-fi0(i));
A(i,i+1)=2/(dx(i)+dx(i-1))/dx(i); % расчет частной производной по fi(i+1) - fi(i)
A(i,i-1)=2/(dx(i)+dx(i-1))/dx(i-1); % расчет частной производной по fi(i) - fi(i-1)
B(i)=2/(dx(i)+dx(i-1))*((fi0(i+1)-fi0(i))/dx(i)-(fi0(i)-fi0(i-1))/dx(i-1))-...
Fn0(i)*exp(fi0(i))+Fp0(i)*exp(-fi0(i))+N(i); % расчет частной производной по потенциалу для формулы градиента потенциала
end
fi=fi0-A^(-1)*B; % расчет разницы между начальными условиями и рассчитанными
% (-1 в степени А - обратная матрица)
clc % стирание из командного окна всей информации
Voltage=U(j)*V0 % денормирование напряжения, вывод на экран
err=max(abs(fi-fi0))/max(abs(fi0)) % заменяем невязку на абсолютную погрешность
ct % вывод счетчика итераций
fi0=fi0+(fi-fi0)./K(j); % расчет новых начальных условий
else % если напряжение не 1 по списку - считаем по этой части
A=zeros(3*I,3*I); % инициализация матрицы Якоби сразу для
% электропотенциала и экспонент квазиуровня ферми
B=zeros(3*I,1); % инициализация матрицы значений функций (fi0, Fn0 и Fp0)
% задание начальных значений для
A(1,1)=1; % матрицы Якоби, отвечающей за fi0
B(1)=fi0(1)-log(N(1)/2+sqrt((N(1)/2)^2+1)); % матрицы значений для fi0
% задание конечных значений для:
A(I,I)=1; % матрицы Якоби, отвечающей за fi0
B(I)=fi0(I)+log(-N(I)/2+sqrt((N(I)/2)^2+1))-U(j); % матрицы значений для fi0
% задание начальных значений для:
for i=2:I-1 % расчет значений для матрицы Якоби
% и расчет значений для матрицы значений функции
% частные производные для формулы градиента потенциала
A(i,i)=2/(dx(i)+dx(i-1))*(-1/dx(i)-1/dx(i-1))-Fn0(i)*exp(fi0(i))-Fp0(i)*exp(-fi0(i)); % расчет частной производной по fi(i)
A(i,i+1)=2/(dx(i)+dx(i-1))/dx(i); % расчет частной производной по fi(i+1) - fi(i)
A(i,i-1)=2/(dx(i)+dx(i-1))/dx(i-1); % расчет частной производной по fi(i) - fi(i-1)
A(i,I+i)=-exp(fi0(i)); % расчет частной производной по Fn0
A(i,2*I+i)=exp(-fi0(i)); % расчет частной производной по Fp0
B(i)=2/(dx(i)+dx(i-1))*((fi0(i+1)-fi0(i))/dx(i)-(fi0(i)-fi0(i-1))/dx(i-1))-...
Fn0(i)*exp(fi0(i))+Fp0(i)*exp(-fi0(i))+N(i); % расчет значения функции
% частные производные для формулы градиента экспоненты
% квазиуровня Ферми для электронов
end
A(I+1,I+1)=1;
B(I+1)=Fn0(1)-1;
A(2*I,2*I)=1;
B(2*I)=Fn0(I)-exp(-U(j));
for i=2:I-1
A(I+i,I+i)=-(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i)-...
(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fn(i)
A(I+i,I+i+1)=(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i); % расчет частной производной по Fn(i+1)-Fn(i)
A(I+i,I+i-1)=(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fn(i)-Fn(i-1)
A(I+i,i)=mun(i)/2*exp(fi0(i))*(Fn0(i+1)-Fn0(i))/dx(i)-...
mun(i)/2*exp(fi0(i))*(Fn0(i)-Fn0(i-1))/dx(i-1); % расчет частной производной по fi(i)
A(I+i,i+1)=mun(i+1)/2*exp(fi0(i+1))*(Fn0(i+1)-Fn0(i))/dx(i); % расчет частной производной по fi(i+1)
A(I+i,i-1)=-mun(i-1)/2*exp(fi0(i-1))*(Fn0(i)-Fn0(i-1))/dx(i-1); % расчет частной производной по fi(i-1)
B(I+i)=(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2*(Fn0(i+1)-Fn0(i))/dx(i)-...
(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2*(Fn0(i)-Fn0(i-1))/dx(i-1); % расчет значения функции
% частные производные для формулы градиента экспоненты
% квазиуровня Ферми для электронов
end
A(2*I+1,2*I+1)=1;
B(2*I+1)=Fp0(1)-1;
A(3*I,3*I)=1;
B(3*I)=Fp0(I)-exp(U(j));
for i=2:I-1
A(2*I+i,2*I+i)=-(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2/dx(i)-...
(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fp(i)
A(2*I+i,2*I+i+1)=(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2/dx(i); % расчет частной производной по Fp(i+1)-Fp(i)
A(2*I+i,2*I+i-1)=(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fp(i)-Fp(i-1)
A(2*I+i,i)=-mup(i)/2*exp(-fi0(i))*(Fp0(i+1)-Fp0(i))/dx(i)+...
mup(i)/2*exp(-fi0(i))*(Fp0(i)-Fp0(i-1))/dx(i-1); % расчет частной производной по fi(i)
A(2*I+i,i+1)=-mup(i+1)/2*exp(-fi0(i+1))*(Fp0(i+1)-Fp0(i))/dx(i); % расчет частной производной по fi(i+1)
A(2*I+i,i-1)=mup(i-1)/2*exp(-fi0(i-1))*(Fp0(i)-Fp0(i-1))/dx(i-1); % расчет частной производной по fi(i-1)
B(2*I+i)=(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2*(Fp0(i+1)-Fp0(i))/dx(i)-...
(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2*(Fp0(i)-Fp0(i-1))/dx(i-1); % расчет значения функции
end
f=f0-A^(-1)*B; % расчет разницы между начальными условиями и рассчитанными
clc % очистка командного окна
Voltage=U(j)*V0 % денормирование напряжения, вывод на экран
err=max(abs(f-f0))/max(abs(f0)) % заменяем невязку на абсолютную погрешность
ct % вывод счетчика итераций
f0=f0+(f-f0)./K(j); % расчет новых начальных условий
fi0=f0(1:I);
Fn0=f0(I+1:2*I);
Fp0=f0(2*I+1:3*I); % замена предыдущих значений новыми
end
n=Fn0.*exp(fi0); % расчет концентрации электронов
p=Fp0.*exp(-fi0); % расчет концентрации дырок
for i=1:I-1
E(i)=-(fi0(i+1)-fi0(i))/dx(i); % расчет напряженности электрического поля
end
E(I)=E(I-1)+dx(I-1)/dx(I-2)*(E(I-1)-E(I-2)); % так как не можем посчитать последнее значение для
% поля ввиду отсутсвия следующего значения потенциала - принимаем
% его за последний рассчитанный
Ro=p-n+N;
end
for i=2:I-1 % так как плотность тока меняться внутри диода не будет -
% принимаем всю плотность тока за рассчитанную
Jn(i)=mun(i)*exp(fi0(i))*(Fn0(i+1)-Fn0(i-1))/(dx(i)+dx(i-1)); % расчет плотности тока электронов
Jp(i)=-mup(i)*exp(-fi0(i))*(Fp0(i+1)-Fp0(i-1))/(dx(i)+dx(i-1)); % расчет плотности ток дырок
end
Jn(I)=Jn(I-1);
Jp(I)=Jp(I-1);
Jn(1)=Jn(2);
Jp(1)=Jp(2);
J=Jn+Jp; % суммарная плотность тока
Current(j)=J(I); % запоминаем плотность тока для построения ВАХ
graphics;% вывод графиков
end
figure % создание отдельного окна
plot(U.*V0,-Current.*(J0*1e-4),'LineWidth',LW) % построение ВАХ
xlabel('Voltage, V','FontSize',FS) % именуем ось абсцисс
ylabel('Current density, A/cm^2','FontSize',FS) % именуем ось ординат
grid on % включаем сетку на графике
print(gcf,'-djpeg','VAH') % сохраняем график ВАХ в виде файла джипег
Результаты расчетов:
Voltage = 0.8000
err = 2.8880e-08
Полученные графики:
Рис. 6. – Полученные программой графики.
Рис. 7. – Полученный программой график.
Текст программы и результаты метода Гуммеля:
clear all
close all
clc
warning off
cconstants; % инициализация необходимых физических констант
parameters; % инициализация шага на сетке, значения относительной погрешности к которой стремимся,
structure; % инициализация структуры, длины p и n областей, концентрации доноров и акцепторов в данных областях,
% шага для подгонки, напряжения приложенного к диоду, расчет концентраций
% собственых носителей заряда, а также параметров необходимых для расчета
% подвижности электронов и дырок
grid_functions; % задание сетки координат, концентраций и подвижностей с разным шагом для каждой области.
% загон конечной длины получившейся сетки в переменную I, а также
% инициализация матрицы разностей двух координат (x(2) - x(1)) и т.д.
normirovka; % нормировка координат, матрицы разностей, напряжения, концентрации, подвижностей
prog_param; % значения для толщин линий и размера шрифта, название зашитое в переменной
% для будущих файлов
fi0(1)=log(N(1)/2+sqrt((N(1)/2)^2+1));% граничные условия для потенциала. В начале
fi0(I)=-log(-N(I)/2+sqrt((N(I)/2)^2+1));% В конце
for i=2:I-1 %проход по всей сетке кроме граничных случаев, расчет электростатических
% потенциалов (далее электропотенциал) в ней
fi0(i)=fi0(1)+(fi0(I)-fi0(1))/(x(I)-x(1))*(x(i)-x(1));
end
fi0=fi0';% транспонирование
E=zeros(I,1);% инициализация матрицы единиц для напряженности электрического
% поля с I строк (кол-во точек в сетке) и 1 столбцом
for j=1:length(U) % проход по всем заданным напряжениям
if j>1 % условие для задания начальных значений, если напряжение является не 1 по списку
for i=1:I
fi0(i)=fi0(i)+(fi0(i)-fi0(1))/(fi0(I)-fi0(1))*(U(j)-U(j-1)); % перерасчет электропотенциалов
end
end
ct=0; % инициализация счетчика
err=100*Delta; % задание ошибки такой, чтоб начал работать цикл
while err>Delta
ct=ct+1; % увеличение счетчика на единицу
if U(j)==0 % проверка для определения 0 ли напряжение по списку, и если 0 то
Fn=ones(I,1); % инициализация матрицы для экспоненты квазиуровня
% ферми для электронов + задание начального значения
Fp=ones(I,1); % тоже самое но для дырок
else
A=zeros(I,I); % инициализация матрицы для частных производных Fn
B=zeros(I,1); % инициализация матрицы для значений функции
A(1,1)=1;
B(1)=1; % начальные условия
A(I,I)=1;
B(I)=exp(-U(j)); % конечные условия
for i=2:I-1
A(i,i)=-(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i)-...
(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); % расчет частной произврдная по Fn(i)
A(i,i+1)=(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i); % расчет частной производной по Fn(i+1)-Fn(i)
A(i,i-1)=(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fn(i)-Fn(i-1)
B(i)=0; % значение функции
end
Fn=A^(-1)*B; % решение уравнения
A=zeros(I,I);
B=zeros(I,1); % сброс матриц для использования в расчетах Fp
A(1,1)=1;
B(1)=1; % начальные условия
A(I,I)=1;
B(I)=exp(U(j)); % конечные условия
for i=2:I-1
A(i,i)=-(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2/dx(i)-...
(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fp(i)
A(i,i+1)=(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2/dx(i); % расчет частной производной по Fp(i+1)-Fp(i)
A(i,i-1)=(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fp(i)-Fp(i-1)
B(i)=0; % значение функции
end
Fp=A^(-1)*B; % решение уравнения
end
A=zeros(I,I);
B=zeros(I,1); % сброс матриц для последующего перерасчет fi
A(1,1)=1;
B(1)=log(N(1)/2+sqrt((N(1)/2)^2+1));
A(I,I)=1;
B(I)=-log(-N(I)/2+sqrt((N(I)/2)^2+1))+U(j);
for i=2:I-1
A(i,i)=-2/(dx(i)+dx(i-1))*(1/dx(i)+1/dx(i-1)); % расчет частной производной по fi(i)
A(i,i+1)=2/(dx(i)+dx(i-1))/dx(i); % расчет частной производной по fi(i+1)
A(i,i-1)=2/(dx(i)+dx(i-1))/dx(i-1); % расчет частной производной по fi(i-1)
B(i)=Fn(i)*exp(fi0(i))-Fp(i)*exp(-fi0(i))-N(i); % расчет значения функции
end
fi=A^(-1)*B; % решение уравнения
clc % очистка командного окна
Voltage=U(j)*V0 % денормирование напряжения, вывод на экран
err=max(abs(fi-fi0))/max(abs(fi0)) % пересчет невязки
fi0=fi0+(fi-fi0)./K(j); % расчет новых начальных условий
n=Fn.*exp(fi0); % расчет концентрации электронов
p=Fp.*exp(-fi0); % расчет концентрации дырок
for i=1:I-1
E(i)=-(fi0(i+1)-fi0(i))/dx(i); % расчет напряженности электрического поля
end
E(I)=E(I-1)+dx(I-1)/dx(I-2)*(E(I-1)-E(I-2)); % так как не можем посчитать последнее значение для
% поля ввиду отсутсвия следующего значения потенциала - принимаем
% его за последний рассчитанный
Ro=p-n+N; % общая концентрация зарядов
end
for i=2:I-1
Jn(i)=mun(i)*exp(fi0(i))*(Fn(i+1)-Fn(i-1))/(dx(i)+dx(i-1)); % расчет плотности тока электронов
Jp(i)=-mup(i)*exp(-fi0(i))*(Fp(i+1)-Fp(i-1))/(dx(i)+dx(i-1)); % расчет плотности ток дырок
end
Jn(I)=Jn(I-1);
Jp(I)=Jp(I-1);
Jn(1)=Jn(2);
Jp(1)=Jp(2);
J=Jn+Jp;
Current(j)=J(I);
graphics;
end
figure % создание отдельного окна
plot(U.*V0,-Current.*(J0*1e-4),'LineWidth',LW) % построение ВАХ
xlabel('Voltage, V','FontSize',FS) % именуем ось абсцисс
ylabel('Current density, A/cm^2','FontSize',FS) % именуем ось ординат
grid on % включаем сетку на графике
print(gcf,'-djpeg','VAH') % сохраняем график ВАХ в виде файла джипег
Результаты расчетов:
Voltage = 0.8000
err = 9.9987e-07
Полученные графики:
Рис. 8. – Полученные программой графики.
Рис. 9. – Полученный программой график.