Содержание к диссертации
Введение
1. Теоретические основы моделирования процессов самоорганизации в явлениях адсорбции и порообразовния 9
1.1 Клеточные автоматы 9
1.2 Автоволновые структуры 14
1.3 Классические модели физической адсорбции 15
1.4 Классические модели порообразования 18
1.5 Постановка задачи 28
1.5.1 Адсорбционные явления на поверхности 28
1.5.2 Процессы порообразования в глубине кристалла 31
2. Дискретная модель физической адсорбции с конечным числом состояний 39
2.1 Двумерная дискретная модель адсорбции с тремя состояниями 39
2.1.1 Формула отображения 42
2.1.2 Реализация алгоритма в вычислительной среде Matlab 49
2.2 Обсуждение результатов моделирования процессов адсорбции 56
3. Дискретная модель порообразования в детерминированных и стоха стических полях 69
3.1 Расчет распределения потенциального поля 69
3.1.1 Решение электростатических задач в случае зависимости потенциала от одного параметра (метод Ламе) 69
3.1.2 Расчет поля при логарифмической, обратно - пропорциональной аппроксимациях скорости раскрытия ветвей гиперболы 73
3.2 Моделирование процесса блуждания частицы 78
3.3 Моделирование кластерообразования 82
3.3.1 Модель образования кластеров в электростатическом поле 85
3.3.2 Модель роста с иссякающим источником 87
3.3.3 Модель роста с корректировкой траектории движения диффундирующей частицы 91
3.3.4 Модель роста кластера в полях неэлектрической природы 95
3.4 Обсуждение результатов моделирования процессов порообразования 100
Заключение 106
- Классические модели физической адсорбции
- Обсуждение результатов моделирования процессов адсорбции
- Расчет поля при логарифмической, обратно - пропорциональной аппроксимациях скорости раскрытия ветвей гиперболы
- Модель роста с корректировкой траектории движения диффундирующей частицы
Введение к работе
Коллективные явления самоорганизации изучаются в стремительно развивающейся междисциплинарной области исследований, закономерности которой обладают достаточно широкой общностью. Сфера применения основных закономерностей синергетики охватывает области многих наук. Исследование явлений самоорганизации в различных физических процессах ведется сравнительно давно. Экспериментально наблюдаемые явления самоорганизации чаще всего выражаются в синхронизации динамических параметров системы, образовании кластерных структур или появлении автоколебательных и автоволновых процессов. К настоящему времени построено большое число математических моделей, с разной степенью адекватности объясняющих подобные явления. Процессы адсорбции и кластерообразо-вания являются яркими представителями процессов самоорганизации. Например, физическая адсорбция при математическом моделировании до на стоящего времени не изучалась в качестве динамического процесса, а построенные математические модели не объясняют хаотического поведения параметров адсорбционной системы. При математическом моделировании процессов порообразования, рассматриваемых как пример кластерной системы, получены различные кластерные структуры, но в этих моделях не учитывается динамика поверхностных явлений. Представленная совокупность экспериментальных фактов из области физики взаимодействия по -5-верхности с активными газовыми и жидкими средами к настоящему времени не находят исчерпывающего объяснения в рамках единой теории.
Приведенные замечания свидетельствуют об актуальности темы данного исследования, и поэтому целью исследования является развитие с помощью методов математического моделирования теории, объясняющей совокупность экспериментальных данных в рассматриваемых областях физики.
В качестве объекта исследования при математическом моделировании кластерообразующей системы в диссертационной работе рассматриваются закономерности формирования пористого пространства. Процесс представляется как совокупность двух взаимосвязанных стадий - стартовыми процессами считаются динамические адсорбционные явления на поверхности, от них зависят явления образования кластеров, протекающие в глубине кристалла под поверхностью.
Обращаясь к проблеме исследования адсорбции различных веществ на поверхности твердых тел, можно сказать, что она привлекает к себе пристальное внимание исследователей на протяжении длительного времени, так как поверхность является единственным каналом проникновения вглубь кристалла [1]. Процессы на поверхности определяют также состояния более сложных систем, состоящих из граничащих друг с другом сред, например, газа или жидкости с твердым телом. Кроме того, модификация свойств поверхности существенно сказывается на изменении объемных свойств твердого тела, особенно в областях с малыми размерами кристаллитов, где объем и площадь поверхности становятся соизмеримы друг с другом. Изменение объемных свойств вещества непосредственно связано с процессами, протекающими на поверхности твердого тела, что определяет как динамику процессов, так и физические свойства конечного продукта.
В ходе процессов анодирования полупроводниковых кристаллов кремния в плавиковой кислоте в определенных режимах было обнаружено периодическое изменение динамических характеристик процесса с изменением времени [2-9]. Осцилляционные явления были обнаружены во время процессов дотравливания пористого кремния. Аналогичные явления были также выявлены в ходе экспериментов с электролюминесцентными устройствами [10] на основе пористого кремния, помещенными в среду поверхностно активных веществ, как правило, содержащих полярные молекулы [11,12]. В работе [13], были зарегистрированы периодические изменения динамических переменных при снятии вольт-амперных характеристик (ВАХ) структур на основе пористого кремния и помещенных в среду с полярными молекулами.
Адсорбция и десорбция молекул (атомов), способных изменять свое зарядовое состояние, также проявляет периодичность развития динамики процесса во времени. В работе [14,15] обнаружены автоколебательные эффекты при полевой десорбции калия с поверхности, на которой находится смесь адсорбатов калий - золото. Эти эффекты состояли в периодических изменениях ионного тока и адсорбционных изображений при постоянных условиях эксперимента [14,15].
Все описанные выше факты свидетельствуют о необходимости развития теории, способной объяснить совокупность экспериментальных данных, наблюдаемых в рассматриваемых областях физики.
Когерентные явления в стохастических системах представляют интерес с точки зрения применения их в различных физических областях. Явления когерентности обуславливают при определенных условиях формирование кластеров в стохастических системах. Кластеризация может быть описана в рамках многих моделей, которые по некоторым аспектам схожи с моделями, описывающими когерентные явления в динамических системах. Одной из задач данного класса является кластеризация примеси в случайных полях скоростей. Принципиальной особенностью этой задачи является существование кластерной структуры распределения концентрации примеси. Аналогичная особенность проявляется в виде возникновения каустической структуры в результате фокусировок и дефокусировок света, распространяющегося в случайной среде [16]. Подобная задача возникает, например, при распространении луча лазера в волноводе и связана с образованием, так называемых, "спеклов" [16]. Простейшим случаем, описывающим движение частицы в поле случайных скоростей, является система обыкновенных уравнений первого порядка [16]:
Еще одним примером явлений, связанных с формированием кластеров, является образование пористого пространства в полупроводниковых кристаллах вследствие различного рода процессов. Представителями веществ с подобными свойствами являются пористый кремний [17], формирующийся в ходе анодирования в растворах плавиковой кислоты, а также пористые вещества, образованные спеканием мелкодисперсных порошков [18].
Компьютерное моделирование позволяет решить достаточно широкий круг алгоритмически сформулированных задач, исследование асимптотического поведения решений которых иногда не представляется возможным иными способами. Эти задачи связаны прежде всего с ситуацией, когда достаточно четко описываются лишь локальные свойства системы, а учет глобальных свойств достигается посредством усреднений. Попытки использования компьютерных моделей для объяснения закономерностей порообразования предпринимались неоднократно [19]. Однако, исчерпывающей модели, описывающей формирование пористых пространств, к настоящему времени не построено, что свидетельствует о необходимости продолжения исследований в данном направлении.
Таким образом, исследуемая проблема находит свое применение в широком спектре прикладных задач современной физики и технологии.
Классические модели физической адсорбции
Явление адсорбции известно очень давно. Такие природные материалы, как песок и почва, использовали для очистки воды еще на заре человеческого общества. В конце XVIII века К. Шееле и одновременно Фонтана обнаружили способность свежепрокаленного угля поглощать различные газы в объемах, в несколько раз превышающих его собственный объем. Вскоре выяснилось, что величина поглощенного газа зависит от типа угля и природы газа. Адсорбция (от лат. ad - на, при и sorbeo - поглощаю) - это преимущественное концентрированное молекул газа или растворенного в жидкости вещества (адсорбата) на поверхности жидкости или твердого тела (адсорбента), а также растворенного в жидкости вещества на границе ее раздела с газовой фазой [34]. Количественной характеристикой адсорбции является величина Г, представляющая собой избыток адсорбата, приходящийся на единицу площади поверхностного слоя, по сравнению с количеством адсорбата в единицу объема фазы адсорбента. Отношение 0 = /р называют степенью (или долей) покрытия поверхности (Гоо - предельно возможная величина монослойной адсорбции для данной системы).
Сегодня адсорбция составляет основу многих промышленных операций и научных исследований. Наиболее важные из них - очистка, выделение и разделение различных веществ. Исследования поверхности тесно связаны с развитием полупроводниковой техники, медицины, строительства, военного дела - одним словом везде, где существенную роль играют поверхностные явления [35].
Следует отметить, что не существует особых сил, вызывающих адсорбцию. Явление адсорбции связано с тем, что силы межмолекулярного взаимодействия на границе раздела фаз не скомпенсированы, и, следовательно, пограничный слой обладает избытком энергии - свободной поверхностной энергией. В результате притяжения поверхностью раздела фаз находящихся вблизи нее молекул адсорбата свободная поверхностная энергия уменьшается, то есть процессы адсорбции энергетически выгодны [1].
В зависимости от характера взаимодействия молекул адсорбата и адсорбента различают физическую адсорбцию и хемосорбцию. Физическая адсорбция обусловлена силами межмолекулярного взаимодействия и не сопровождается существенным изменением электронной структуры молекул адсорбата. Физическая адсорбция может быть как монослойной (с образованием мономолекулярного слоя), так и полимолекулярной (многослойной). При физической адсорбции адсорбированные молекулы обычно обладают поверхностной подвижностью.
При хемосорбции между атомами (молекулами) адсорбата и адсорбента образуется химическая связь, таким образом хемосорбцию можно рассматривать как химическую реакцию, область протекания которой ограничена поверхностным слоем. В некоторых случаях на одной поверхности могут протекать оба типа адсорбции одновременно [34].
Обратный адсорбции процесс, при котором адсорбированные частицы покидают поверхность адсорбента, называется десорбцией. Десорбция происходит в результате колебательного движения адсорбированной молекулы вдоль направления действия силы притяжения между адсорбатом и адсорбентом [34].
Для данного адсорбата и адсорбента равновесная величина 9 для газа или пара определяется температурой Т и давлением р. При исследовании адсорбции одну из этих величин обычно поддерживают постоянной. Чаще всего в результате непосредственных измерений получают зависимость в от р при постоянной температуре, которая носит название изотермы адсорбции. В принципе в изотермах адсорбции содержится вся информации об адсорбционной системе. Так, измеряя изотермы адсорбции при разных температурах можно определить зависимость р от Т для данной величины адсорбции. Такие зависимости называются изостерами адсорбции. Из уравнения изостеры где R - газовая постоянная, а А - константа, вытекает способ вычисления теплоты адсорбции q из изотерм для данной величины адсорбции. Теплота адсорбции является основной энергетической характеристикой, определяющей природу и величину сил, вызывающих адсорбцию. Исследуя зависимость теплоты адсорбции от количества адсорбированного вещества можно сделать определённые выводы о характере поверхности адсорбента и величине взаимодействия адсорбат - адсорбат [35].
Единой теории, которая описывала бы любые процессы адсорбции, пока не создана; существующие частные теоретические разработки основываются на различных моделях [1,34]. Модель локализованной (или центровой) адсорбции предполагает наличие на поверхности адсорбента так называемых центров адсорбции, представляющих собой либо строго определённые участки поверхности, на которых образуется сильная адсорбционная связь, либо распределенные по поверхности двумерные ячейки со слабым адсорбционным полем (полем сил межмолекулярного взаимодействия). К этому классу моделей относится знаменитая полимолекулярная теория Брунауэ-ра - Эмметта - Теллера (БЭТ) [34,36], в которой не учитывается латеральное взаимодействие. В последнем случае предполагается наличие плотной упаковки молекул адсорбата на поверхности в пределах рассматриваемой ячейки. В основе модели двумерной фазы лежит положение о том, что адсорбированный монослой представляет собой не идеальный двумерный газ, однако полимолекулярное покрытие в данной модели не рассматривается. И, наконец, потенциальная модель адсорбции базируется на представлениях о потенциальном поле поверхности твердого тела, в котором адсорбированный газ сжат вблизи поверхности и разрежен в наружных слоях [34].
Обсуждение результатов моделирования процессов адсорбции
Отметим, что общее отображение в виде (2.5) - (2.7) дает настолько огромный спектр динамических структур различных конфигураций и цветовых сочетаний, что исчерпать на данный момент все возможные исходы отображения (2.5) - (2.7) не представляется возможным. Остановимся, однако, на некоторых закономерностях, обнаруженных нами на данный момент. Ряд характерных особенностей нашей модели рассмотрен в работе [44]. В этой работе [44] приведены характерные динамические структуры, образующиеся в режимах формирования ведущих центров, а также зависимость черных клеток от времени (тактов), отражающая хаотическую, периодическую динамику системы, и динамику в режиме ведущих центров для некоторых значений параметров системы.
Приведем несколько рисунков динамических структур, формирующихся в результате отображения (2.5) - (2.7), а также зависимости числа черных (В), серых (G) и белых (W) клеток от времени (тактов). На рисунках А.1 -А. 12 приведены снимки спиральных волн, состоящих как из серых на фоне белых (рис. А.З, А.4), из черно-серых на фоне белых (рис. А.1 - А.2), из черных на фоне серых (рис. А.5), или из серых на фоне черных (рис. А.б - А.8, АЛО) и тоже, но с участием белых (рис. А.9 (а)), из серо-черно-серых на белом фоне (рис. А.11) и соответствующие зависимости от времени (тактов) числа клеток (ячеек) в соответствующих состояниях. Спиральные волны состоят из двух участков с различной киральностью, состав (наполненность клетками разных цветов) же их может быть различным, как это видно из рисунков.
Часто солитоны обладают большой скважностью, то есть расстояние (во времени) между соседними волнами оказывается во много раз большим, чем длительность самого колебания. В таких случаях энергия колебаний сосредотачивается на гораздо более высоких частотах, чем частота следования волн (рис. А.2, А.З, А.5).
Наличие или отсутствие члена Pm (t),mtJ(t+i) позволяет нам исследовать поведение системы с учетом стохастического характера переходов или без него. Если алгоритм носит детерминированный характер, то есть член Р не учитывается, то для начала развития системы мы случайным образом распределяем по решетке несколько серых клеток. При этом, изменяя локальные правила, мы получаем решения, характер поведения которых зависит или не зависит от распределения стартовых клеток. Поскольку количество соседей, чье состояние учитывается в алгоритме, равно 8 (область Мура), некоторые примеры автоволн демонстрируют восьмиугольную или прямоугольную симметрию.
Отметим, что структуры с резко обозначенными углами зарождаются при отсутствующем вероятностном факторе Pmij{t),midt+i)i в т0 время как при учете его влияния могут быть получены структуры с круговой симметрией, то есть симметрией более высокого порядка (см. рис. 2.6).
На рисунках (А.14 - А.24) приведены мгновенные снимки и динамические зависимости числа различных точек от времени (тактов) для систем, находящихся в различных хаотических режимах. Динамическое поведение в этом случае характеризуется большим разнообразием.
Подчеркнем, что предметом исследования, рассматриваемым в настоящей главе, являются динамические свойства адсорбционных систем с учетом возможности изменения состояний основных, базовых, структурных элементов (атомов, молекул), тогда как традиционный подход концентрирует внимание на изучении, главным образом, некоторых усредненных глобальных характеристик системы (см. например, [1,51,52]).
Разработанная модель объясняет достаточно широкий круг физических и физико-химических явлений до настоящего времени представлявших собой лишь совокупность разрозненных фактов из различных областей. Рис. 2.7. Полевые десорбционные изображения непрерывного режима при разрушении пленки соединения КАи электрическим полем [15]. 1. Становиться понятным качественное поведение адсорбционной системы КАи па поверхности вольфрамовой иглы [14,15]. На рис. 2.7 изображено временное поведение этой системы, полученное в работе [14,15]. Наша модель описывает подобное поведение системы, когда имеют место расходящиеся из некоторого центра структуры (см. рис. 2.8).
Расчет поля при логарифмической, обратно - пропорциональной аппроксимациях скорости раскрытия ветвей гиперболы
Считаем, что параметр b, входящий в формулу (3.2), не зависит от Л. Для реализации алгоритма используются два вложенных цикла, в теле которых последовательно изменяются значения переменных А& и хг. Зафиксировав эти числа, вычисляется соответствующее значение у по формуле (3.2). Далее находим такое j, что уэ-\ у у3+\- В итоге получим пространственную матрицу (хі,уу) и соответствующие этим координатам значения потенциала поля (р 3 и значения А&. Индексы i,j, к соответствуют геометрии половины поля одной затравки.
Затем, используя принцип суперпозиции полей и симметрию задачи, формируется суммарное поле трех эквидистантных затравок, то есть полученную матрицу, содержащую значения потенциала поля ср обрезаем до индекса, соответствующего расстоянию dmax. Зеркально отобразим относительно вертикальной центральной оси полученный результат, и после конкатенации с матрицей ip y получим матрицу с распределением потенциала поля крайней левой затравки. Матрица, содержащая распределение для правой крайней затравки, получается зеркальным отображением полученной матрицы. Если аналогичную процедуру проделать, обрезая матрицу p j до индекса соответствующего расстоянию D-\-dmax, получим распределение потенциала поля центральной затравки. Суммируя почленно полученные три матрицы, в итоге получаем дискретное распределение потенциала поля трех затравок (ptj. Проделав аналогичную процедуру, можно получить матрицу значений параметра А,у. Здесь каждой клетке (хиУ3) ставится в соответствие значение рассчитанного поля и параметра А. Причем, индексы г, j здесь соответствуют геометрии всей рабочей области. Описанным образом можно рассчитать электростатическое поле необходимого количества затравок. В данной работе используется поле трех затравок, чтобы иметь возможность анализировать влияние полей крайних затравок на поле цен-тральной.
Разность потенциалов между иглой и нижней границей равна cpr = ZSiVo, Si - диэлектрическая проницаемость кремния, щ - параметр размерности разности потенциалов, используемый при моделировании.
Электростатическое поле трех затравок, полученное описанным выше способом, показано на рисунке (3.3).
Для моделирования медленно спадающих на больших расстояниях потен-циалов параметр (б2 (А) — А2) аппроксимировался логарифмической и обратно-пропорциональной зависимостью, то есть полагалось (б2 (А) — А2) = / (А) при / (А) = р In (A) 4- q, f (A) = р/Х + q. При этом "разворот" ветвей гипербол происходит более медленно с расстоянием, то есть на более далеких расстояниях от верхней границы.
При рассмотрении двумерного случая с одной игольчатой затравкой эквипотенциальные кривые описываются уравнением гиперболы с параметрами А и (б2 —А2)1 2 [58] (см. рисунок 3.4). В предыдущем разделе (стр. 70), параметр Ь в формуле (3.2) мы считали не зависящим от параметра А. Если величину (2 — А2)1 2 считать линейно зависимой от А, то можно положить, что
Коэффициенты p,q вычисляются исходя из следующих соображений. При стремлении А к своему минимальному значению Amjn ветви гиперболы разворачиваются в прямую, то есть (б2 — А2) — со. Поэтому можно ввести достаточно большое число dmax такое, что [b2 (Лщш) ] = тах-Второе условие состоит в том, что при А — Атах ветви гиперболы прижимаются к игле, то есть \Ь2 (Amax) — A J = (twin (см. рис. 3.4). Решая систему из двух уравнений, находим параметры p,q (3.8, 3.9).
При моделировании потенциальных полей также использовались другие способы задания скорости раскрытия ветвей гиперболы, отличные от линейной аппроксимации: логарифмическая аппроксимация и обратно-пропорциональная
Необходимость таких полей будет объяснена ниже (п. 3.3.4). Коэффициенты р , q , р", q" вычисляются аналогично, как и в случае линейной аппроксимации (см. (3.10)), подстановкой логарифмической или обратно-пропорциональной зависимостей.
Для учета описанных выше правил изменения Ъ (А) (3.11, 3.14) в алгоритме, показанном на рисунке 3.2, нужно внести дополнительную операцию: в вычисляемую зависимость у (3.2) добавить расчет Ъ (А) в зависимости от способа аппроксимации скорости раскрытия ветвей гиперболы (3.11, 3.14).
Формирующиеся кластерные структуры становятся более вытянутыми. Расчет потенциала с учетом аппроксимации скорости раскрытия ветвей гипербол логарифмической и обратно пропорциональной зависимостями представлен на рисунке 3.5 (а, Ь).
Модель роста с корректировкой траектории движения диффундирующей частицы
Детальное исследование построения результирующего потенциала и вычисленных вероятностей перехода в различные состояния обнаруживает следующий эффект. На вертикальной линии строго под любой из затравок вероятность перехода в сторону выше вероятности переходов в вертикальном направлении. Эти вертикальные линии не являются линиями большего притяжения, чем ожидалось из качественных рассмотрений, в виду максимального значения электрического поля вдоль них, а имеют характер линий локального рассеяния, вероятность ухода с которых выше вероятности остаться на этой линии. Этот эффект связан с тем, что на вертикальной линии потенциал достаточно высок, но является слабо меняющейся величиной от клетки к клетке вдоль этой линии, тогда как при переходе с этой линии в боковую клетку потенциал меняется более значительно. Этот эффект приводит к тому, что несмотря на потенциалы, меняющиеся в достаточно узкой области, происходит уход блуждающей частицы с вертикали, в результате чего формируются более разветвленные кластеры, чем при подавлении этого эффекта.
С целью подавления отмеченного выше эффекта и анализа условий, необходимых для формирования неветвящихся, изолированных, вертикальных кластеров была проделана следующая процедура. Для объяснения ее принципов работы определим функцию распределения вероятности (ФРВ). / = /о + &! гДе /о невозмущенная стационарная изотропная и пространственно однородная функция распределения, Sf - добавка, описы Эта процедура применяется в двух вариантах. Первый - вероятность прыжка корректируется независимо от положения частицы. Второй - модуляция вероятности производится только напротив затравок. Если Z% -абсциссы затравок, где і = 1 : 3, то вероятность прыжка корректируется в области Zt ± Az.
Необходимо отметить, что средняя длина траектории и средняя скорость смещения в вертикальном направлении при данных условиях больше по сравнению с блужданием частицы по основному алгоритму (см. п. 3.3.1) и алгоритму с иссякающим источником (см. п. 3.3.2). Рисунок 3.16 демонстрирует зависимость среднего отклонения частицы от времени роста кластера.
На рисунке 3.17 изображены примеры кластерных структур, полученных при первой и второй стадиях возмущения функции распределения вероятности. Pk = 0.1. Как видно из рисунков при увеличении параметра pk расстояние между ветвями кластера значительно уменьшается. Это можно объяснить уменьшением времени блуждания частицы в непосредственной близости от кластера. В следствие чего, блуждание частиц становится больше похоже на их падение на верхнюю границу образца. і Структуры В.2(a) и В.2(Ь,с) получены при различных параметрах при s соединения частицы к кластеру, rj = \/2Дж/8 и 7) = 2Лж/8, соответственно.
На рисунках (В.З, В.4) показаны кластерные структуры, полученные при первом способе корректировки вероятности прыжка частицы. Блуждание происходит, используя основной алгоритм 3.3.1. Рисунки отличаются шириной коридора Az и величиной j &.
На рисунках (В.5 - В.8) представлены структуры, полученные, используя алгоритм иссякающего источника (3.3.2).
Применение аппроксимаций, обуславливающих различные способы раскрытия ветвей гипербол и тем самым моделирующих поля неэлектрической природы, связано со следующим фактом. Задача по точному расчету деформационного потенциала и его градиента от трех описанных выше затравок достаточно сложна. Решение задачи методом конформных отображений в случае полей, отличных от электростатических, становится невозможным. Например, распределение полей, полученное для простейших конфигураций электродов, можно было бы с помощью соответствующего конформного преобразования отобразить в распределение полей при более сложной конфигурации электродов. В случае сил Казимира (F г-4) подобная ситуация разрешается достаточно нетривиальным способом, например, в случае взаимодействия между двумя плоскостями применяется так называемая "теорема близости" ("proximity force theorem"), что позволяет рассчитать взаимодействие между сферической поверхностью и плоскостью на близких расстояниях [65], исходя из известного решения для взаимодействующих плоскостей.
Как показывают исследования компьютерных моделей, для получения изолированных, неветвящихся, вертикальных кластеров необходим потенциал, имеющий больший градиент, чем тот, который соответствует распределению электрического напряжения в пластине [67].
