В конечных разностях первую производную по времени в уравнении (26) запишем
@ , (29)
а вторые производные по пространственным координатам представим в виде
@ и (30)
@ . (31)
В (29) – (31) посредством τ обозначен шаг по времени, а h1 и h2 – по пространственным координатам и соответственно. С учетом (29) – (31) уравнение (26) с граничными и начальными условиями (27) – (28) примет вид
=D ( + ), (32)
c0i,j=0, i= , j= , (33)
ck0,j=f(j h2), k= , j= , (34)
где посредством |...| обозначена целая часть числа, а Tп – время диффузии.
Основной проблемой при реализации численного решения (32) – (34) является выбор величин τ, h1 и h2, то есть проблема устойчивости решения. Для уравнения диффузии, как показано в работе [5], эта задача решается выполнением следующих условий:
и (35)
. (36)
Рис. 2. Шаблон решения задачи (26) – (28)
С целью уменьшения времени вычислительного процесса предположим, что , а сами величины (35) – (36) –равными 0.25. В этом случае в уравнении (32) слагаемые, содержащие , уничтожатся и выражения (32) – (34) примут вид
=0.25 ( + ), (37)
c0i,j=0, i= , j= , (38)
ck0,j=f(j h), k= , j= , (39)
Таким образом, чтобы определить концентрацию диффузанта в узле с пространственными координатами ( , ) в момент времени (k + 1)∙τ, необходимо лишь просуммировать концентрации диффузанта в сопряженных узлах в момент времени kτ и полученную сумму разделить на четыре.
В качестве иллюстрации на рис. 3 показаны решения задачи (26) – (28) в моменты времени t=0, t, 2tи3t при следующих допущениях:
L1/h=L2/h=4;
f(j h)=const=1.6 – постоянный источник;
Tп /t=4.
| k=0 а. | k=1 б. |
| k=2 в. | k=3 г. |
Рис. 3. Решения задачи (26) – (28)
а. – t=0 (начальные и граничные условия), б. – t=t (первый шаг по времени), в. – t=2t(второй шаг по времени) иг. – t=3t(третий шаг по времени)
§3 Решение уравнения диффузии (разбор полетов)
Основная цель, которую автор поставил перед собой в предыдущем разделе, заключалась в том, чтобы показать читателю легкость, с которой численными методами может быть получено решение уравнения диффузии. Надеюсь, что цель достигнута и у читателя не возникли какие-либо трудности с восприятием алгоритма решения задачи о диффузии, а сам читатель стал в ряды сторонников решения ее численными методами. Дело все в том, что к сегодняшним практическим задачам полупроводниковой электроники [6] имеющиеся классические аналитические решения уравнения диффузии имеют весьма отдаленное отношение, поскольку в основу последних в большинстве своем заложены макросистемы, например, бесконечное или полубесконечное тело, источники диффузанта с неограниченной емкостью и т.д., что, скажем, не имеет места в технологии производства современных интегральных схем, где размеры каждого элемента соизмеримы с линейными параметрами таких величин, как дебаевская длина экранирования и толщина поверхностной фазы. Более того, явления, происходящие в областях, соизмеримых с дебаевской длиной экранирования, играют существенную роль и в процессе эксплуатации приборов, созданных на современной элементной базе. Например, актуальнейшей проблемой нынешней электроники является деградация электрофизических параметров приборов. Представляется, что решение указанной задачи в первую очередь напрямую связано и с пониманием диффузионных процессов, обусловленных электрическими полями, возникающими из-за наличия контакта разнородных материалов. Поэтому далее в настоящей работе мы получим уравнение диффузии, учитывающее влияние внутренних электрических полей на процессы переноса вещества, и покажем некоторые его решения.
Ну а сейчас вернемся к началу настоящего параграфа и обсудим, за счет чего же была достигнута простота решения уравнения диффузии? Во-первых, за счет сильного упрощения (идеализации) самой задачи, а, во-вторых, не обошлось и без некоторого лукавства. Если еще раз более внимательно проанализировать предыдущий раздел, то легко можно обнаружить, что и описание задачи (26)-(28), и рис. 1, и решения, представленные на рис. 3, содержат некоторые недомолвки - "параметры по умолчанию", но при этом вроде бы, дополняя друг друга, создают образ единой задачи. А, по сути, в (26)-(28) просто отсутствуют граничные условия за исключением области окна в маске. Рисунок 1 несколько восполняет этот недостаток своей волнистой обрамляющей, намекающей на бесконечную протяженность объекта исследования, и словом "маска", говорящему всякому, кто имеет хотя бы какое-либо отношение к технологии полупроводниковых приборов, что поток диффузанта через маску, по крайней мере, значительно ограничен по сравнению с областью окна или даже равен нулю. Предположив, что концентрация диффузанта в области маски всегда равна нулю, мы, решая уравнение диффузии, вроде бы использовали последнее более сильное условие. И в результате решили … совершенно не ту задачу. Оказывается, что область маски в решенной задаче не только препятствует проникновению диффузанта в кристалл, но и одновременно является местом стока диффузанта неограниченной емкости. Действительно, обратившись к рисунку 3, видим, что изоконцентраты закреплены у краев маски и с течением времени упорно "не хотят" перемещаться вдоль границы "маска – кристалл" несмотря на то, что градиенты концентрации у краев маски в этих направлениях, а, следовательно, и движущие силы диффузии, максимальны. А ларчик тут открывается просто. Поговорить то мы поговорили об отсутствии потоков диффузанта через маску, но ничего для того, чтобы их не было, не предприняли.
В целом наши проблемы связаны с тем, что мы неверно решали задачу. Подход к методике решения настолько плох, что, как я надеюсь, читатель и сам уже пришел к выводу, так задачи решать нельзя. За многие тысячелетия человечество сформулировало подходы, методики, алгоритмы решения проблем. Они просты и очевидны и, наверное, именно поэтому мы зачастую ими не пользуемся, не придаем им должного значения и как следствие затрачиваем огромные усилия на то, что делать и не следовало бы, но, к сожалению, об этом узнаем, когда работа уже выполнена или не узнаем вовсе. Современные представления о том, как подходить к решению той или иной проблемы можно найти в прекрасной книге профессора МИСиС М.А. Штремеля [7] и в литературе, указанной в библиографическом списке к ней. Для тех, кто занимается моделированием можно порекомендовать [8]. И, тем не менее, очень кратко сформулируем основные этапы решения задач, аналогичных решенной нами:
1. Постановка задачи и определение конечных целей – заключается в содержательной (физической) постановке задачи: выбор общего подхода, определение совокупности критериев, которым должна удовлетворять система, задание условий ее реализации (работы), выбор системы координат и масштабов, выбор системы единиц.
2. Построение математической модели – математическая формулировка задачи. Модель должна правильно (адекватно) описывать основные законы физики, функционирования объекта (системы, устройства), влияние параметров на его свойства и т.д.
3. Разработка (выбор) численного метода. Математическая формулировка задачи может оказаться непереводимой непосредственно на язык ЭВМ. Для того, чтобы задача была правильно решена, необходимо уравнения, функции, интегралы выразить через элементарные функции. Более того, необходимо убедиться, что никакие погрешности, содержащиеся в исходных данных или внесенные в процессе вычислений, не повлияют сколько-нибудь заметным образом на точность результатов.
4. Разработка алгоритма – выражение решения задачи в виде точно определенной последовательности операций вычислительной машины. Раньше алгоритм представляли графически в виде блок-схемы. Сегодня от графического представления чаще всего отказываются (особенно в случае больших задач), но словесный портрет задачи в виде последовательности операций с указанием условий переходов от пункта к пункту записывают всегда.
5. Программирование – изложение алгоритма на языке "понятном" ЭВМ.
6. Отладка программы. Очень важный этап в процессе решения задачи. С этой целью используют тестовые задачи, пошаговое выполнение программы.
7. Выполнение расчетов. Если это возможно, то выполняют для нескольких наборов данных.
8. Анализ и оформление результатов.
Самое главное в этом процессе заключается в том, что на каждом этапе необходимо постоянно обращаться к предыдущим и, в особенности к первому, задавая себе вопрос "А ту ли задачу я решаю?"
§4 Решение уравнения диффузии (продолжение)
Задача 1. Полное повторение предыдущей задачи с добавлением условия отсутствия потоков на границе "маска-кристалл" и уточненными начальными и граничными условиями:
= (40)
с начальными и граничными условиями:
c(x,y,0)=0, xÎ(0,+L1), yÎ[-L2,+L2], (41)
c(0,y,0)=0, (42)
c(0,y,t)=f(y), yÎ(-l, +l), (43)
(44)
Решение (изоконцентраты диффузанта) задачи (40) – (44) показано на рис. 5 (для сравнения на рис. 4 изображено решение этой же задачи, но без учета (44)), а все необходимые замечания к решению представлены в комментариях в программе для ЭВМ и к программе.
Решение задачи 1 без учета (44)
Рис. 4. Решение задачи 1 без учета (44)
Рис. 5. Решение задачи 1
Программа для решения задачи 1
01 # include <stdio.h>
02 # define SizeX (23) // Размер диффузионной структуры в направлении оси X
03 # define SizeZ (12) // Размер диффузионной структуры в направлении оси Z
04 # define SizeT (12) // Количество шагов по времени
05 # define SizeXWin (5) // Размер окна в направлении оси X
06 int main(void)
07 {
08 int iZ,iX,k;
09 double *Ck, *Ck1;
10 FILE *dif_dat;
11 dif_dat = fopen(" … :\\ … \\ … \\ dif_dat.dat", "w"); // Открытие файла для записи
// данных
12 int aX=(SizeX-SizeXWin)/2; // Координата начала окна в направлении оси X
13 int bX=aX+SizeXWin; // Координата конца окна в направлении оси X
14 Ck=new double [300]; // Выделение памяти для хранения данных на i-том шаге по
// времени
15 Ck1=new double [300]; // Выделение памяти для хранения данных на (i+1)-м шаге
// по времени
16 for (iX=0; iX<(aX-1); iX++) { // Граничные условия в области маски слева от окна
17 *(Ck+iX)=0;
18 *(Ck1+iX)=0;
19 }
20 for (iX=(aX-1); iX<(bX+1); iX++) { // Граничные условия в области окна
21 *(Ck+iX)=16;
22 *(Ck1+iX)=16;
23 }
24 for (iX=(bX+1); iX<SizeX; iX++) { // Граничные условия в области маски справа от окна
25 *(Ck+iX)=0;
26 *(Ck1+iX)=0;
27 }
28 for (iZ=1; iZ<SizeZ; iZ++) { // Начальные условия в области кристалла
29 for (iX=0; iX<SizeX; iX++) {
30 *(Ck+iX+iZ*SizeX)=0;
31 *(Ck1+iX+iZ*SizeX)=0;
32 }
33 }
// Решение уравнения
34 for (k=0; k<SizeT; k++){
35 for (iZ=1; iZ<(SizeZ-1); iZ++)
36 for (iX=1; iX<(SizeX-1); iX++)
37 *(Ck1+iX+iZ*SizeX)=*(Ck+iX+iZ*SizeX)+0.25*(*(Ck+iX+
(iZ+1)*SizeX)+(*(Ck+iX+(iZ-1)*SizeX))+
(*(Ck+iZ*SizeX+iX+1))+(*(Ck+iZ*SizeX+
iX-1))-(*(Ck+iX+iZ*SizeX))*4);
38 for (iZ=0; iZ<SizeZ; iZ++) // Копирование массива данных
39 for (iX=0; iX<SizeX; iX++)
40 *(Ck+iX+iZ*SizeX)=*(Ck1+iX+iZ*SizeX);
41 for (iX=0; iX<(aX-1); iX++) // Выполнение условия (44) слева от окна
42 *(Ck+iX)=*(Ck1+iX+SizeX);
43 for (iX=(bX+1); iX<SizeX; iX++) // Выполнение условия (44) справа от окна
44 *(Ck+iX)=*(Ck1+iX+SizeX);
45 }
46 for (iZ=0; iZ<SizeZ; iZ++){ // Сохранение результатов счета
47 for (iX=0; iX<SizeX; iX++)
48 fprintf(dif_dat," %7.3f",*(Ck+iX+iZ*SizeX));
49 fprintf(dif_dat,"\n");
50 }
51 delete Ck; // Освобождение выделенной памяти
52 delete Ck1;
53 fclose(dif_dat); // Закрытие файла
56 return 0;
57 }
Комментарии к программе:
· Нумерация строк (цифры слева) не является атрибутом языка программирования СИ.
· По сравнению с рисунком 1 в программе изменен выбор наименования осей системы координат, а само начало координат перенесено для того, чтобы отчетливей показать последовательность решения задачи. Ось X заменена осью Z, а ось Y – на X. Начало координат перенесено в левый верхний угол рисунка. В действительности выбор расположения начала координат так, как представлено на рис. 1, предпочтительнее в силу симметрии задачи относительно плоскости Y=0, что в свою очередь дает возможность решать "половину" задачи, сокращая при этом объем выделенной оперативной памяти ЭВМ и время счета соответственно в два раза. Правда, при этом необходимо несколько изменить граничные условия, не затрагивая существа задачи. Мы этого делать ни в данном случае, ни в следующих задачах не станем, оставляя возможность читателю разрешить указанную проблему самостоятельно, с которой он легко справится, внимательно проанализировав представленную выше программу.
· В строке 11 открывается файл для записи результатов счета. Если читатель захочет выполнить расчеты с использованием предложенной программы, то ему необходимо вместо многоточий указать место хранения файла (путь) применительно к конкретной ЭВМ.
· Следует обратить внимание на строку 37, где используется не частный случай представления уравнения диффузии в конечных разностях, формула (37), а более общий случай, который следует из (32) при h1 = h2 = h:
= + ( + ).
При этом численное значение величины D∙τ / h2, равное 0.25, сохранено.
· В строках 41 – 44 реализовано выполнение условия отсутствия потока на границе "маска – кристалл". Концентрация примеси во всех узлах в области маски приравнивается к концентрации примеси в области кристалла в узлах, расположенных непосредственно под первыми, то есть c(xi, 0)= c(xi,1), где i – номера узлов в области маски вдоль оси X -ов.
Имея ввиду, как отмечалось выше, малость объектов современной электроники и как следствие необходимость учета геометрической формы этих объектов и явлений, сопровождающих диффузионные процессы в очень тонких слоях, рассмотрим задачу о диффузии, где будет использован тот факт, что в различных материалах коэффициент диффузии одной и той же примеси неодинаковый и на границах раздела фаз может быть в десятки раз больше [1], чем в объеме материала. Кроме того, в отличие от предыдущей задачи, заменим источник диффузанта бесконечной емкости на источник ограниченной емкости (например, диффузия из силикатного стекла).
Задача 2. Решить уравнение диффузии
= (45)
где
D1 – коэффициент диффузии в области диффузанта,
D2 – коэффициент диффузии на границе "диффузант-маска",
D3 – коэффициент диффузии в области кристалла,
D4 – коэффициент диффузии на границе "маска-кристалл",
с начальными и граничными условиями:
c(x,0,t)=0, xÎ[0,L1], tÎ[0, Tпр.], (46)
c(x,z,0)=c0д, xÎ[0,L1], zÎ(0,(L2-l2)*0.5], (47)
c(x,z,0)=c0д, xÎ[(L1-l1)*0.5, (L1+l1)*0.5), zÎ((L2-l2)*0.5,(L2+l2)*0.5], (48)
c(x,(L2+l2)*0.5,0)=0, xÎ[0, (L1-l1)*0.5), (49)
c(x,(L2+l2)*0.5,0)=0, xÎ[(L1+l1)*0.5), L1), (50)
c(x,z,0)=0, xÎ[0,L1], zÎ((L2+l2)*0.5), L2] (51)
zÎ[0, L2],(52)
zÎ[0, L2], (53)
zÎ((L2-l2)*0.5,(L2+l2)*0.5], (54)
zÎ((L2-l2)*0.5,(L2+l2)*0.5], (55)
xÎ[0,(L1-l1)*0.5], (56)
xÎ[0,(L1-l1)*0.5], (57)
xÎ[(L1+l1)*0.5,L1], (58)
xÎ[(L1+l1)*0.5,L1], (59)
где
Tпр. – время диффузии,
L1 – линейный размер структуры в направлении оси x,
L2 – линейный размер структуры в направлении оси Z,
l1 – линейный размер щели маски (ширина окна) в направлении оси x,
l2 – линейный размер маски (толщина маски) в направлении оси Z,
c0д – начальная концентрация примеси в диффузанте.
Схема процесса диффузии показана на рисунке 6.
Решение (изоконцентраты диффузанта) задачи 2 показано на рис. 8. Для сравнения на рис. 7 изображено решение аналогичной задачи, но без учета повышенной скорости диффузии на границах раздела, а на рис. 9 – решение уравнения диффузии с учетом повышенной скорости диффузии на границах раздела при наличии в кристалле незаряженной стенки. В комментариях полученные решения, как мне представляется, не нуждаются, но все же следует напомнить должно быть известное читателю из курса общей физики положение о том, что повышенная плотность изоконцентрат свидетельствует о более высоких градиентах концентрации и следовательно о повышенных скоростях диффузии в этих областях. На рисунках 7-9 не показаны решения в области диффузанта выше маски, поскольку для данных задач они не представляют значительного интереса (изоконцентраты – прямые линии, которые располагаются перпендикулярно оси Z) и при этом существенно увеличивают объем рисунков. Все необходимые замечания к решению представлены в комментариях в программе для ЭВМ и к программе.
Рис. 6. Схема процесса диффузии. Задача 2
Рис. 7. Решение задачи 2 без учета повышенной скорости диффузии на границах раздела
Рис. 8. Решение задачи 2
Рис. 9. Решение задачи, аналогичной задаче 2 при наличии стенки в кристалле в области окна
Программа для решения задачи 2
01 # include <stdio.h>
02 # define SizeX (46) // Размер диффузионной структуры в направлении оси X
03 # define SizeZ (46) // Размер диффузионной структуры в направлении оси Z
04 # define SizeT (128) // Количество шагов по времени
05 # define SizeXWin (14) // Размер окна в направлении оси X
06 # define SizeZWin (6) // Размер окна в направлении оси Z
07 # define C0d (104) // Начальная концентрация примеси в диффузанте
08 # define C0k (0) // Начальная концентрация примеси в кристалле
09 # define DtauD (0.625) // Параметр DtauD=D1*tau/(h*h) в области диффузанта
10 # define DtauDW (0.0625) // Параметр DtauDW=D2*tau/(h*h) на границе "диффу-
// зант - маска"
11 # define DtauK (0.025) // Параметр DtauK=D3*tau/(h*h) в области кристалла
12 # define DtauKW (0.25) // Параметр DtauKW=D4*tau/(h*h) на границе "маска-
// кристалл"
13 int main(void)
14 {
15 int iZ,iX,k;
16 double *Ck, *Ck1;
17 FILE *dif_dat;
18 dif_dat = fopen(" … :\\ … \\ … \\ dif_dat.dat", "w");
19 if (dif_dat == NULL){
20 printf("Файл dif_dat.dat не открыт.\n");
21 return 1;
22 }
[1] От латинского diffusio - распространение, растекание.
[2] Векторные величины в работе выделены "жирным" шрифтом.
[3] Здесь и далее для краткости записи аргументы функций опущены за исключением тех случаев, когда это может, по мнению автора, привести к недоразумениям.
[4] Во всех последующих формулах посредством kБ обозначена постоянная Больцмана