В конечных разностях первую производную по времени в уравнении (26) запишем

Основная цель, которую автор поставил перед собой в предыдущем разделе, заключалась в том, чтобы показать читателю легкость, с которой численными методами может быть получено решение уравнения диффузии. Надеюсь, что цель достигнута и у читателя не возникли какие-либо трудности с восприятием алгоритма решения задачи о диффузии, а сам читатель стал в ряды сторонников решения ее численными методами. Дело все в том, что к сегодняшним практическим задачам полупроводниковой электроники [6] имеющиеся классические аналитические решения уравнения диффузии имеют весьма отдаленное отношение, поскольку в основу последних в большинстве своем заложены макросистемы, например, бесконечное или полубесконечное тело, источники диффузанта с неограниченной емкостью и т.д., что, скажем, не имеет места в технологии производства современных интегральных схем, где размеры каждого элемента соизмеримы с линейными параметрами таких величин, как дебаевская длина экранирования и толщина поверхностной фазы. Более того, явления, происходящие в областях, соизмеримых с дебаевской длиной экранирования, играют существенную роль и в процессе эксплуатации приборов, созданных на современной элементной базе. Например, актуальнейшей проблемой нынешней электроники является деградация электрофизических параметров приборов. Представляется, что решение указанной задачи в первую очередь напрямую связано и с пониманием диффузионных процессов, обусловленных электрическими полями, возникающими из-за наличия контакта разнородных материалов. Поэтому далее в настоящей работе мы получим уравнение диффузии, учитывающее влияние внутренних электрических полей на процессы переноса вещества, и покажем некоторые его решения.

Ну а сейчас вернемся к началу настоящего параграфа и обсудим, за счет чего же была достигнута простота решения уравнения диффузии? Во-первых, за счет сильного упрощения (идеализации) самой задачи, а, во-вторых, не обошлось и без некоторого лукавства. Если еще раз более внимательно проанализировать предыдущий раздел, то легко можно обнаружить, что и описание задачи (26)-(28), и рис. 1, и решения, представленные на рис. 3, содержат некоторые недомолвки - "параметры по умолчанию", но при этом вроде бы, дополняя друг друга, создают образ единой задачи. А, по сути, в (26)-(28) просто отсутствуют граничные условия за исключением области окна в маске. Рисунок 1 несколько восполняет этот недостаток своей волнистой обрамляющей, намекающей на бесконечную протяженность объекта исследования, и словом "маска", говорящему всякому, кто имеет хотя бы какое-либо отношение к технологии полупроводниковых приборов, что поток диффузанта через маску, по крайней мере, значительно ограничен по сравнению с областью окна или даже равен нулю. Предположив, что концентрация диффузанта в области маски всегда равна нулю, мы, решая уравнение диффузии, вроде бы использовали последнее более сильное условие. И в результате решили … совершенно не ту задачу. Оказывается, что область маски в решенной задаче не только препятствует проникновению диффузанта в кристалл, но и одновременно является местом стока диффузанта неограниченной емкости. Действительно, обратившись к рисунку 3, видим, что изоконцентраты закреплены у краев маски и с течением времени упорно "не хотят" перемещаться вдоль границы "маска – кристалл" несмотря на то, что градиенты концентрации у краев маски в этих направлениях, а, следовательно, и движущие силы диффузии, максимальны. А ларчик тут открывается просто. Поговорить то мы поговорили об отсутствии потоков диффузанта через маску, но ничего для того, чтобы их не было, не предприняли.

В целом наши проблемы связаны с тем, что мы неверно решали задачу. Подход к методике решения настолько плох, что, как я надеюсь, читатель и сам уже пришел к выводу, так задачи решать нельзя. За многие тысячелетия человечество сформулировало подходы, методики, алгоритмы решения проблем. Они просты и очевидны и, наверное, именно поэтому мы зачастую ими не пользуемся, не придаем им должного значения и как следствие затрачиваем огромные усилия на то, что делать и не следовало бы, но, к сожалению, об этом узнаем, когда работа уже выполнена или не узнаем вовсе. Современные представления о том, как подходить к решению той или иной проблемы можно найти в прекрасной книге профессора МИСиС М.А. Штремеля [7] и в литературе, указанной в библиографическом списке к ней. Для тех, кто занимается моделированием можно порекомендовать [8]. И, тем не менее, очень кратко сформулируем основные этапы решения задач, аналогичных решенной нами:

1. Постановка задачи и определение конечных целей – заключается в содержательной (физической) постановке задачи: выбор общего подхода, определение совокупности критериев, которым должна удовлетворять система, задание условий ее реализации (работы), выбор системы координат и масштабов, выбор системы единиц.

2. Построение математической модели – математическая формулировка задачи. Модель должна правильно (адекватно) описывать основные законы физики, функционирования объекта (системы, устройства), влияние параметров на его свойства и т.д.

3. Разработка (выбор) численного метода. Математическая формулировка задачи может оказаться непереводимой непосредственно на язык ЭВМ. Для того, чтобы задача была правильно решена, необходимо уравнения, функции, интегралы выразить через элементарные функции. Более того, необходимо убедиться, что никакие погрешности, содержащиеся в исходных данных или внесенные в процессе вычислений, не повлияют сколько-нибудь заметным образом на точность результатов.

4. Разработка алгоритма – выражение решения задачи в виде точно определенной последовательности операций вычислительной машины. Раньше алгоритм представляли графически в виде блок-схемы. Сегодня от графического представления чаще всего отказываются (особенно в случае больших задач), но словесный портрет задачи в виде последовательности операций с указанием условий переходов от пункта к пункту записывают всегда.

5. Программирование – изложение алгоритма на языке "понятном" ЭВМ.

6. Отладка программы. Очень важный этап в процессе решения задачи. С этой целью используют тестовые задачи, пошаговое выполнение программы.

7. Выполнение расчетов. Если это возможно, то выполняют для нескольких наборов данных.

8. Анализ и оформление результатов.

Самое главное в этом процессе заключается в том, что на каждом этапе необходимо постоянно обращаться к предыдущим и, в особенности к первому, задавая себе вопрос "А ту ли задачу я решаю?"

Задача 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) при h₁ = h₂ = h:

= + ( + ).

При этом численное значение величины D ∙τ / h ², равное 0.25, сохранено.

· В строках 41 – 44 реализовано выполнение условия отсутствия потока на границе "маска – кристалл". Концентрация примеси во всех узлах в области маски приравнивается к концентрации примеси в области кристалла в узлах, расположенных непосредственно под первыми, то есть c (x_i, 0)= c (x_i, 1), где i – номера узлов в области маски вдоль оси X -ов.

Имея ввиду, как отмечалось выше, малость объектов современной электроники и как следствие необходимость учета геометрической формы этих объектов и явлений, сопровождающих диффузионные процессы в очень тонких слоях, рассмотрим задачу о диффузии, где будет использован тот факт, что в различных материалах коэффициент диффузии одной и той же примеси неодинаковый и на границах раздела фаз может быть в десятки раз больше [1], чем в объеме материала. Кроме того, в отличие от предыдущей задачи, заменим источник диффузанта бесконечной емкости на источник ограниченной емкости (например, диффузия из силикатного стекла).