Уравнения диффузии при наличии электрического поля
I-том шаге по времени
26 Ck1=new double [(SizeX+1)*(SizeZ+1)]; // Выделение памяти для хранения данных
// на (i+1)-м шаге по времени
27 // Начальные и граничные условия
28 for (iX=0; iX<SizeX+1; iX++) { // Граница "вакуум (среда) – диффузант"
29 *(Ck+iX)=0;
30 *(Ck1+iX)=0;
31 }
32 for (iZ=1; iZ<aZ+1; iZ++) { // Область диффузанта
33 for (iX=0; iX<SizeX+1; iX++) {
34 *(Ck+iZ*(SizeX+1)+iX)=C0d;
35 *(Ck1+iZ*(SizeX+1)+iX)=C0d;
36 }
37 }
38 for (iZ=aZ+1; iZ<bZ+1; iZ++) { // Область окна
39 for (iX=aX; iX<bX+1; iX++) {
40 *(Ck+iZ*(SizeX+1)+iX)=C0d;
41 *(Ck1+iZ*(SizeX+1)+iX)=C0d;
42 }
43 }
44 for (iZ=aZ+1; iZ<bZ+1; iZ++) { // Левая область маски
45 for (iX=0; iX<aX; iX++) {
46 *(Ck+iZ*(SizeX+1)+iX)=C0k;
47 *(Ck1+iZ*(SizeX+1)+iX)=C0k;
48 }
49 }
50 for (iZ=aZ+1; iZ<bZ+1; iZ++) { // Правая область маски
51 for (iX=bX+1; iX<SizeX+1; iX++) {
52 *(Ck+iZ*(SizeX+1)+iX)=C0k;
53 *(Ck1+iZ*(SizeX+1)+iX)=C0k;
54 }
55 }
56 for (iZ=bZ+1; iZ<SizeZ+1; iZ++) { // Область кристалла
57 for (iX=0; iX<SizeX+1; iX++) {
58 *(Ck+iZ*(SizeX+1)+iX)=C0k;
59 *(Ck1+iZ*(SizeX+1)+iX)=C0k;
60 }
61 }
62 // Расчет диффузии
63 for (k=0; k<SizeT; k++){
64 for (iZ=1; iZ<aZ-1; iZ++) // Область диффузанта
65 for (iX=1; iX<SizeX; iX++)
*(Ck1+iX+iZ*(SizeX+1))=*(Ck+iX+iZ*(SizeX+1))+DtauD*
(*(Ck+iX+(iZ+1)*(SizeX+1))+(*(Ck+iX+(iZ-1)*(SizeX+1)))+
(*(Ck+iZ*(SizeX+1)+iX+1))+(*(Ck+iZ*(SizeX+1)+iX-1))-
(*(Ck+iX+iZ*(SizeX+1)))*4);
66 for (iX=1; iX<aX+1; iX++) // Область диффузанта, левая граница
// "диффузант-маска"
*(Ck1+iX+(aZ-1)*(SizeX+1))=*(Ck+iX+(aZ-1)*(SizeX+1))+DtauDW*
(*(Ck+iX+aZ*(SizeX+1))+(*(Ck+iX+(aZ-2)*(SizeX+1)))+(*(Ck+
(aZ-1)*(SizeX+1)+iX+1))+(*(Ck+(aZ-1)*(SizeX+1)+iX-1))-(*(Ck+
iX+(aZ-1)*(SizeX+1)))*4);
67 for (iX=aX+1; iX<bX; iX++) // Область диффузанта, центр (продолжение
// строки 66)
*(Ck1+iX+(aZ-1)*(SizeX+1))=*(Ck+iX+(aZ-1)*(SizeX+1))+DtauD*
(*(Ck+iX+aZ*(SizeX+1))+(*(Ck+iX+(aZ-2)*(SizeX+1)))+(*(Ck+
(aZ-1)*(SizeX+1)+iX+1))+(*(Ck+(aZ-1)*(SizeX+1)+iX-1))-(*(Ck+
iX+(aZ-1)*(SizeX+1)))*4);
68 for (iX=bX; iX<SizeX; iX++) // Область диффузанта, правая граница
"диффузант-маска" (продолжение строки 67)
*(Ck1+iX+(aZ-1)*(SizeX+1))=*(Ck+iX+(aZ-1)*(SizeX+1))+DtauDW*
(*(Ck+iX+aZ*(SizeX+1))+(*(Ck+iX+(aZ-2)*(SizeX+1)))+(*(Ck+
(aZ-1)*(SizeX+1)+iX+1))+(*(Ck+(aZ-1)*(SizeX+1)+iX-1))-
(*(Ck+iX+(aZ-1)*(SizeX+1)))*4);
69 for (iZ=aZ; iZ<bZ+1; iZ++) // Область окна, левая вертикальная граница
// "диффузант-маска"
*(Ck1+aX+1+iZ*(SizeX+1))=*(Ck+aX+1+iZ*(SizeX+1))+DtauDW*
(*(Ck+aX+1+(iZ+1)* (SizeX+1))+(*(Ck+aX+1+(iZ-1)*(SizeX+1)))+
(*(Ck+iZ*(SizeX+1)+aX+2))+(*(Ck+iZ*(SizeX+1)+aX))-(*(Ck+aX+
1+iZ*(SizeX+1)))*4);
70 for (iZ=aZ; iZ<bZ+1; iZ++) // Область окна, правая вертикальная граница
// "диффузант-маска"
*(Ck1+bX-1+iZ*(SizeX+1))=*(Ck+bX-1+iZ*(SizeX+1))+DtauDW*
(*(Ck+bX-1+(iZ+1)*(SizeX+1))+(*(Ck+bX-1+(iZ-1)*(SizeX+1)))+
(*(Ck+iZ*(SizeX+1)+bX))+(*(Ck+iZ*(SizeX+1)+bX-2))-(*(Ck+
bX-1+iZ*(SizeX+1)))*4);
71 for (iZ=aZ; iZ<bZ+1; iZ++) // Область окна
for (iX=aX+2; iX<bX-1; iX++)
*(Ck1+iX+iZ*(SizeX+1))=*(Ck+iX+iZ*(SizeX+1))+DtauD*
(*(Ck+iX+(iZ+1)*(SizeX+1))+(*(Ck+iX+(iZ-1)*(SizeX+1)))+
(*(Ck+iZ*(SizeX+1)+iX+1))+(*(Ck+iZ*(SizeX+1)+iX-1))-
(*(Ck+iX+iZ*(SizeX+1)))*4);
72 for (iX=1; iX<aX+1; iX++) // Область кристалла, левая граница
// "маска-кристалл"
*(Ck1+iX+(bZ+1)*(SizeX+1))=*(Ck+iX+(bZ+1)*(SizeX+1))+DtauKW*
(*(Ck+iX+(bZ+2)*(SizeX+1))+(*(Ck+iX+bZ*(SizeX+1)))+(*(Ck+(bZ+1)*
(SizeX+1)+iX+1))+(*(Ck+(bZ+1)*(SizeX+1)+iX-1))-
(*(Ck+iX+(bZ+1)*(SizeX+1)))*4);
73 for (iX=aX+1; iX<bX; iX++) // Область кристалла, центр (продолжение
// строки 72)
*(Ck1+iX+(bZ+1)*(SizeX+1))=*(Ck+iX+(bZ+1)*(SizeX+1))+DtauK*
(*(Ck+iX+(bZ+2)*(SizeX+1))+(*(Ck+iX+bZ*(SizeX+1)))+(*(Ck+
(bZ+1)*(SizeX+1)+iX+1))+(*(Ck+(bZ+1)*(SizeX+1)+iX-1))-
(*(Ck+iX+(bZ+1)*(SizeX+1)))*4);
74 for (iX=bX; iX<SizeX; iX++) // Область кристалла, правая граница
// "маска-кристалл"
*(Ck1+iX+(bZ+1)*(SizeX+1))=*(Ck+iX+(bZ+1)*(SizeX+1))+DtauKW*
(*(Ck+iX+(bZ+2)*(SizeX+1))+(*(Ck+iX+bZ*(SizeX+1)))+(*(Ck+
(bZ+1)*(SizeX+1)+iX+1))+(*(Ck+(bZ+1)*(SizeX+1)+iX-1))-(*(Ck+
iX+(bZ+1)*(SizeX+1)))*4);
75 for (iZ=bZ+2; iZ<SizeZ; iZ++) // Область кристалла
76 for (iX=1; iX<SizeX; iX++)
*(Ck1+iX+iZ*(SizeX+1))=*(Ck+iX+iZ*(SizeX+1))+DtauK*
(*(Ck+iX+(iZ+1)*(SizeX+1))+(*(Ck+iX+(iZ-1)*(SizeX+1)))+
(*(Ck+iZ*(SizeX+1)+iX+1))+(*(Ck+iZ*(SizeX+1)+iX-1))-
(*(Ck+iX+iZ*(SizeX+1)))*4);
77 for (iZ=0; iZ<SizeZ+1; iZ++) // Копирование данных (i+1)-го шага по
Времени в массив i-го
78 for (iX=0; iX<SizeX+1; iX++)
79 *(Ck+iX+iZ*(SizeX+1))=*(Ck1+iX+iZ*(SizeX+1));
80 for (iX=0; iX<aX; iX++) { // Выполнения требования отсутствия потоков
// на левой:
81 *(Ck+iX+aZ*(SizeX+1))=*(Ck1+iX+(aZ-1)*(SizeX+1)); // верхней
// границе "диффузант-маска",
82 *(Ck+iX+bZ*(SizeX+1))=*(Ck1+iX+(bZ+1)*(SizeX+1)); // нижней
// границе "маска-кристалл"
83 }
84 for (iX=bX+1; iX<SizeX+1; iX++) { // Выполнения требования отсутствия
// потоков на правой:
85 *(Ck+iX+aZ*(SizeX+1))=*(Ck1+iX+(aZ-1)*(SizeX+1)); // верхней
// границе "диффузант-маска",
86 *(Ck+iX+bZ*(SizeX+1))=*(Ck1+iX+(bZ+1)*(SizeX+1)); // нижней
// границе "маска-кристалл"
87 }
88 for (iZ=aZ; iZ<bZ+1; iZ++) { // Выполнения требования отсутствия
// потоков на вертикальной:
89 *(Ck+aX+iZ*(SizeX+1))=*(Ck1+aX+1+iZ*(SizeX+1)); // левой
// границе "диффузант-маска",
90 *(Ck+bX+iZ*(SizeX+1))=*(Ck1+bX-1+iZ*(SizeX+1)); // правой
// границе "диффузант-маска"
91 }
92 for (iZ=0; iZ<aZ; iZ++) { // Выполнения требования отсутствия потоков
// вдоль направления оси X-ов в области диффузанта на:
93 *(Ck+iZ*(SizeX+1))=*(Ck1+1+iZ*(SizeX+1)); // левой вертикальной
// границе
94 *(Ck+SizeX+iZ*(SizeX+1))=*(Ck1+SizeX-1+iZ*(SizeX+1)); // правой
// вертикальной границе
95 }
96 }
97 for (iZ=0; iZ<SizeZ+1; iZ++){ // Печать результатов счета в файл "dif_dat"
98 for (iX=0; iX<SizeX+1; iX++)
99 fprintf(dif_dat," %7.3f",*(Ck+iX+iZ*(SizeX+1)));
100 fprintf(dif_dat,"\n");
101 }
102 delete Ck; // Освобождение выделенной оперативной памяти
103 delete Ck1;
104 fclose(dif_dat);
105 return 0;
106 }
Комментарии к программе:
· Нумерация строк (цифры слева) не является атрибутом языка программирования СИ.
· В строке 18 открывается файл для записи результатов счета. Если читатель захочет выполнить расчеты с использованием предложенной программы, то ему необходимо вместо многоточий указать место хранения файла (путь) применительно к конкретной ЭВМ.
· В строках 28–30 задаются граничные условия на границе раздела "диффузант - среда". В реальных условиях эксперимента (и в нашей задаче тоже) эта граница является местом стока диффузанта за исключением случая, когда над поверхностью диффузанта специально поддерживают равновесное давление пара диффундирующей примеси.
· В строках 44–55 задаются начальные условия в области маски, которые необходимы в случае рассмотрения диффузионных процессов в этой области (диффузия через маску). Для нашей задачи строки 44, 47, 52, 53 можно исключить, не забыв при этом в строках 46, 47, 52, 53 произвести замену переменной iZ на bZ.
· В строке 66 выполняется решение уравнения диффузии на левой верхней границе раздела "диффузант – маска".
· Строка 67 является продолжением строки 66, но при этом выполняется решение уравнения диффузии в области окна.
· В строке 68 выполняется решение уравнения диффузии на правой верхней границе раздела "диффузант – маска". Эта строка является продолжением строки 67.
· В строке 72 выполняется решение уравнения диффузии на левой нижней границе раздела "маска – кристалл".
· Строка 73 является продолжением строки 72, но при этом выполняется решение уравнения диффузии в области окна.
· В строке 74 выполняется решение уравнения диффузии на правой нижней границе раздела "маска – кристалл". Эта строка является продолжением строки 73.
· В строках 89 – 90 осуществляется выполнение требования отсутствия потоков на вертикальных границах раздела "диффузант – маска" в области окна.
И вот теперь, после решения последней задачи, настала пора перейти изучению влияния электрических полей на диффузионные процессы. Прежде всего, кратко остановимся на выводе уравнения Пуассона, которое необходимо решать совместно с уравнением диффузии для достижения поставленной цели.
§5 Уравнение Пуассона
Рассчитаем поток вектора индукции электрического поля D, созданного положительным точечным зарядом q, через поверхность сферы радиуса R с центром, совпадающим с положением заряда в пространстве, рис.10.
По определению, поток вектора А через поверхность S рассчитывают, применяя интеграл по поверхности, следующим образом
(60)
где использовано общепринятое обозначение dSтаким же образом, как и в (8). В нашем случае
(61)
В системе единиц CGS вектор индукции электрического поля, создаваемого точечным зарядом,
(62)
В (62) ε- относительная диэлектрическая проницаемость среды, в которую помещен заряд, а E –вектор напряженности электрического поля. Поскольку вектор D
Рис.10. Схема расчета потока вектора электрической индукции
нормален во всех точках рассматриваемой поверхности и модуль его постоянен на этой поверхности, а R/R есть ни что иное, как n, то скалярное произведение в (61) равно
(63)
а поток вектора Dравен
(64)
То есть, поток вектора индукции электрического поля через замкнутую поверхность равен произведению 4π на величину заряда, размещенного в пространстве, ограниченном этой поверхностью. Полученное утверждение называют законом или теоремой Гаусса.
Некоторую неудовлетворенность у читателя при выводе (64) может вызвать простота изначальных посылок. Часто возникает вопрос: "А адекватно ли (64) будет выражение, полученное в случае выбора более сложной пространственной конфигурации ограничивающей поверхности?" При условии; что величина модуля вектора, поток которого мы вычисляем, изменяется в зависимости от расстояния пропорционально квадрату этого расстояния в минус первой степени, ответ на вопрос следует дать утвердительный. Действительно, всякую "сложную" поверхность можно представить в виде суммы ее элементарных участков, одни из которых располагаются перпендикулярно к радиус-векторам, начинающимся на заряде q (в этом случае радиус-вектор и нормаль к поверхности компланарны, то есть скалярное произведение R/R∙n равно единице), а с нормалями к другим поверхностям, указанные выше радиус-векторы, ортогональны (в этом случае скалярное произведение и, соответственно, потоки через эти поверхности равны нулю). На рис.11 показано сечение произвольной поверхности некоторой плоскостью, проходящей через центр заряда. Для первых элементарных поверхностей в сферической системе координат
(65)
а вклад в величину потока вторых поверхностей равен нулю в силу обстоятельств, отмеченных в предыдущем предложении. Более подробное исследование областей применения и доказательство теоремы Гаусса заинтересованный читатель может найти в учебной и специальной литературе по общей и теоретической физике и математике.
В случае заряда q, непрерывно распределенного с плотностью ρ(x,y,z)в объеме ΔV, ограниченном поверхностью ΔS, уравнение (64) преобразуется к виду
(66)
Воспользовавшись теоремой Остроградского-Гаусса, из (66) получим
(67)
Поскольку области интегрирования произвольны и одинаковы для левой и правой частей (67), то для выполнения равенства необходимо, чтобы
(68)
Введем, учитывая (62), в (68) вместо вектора индукции электрического поля D скалярную величину - потенциал электрического поля φ, которые связаны между собой
Рис. 11. Сечение произвольной поверхности
следующим образом
. (69)
Имея в виду, что
, (70)
получим из (70) с учетом (69) уравнение Пуассона для однородной и изотропной среды ( )
(71)
Очень важно обратить внимание на то, что при выводе (71) все рассуждения велись применительно к положительному заряду и только в этом случае в правой части уравнения Пуассона стоит знак минус. Более того, необходимо всегда иметь в виду, что ρ(x,y,z)это плотность заряда, а не концентрация заряженных частиц.
Сложностей с аналитическими решениями уравнения Пуассона [9] не меньше, чем с уравнением диффузии, а поскольку во всех интересующих нас задачах уравнение (71) не линейно (правая часть содержит экспоненциальные функции от электрического потенциала), то и при его численном решении возникает достаточно много трудностей. В учебниках и задачниках [10] аналитические решения уравнения Пуассона приводятся чаще всего для случаев, когда разность потенциалов между материалами, находящимися в контакте, такова, что аргумент в экспоненте меньше единицы, что, в свою очередь, позволяет, разложив в ряд правую часть, привести уравнение к решабельному виду. Проблема, однако, заключается в том, что на практике этот случай не всегда реализуется. Действительно, запишем уравнение Пуассона для решения задачи об одномерном пространственном распределении электрического потенциала между двумя бесконечными в двух других пространственных измерениях полупроводниках p и n типов с полностью ионизованными примесями, полагая при этом, что никакой деформации ионных подсистем полупроводников на границе их контакта не происходит,
(72)
(73)
(74)
где
e=4,8032 10-10 CGSEq (1,6022 10-19 Кл) – величина заряда электрона,
p0 – концентрация дырок в валентной зоне вдали от контакта, см-3,
n0 – концентрация электронов в зоне проводимости вдали от контакта, см-3,
NA – концентрация акцепторной примеси в полупроводнике p -типа, см-3,
ND – концентрация донорной примеси в полупроводнике n -типа, см-3,
– величина электрического потенциала в полупроводнике n –типа вдали от контакта, CGSEV (1 В = 3,34 10-3 CGSEV),
– величина электрического потенциала в полупроводнике p –типа вдали от контакта, CGSEV,
x0– координата металлургической границы между p и n областями, см.
Абсолютное значение величин показателей экспонент в (72) при комнатной температуре равно
(75)
то есть даже при контактной разности потенциалов в одну десятую Вольта составляет 1,16·104·0,1·3,34·10-3=3,87 относительных единиц, что почти в четыре раза больше единицы и, соответственно, о замене ex на (1+x) не может быть и речи, поскольку e3,855 ≈47, а (1+x)=4,855 (относительная ошибка составляет ~860 %). И только при eφ < kБT (для температуры в 300 К контактная разность потенциалов должна быть меньше 26 мВ) выражение ex ≈ 1+x справедливо.
Введя безразмерный потенциал
(76)
и полагая, что ni << NA < ND, получим из (72)-(74)
(77)
(78)
(79)
где
– дебаевская длина экранирования в полупроводнике n –типа, см,
– дебаевская длина экранирования в полупроводнике p -типа, см,
– концентрация собственных носителей в полупроводнике, см-3,
, относ. ед.
В конечных разностях (77)-(79) примут вид
(80)
(81)
(82)
Решить задачу (80)-(82) можно было бы методом стрельбы [11] при условии, что мы достаточно просто сможем определить границы существования (81)-(82), то есть точно указать координаты точек, между которыми происходит изменение электрического потенциала между p и n областями (другими словами необходимо определить координаты точек, где "заканчивается" -∞ и "начинается" +∞). Эта проблема в совокупности с нелинейностью уравнения является достаточно трудоемкой, поэтому, приняв за единицу измерения пространственных координат наименьшую из дебаевских длин экранирования (в нашем случае - ), перепишем (80) для метода итераций
(83)
с начальными условиями:
(84)
где и все величины в левых частях равенств относятся к (k+1)-ой итерации, а все величины в правых частях равенств – к k-ой. Далее представлена программа для решения (83)-(84) на ЭВМ в несколько расширенном для учебных целей варианте (предусмотрен вакуумный зазор между полупроводниками), которая решает следующую задачу:
(85)
, (86)
где
x1 – координата границы раздела "полупроводник p-типа – вакуум", см,
x2– координата границы раздела "полупроводник n-типа – вакуум", см,
– работа выхода электрона полупроводника p-типа, отн.ед.,
– работа выхода электрона полупроводника p-типа, отн.ед.
Решение задачи (83) – (84) является частным случаем решения (85) – (86), когда (x2-x1)→0.
Программа для решения уравнения Пуассона
01 # include <stdio.h>
02 # include <math.h>
03 #include <conio.h>
04 #include <graphics.h>
05 int main()
06 {
07 int i,j;
08 double v,d,M;
09 double *fi1; //Электрический потенциал
10 double *Е; //Напряженность электрического поля
11 double *ro; //Плотность заряда
12 double r=1; //масштаб ro по вертикали
13 double *h; //Величина шага
14 int size=1000; //Размер области интегрирования уравнения (количество узлов)
15 int l=2; //Масштабирование графиков на экране дисплея по горизонтали от-
//носитель но области с меньшей концентрацией заряда
16 double f=5; //масштаб fi1 по вертикали
17 double e=1; //масштаб D по вертикали
18 double m=0.01; //Масштабирование шага в области с большей концентрацией заряда
19 double eps=0.001; //Модуль величины максимальной разности абсолютных зна-
//чений электрического потенциала между k–ой и (k+1)–ой
//итерациями на всей области интегрирования
20 double h0=0.001; //Начальное значение величины шага в области с меньшей
//концентрацией заряда
21 double An=0.5; //
22 double border1=0.3; //x1в долях от size
23 double border2=0.5; //x2в долях от size
24 FILE *dif_dat;
25 dif_dat = fopen("…:\\…\\…\\dif_dat.dat", "w");
26 if (dif_dat == NULL){
27 printf("Файл dif_dat.dat не открыт.\n");
28 return 1;
29 }
30 int gdriver = DETECT, gmode, errorcode;
31 initgraph(&gdriver,&gmode,"…:\\…\\BGI\\");
32 errorcode = graphresult();
33 if (errorcode != grOk) {
34 printf("Graphics error: %s\n", grapherrormsg(errorcode));
35 printf("Press any key to halt:");
36 getch();
37 return 1;
38 }
39 fi1=new double [size];
40 Е=new double [size];
41 ro=new double [size];
42 h=new double [size];
43 for (j=0; j<size*border1; j++) { //Начальные условия
44 *(fi1+j)=3;
45 *(h+j)=h0;
46 }
47 for (j=size*border1; j<size*border2; j++) {
48 *(fi1+j)=0;
49 *(h+j)=h0;
50 }
51 for (j=size*border2; j<size; j++) {
52 *(fi1+j)=10;
53 *(h+j)=m*h0;
54 }
55 int scrY=getmaxy(); //Вывод на экран начальных условий ψ(x)
56 for (i=0; i<size; i+=l)
57 putpixel(ceil(i/l), scrY*0.30-ceil(f*(*(fi1+i))), 5);
58 i=0;
59 do {
60 M=0;
61 for (j=1; j<size*border1; j++) { //Расчет ψ(x) в области x<x1
62 v=(*(fi1+j+1)+(*(fi1+j-1))+(*(h+j))*(*(h+j))*An*
(expl(-(*(fi1+j)-3))-1))/2;
63 d=fabs(*(fi1+j)-v);
64 *(fi1+j)=v;
65 if (M<d)
66 M=d;
67 else if (d>eps) //Изменение величины шага
68 (*(h+j))*=(1-eps/d);
69 }
70 for (j=size*border1; j<size* (border1+border2)/2; j++) {//Расчет ψ(x) в области
x1≤x<(x1+x2)/2
71 v=(*(fi1+j+1)+(*(fi1+j-1)))/2;
72 d=fabs(*(fi1+j)-v);
73 *(fi1+j)=v;
74 if (M<d)
75 M=d;
76 else if (d>eps) //Изменение величины шага
77 (*(h+j))*=(1-eps/d);
78 for (j=size*border2; j<size-1; j++) { //Расчет ψ(x) в области x≥x2
79 v=(*(fi1+j+1)+(*(fi1+j-1))-(*(h+j))*(*(h+j))*(expl((*(fi1+j)-10))-1))/2;
80 d=fabs(*(fi1+j)-v);
81 *(fi1+j)=v;
82 if (M<d)
83 M=d;
84 else if (d>eps) //Изменение величины шага
85 (*(h+j))*=(1-eps/d);
86 }
87 i++;
88 } while (M>eps);
89 printf("%d\n",i);
90 for (j=1; j<11; j++) { // Сохранение рассчитанных величин ψ(x) в файл dif_dat
91 for (i=size*(j-1)/10; i<size*j/10; i++)
92 fprintf(dif_dat," %7.3f",*(fi1+i));
93 fprintf(dif_dat,"\n");
94 }
95 for (i=0; i<size-1; i++) //Расчет Е(x)
96 *(Е+i)=-(*(fi1+i+1)-(*(fi1+i)))/(*(h+j));
97 for (i=1; i<size-1; i++) //Расчет ρ(x)
98 *(ro+i)=-(*(fi1+i+1)-2*(*(fi1+i))+(*(fi1+i-1)))/(*(h+j))/(*(h+j));
99 setcolor(5);
100 for (i=0; i<size; i+=l) //Вывод на экран рассчитанных величин ψ(x)
101 putpixel(ceil(i/l), scrY*0.30-ceil(f*(*(fi1+i))), 2);
102 for (i=0; i<size-1; i+=l) //Вывод на экран рассчитанных величин D(x)
103 putpixel(ceil(i/l), scrY*0.5+ceil(e*(*(Е+i))), 3);
104 for (i=1; i<size-1; i+=l) //Вывод на экран рассчитанных величин ro(x)
105 putpixel(ceil(i/l), scrY*0.8+ceil(r*(*(ro+i))), 4);
106 getch();
107 closegraph();
108 delete fi1;
109 delete E;
110 delete ro;
111 fclose(dif_dat);
112 return 0;
113 }
В тексте программы даны достаточно подробные комментарии, поэтому в дополнительных нет никакой необходимости.
На рисунках 12 – 15 показаны решения уравнения Пуассона для случаев: (x2-x1)= 0.5*size, (x2-x1)=0.3*size, (x2-x1)= 0.2*size и (x2-x1)= 0.0 соответственно. Следует отметить, что на всех рисунках правее вертикальных стрелок величина шага по пространственной координате примерно в 10 раз меньше, чем в областях до стрелок. На рисунке 12 представлен случай, когда расстояние между материалами велико, то есть между ними нет взаимодействия. Поэтому оба образца ведут себя как самостоятельные объекты. На границах раздела "материал – вакуум" каждого из них образуются по две заряженные стенки разного знака, причем внешние относительно объемов этих объектов слои отрицательно заряжены, так как они состоят из электронов вне зависимости от типа проводимости полупроводников. По мере сближения материалов в результате их взаимодействия происходит деформация заряженных слоев, рис. 13 – 14, при этом абсолютные значения величин зарядов в этих слоях уменьшаются. И, наконец, на рисунке 15 показан результат полного контакта двух полупроводников разного типа проводимости (без учета деформации пограничных слоев каждого из них).
Рис. 12. Решение уравнения Пуассона при (x2-x1) = 0.5*size
Рис. 13. Решение уравнения Пуассона при (x2-x1) = 0.3*size
Рис. 14. Решение уравнения Пуассона при (x2-x1) = 0.2*size
Рис. 15. Решение уравнения Пуассона при (x2-x1) = 0.0
Уравнения диффузии при наличии электрического поля
Для того, чтобы получить уравнение диффузии, которое учитывало бы влияние электрических полей в материале на диффузионные процессы, достаточно в соответствующих уравнениях §1 добавить к химическому потенциалу μi заряженной частицы qi, находящейся в электрическом поле с потенциалом φ, ее энергию, обусловленную этим полем,
. (87)
В этом случае уравнение (7) примет вид
(88)
или с учетом (19)
, (89)
а уравнение (13) преобразуется в выражение
. (90)
Имея в виду (62), (69) и (71) и полагая, что все моменты, связанные с неидеальностью системы, перенесены на коэффициент диффузии, который при этом не зависит от пространственных координат, из (89) и (90) получим соответственно первый и второй законы Фика в присутствии электрического поля:
, (91)
. (92)
Представляется, не лишним будет напоминание о том, что ci ≡ ci(x,y,z,t), E ≡ E(x,y,z,t) и ρ ≡ ρ(x,y,z,t), то есть для краткости записи у переменных величин