Отображение физических процессов на их Клеточно-автоматные модели Текст научной статьи по специальности «Математика»

2008

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Управление, вычислительная техника и информатика

№ 2(3)

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 539.2: 536.1

О.Л. Бандман ОТОБРАЖЕНИЕ ФИЗИЧЕСКИХ ПРОЦЕССОВ НА ИХ КЛЕТОЧНО-АВТОМАТНЫЕ МОДЕЛИ

Задача сопоставления физических величин модельным значениям (масштабирование) возникает при клеточно-автоматном (КА) моделировании дважды: при построении КА-модели по заданному описанию процесса и при интерпретации результатов моделирования. В статье предлагается систематический подход к решению этой задачи: формулируются общие принципы построения КА-моделей, на основе которых детально разработаны методы масштабирования для трех основных типов КА-моделей (диффузионных КА, КА-гидродинамики, асинхронных кинетических КА).

Ключевые слова: клеточно-автоматное моделирование, пространственная динамика, реакционно-диффузионные процессы, клеточно-автоматная гидродинамика, поверхностная химия.

Клеточный автомат (КА) в последние годы привлекает большое внимание как математическая модель пространственной динамики. Поскольку КА способны моделировать нелинейные и разрывные процессы, их иногда считают дополнением к дифференциальным уравнениям с частными производными [1, 2]. В настоящее время уже известно довольно много разных КА, эволюция которых моделирует пространственную динамику в физике, биологии, химии, а также поведение колоний животных и толпы людей [3]. Наиболее известные КА-модели можно подразделить на три типа: 1) КА-модели процессов, в которых диффузия присутствует в явном виде [4], 2) КА-модели, получившие названия «решеточный газ» (неудачный перевод слов «Lattice-Gas models»), которые далее называются КА-гидродинамикой [5 - 7] и 3) асинхронные вероятностные КА-модели кинетических взаимодействий наночастиц (называемые химиками методами Монте-Карло) [8,9]. Хотя эти типы КА достаточно хорошо изучены, высокий уровень абстракции клеточного автомата вызывает иногда большие трудности в сопоставлении КА-модельных величин соответствующим реальным физическим, а именно размеров КА, его начального состояния, правил переходов, вероятностей их применения. Обратная проблема возникает при интерпретации результатов моделирования.

В статье делается попытка сформулировать общие принципы сопоставления физических величин модельным и показать их применение на трех типовых примерах КА-моделирования.

1. Общий подход к построению КА-моделей

Задача моделирования какого-либо естественного или лабораторного процесса ставится обычно следующим образом. Известны следующие величины: 1) разме-

ры области, в которой происходит моделируемый процесс, 2) свойства среды, в которой он происходит, 3) начальные и краевые условия, 4) внешние воздействия.

Главный результат моделирования - это скалярная u(x,t) или векторная u(x,t) функция, представляющая пространственную динамику изучаемого процесса (распределение концентраций веществ, участвующих в процессе, интенсивность магнитного или электрического поля, поле скоростей потока и т.д.).

Чтобы достигнуть цели моделирования, первое, что надо сделать, - это правильно построить КА-модель, для чего необходимо соблюдать следующие принципы.

1) КА-модель должна быть справедливой для всей области значений моделируемых величин. Например, при моделирования потока жидкости с помощью КА необходимо убедиться, что коэффициент Рейнольдса потока не превысит ограничений КА-гидродинамики [5, 7].

2) Размеры КА-модели должны быть такими, чтобы она была реализуема на имеющемся оборудовании в приемлемое время.

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

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

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

2) Опыт применения КА-моделирования еще не достаточно богат, чтобы считать некоторые масштабы достаточно достоверными и проверенными на практических задачах или в натурных экспериментах.

3) КА-модели разнообразны, большинство из них существенно нелинейны с совершенно непредсказуемым поведением, присущим так называемым сложным системам [2].

Физические величины далее будут выражаться в системе МКС (метр, килограмм массы, секунда). Модельные величины безразмерны: они имеют смысл количества клеток с длиной стороны l = 1, числа частиц N с массой m = 1, вектора скорости с модулем |v| = 1. Масштабные коэффициенты (масштабы) определяются как отношения физической величины Z к соответствующей модельной z и обозначаются как pz = Z/z, где z - любая величина, участвующая в процессе. Масштабы в КА-моделировании разделены на три группы.

Первая группа включает в себя основные масштабы, которые применяются во всех методах численного моделирования. Это масштаб времени р, = т и масштаб длины р/ = h. В КА-моделях р, соответствует времени смены глобального состояния КА в секундах, а р/ - линейному размеру клетки в метрах.

Вторая группа включает в себя масштабы свойств веществ, которые участвуют в моделируемом процессе (вязкость, плотность, проводимость и др.). Эти свойства являются характеристиками (параметрами) модели. Они обычно выводятся из

правил перехода КА или получаются путем выполнения вычислительных экспериментов над заранее известными процессами.

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

Для выражения коэффициентов масштабирования и иллюстрации их применения используется формализм алгоритма параллельных подстановок [10]. В соответствии с ним КА-модель К представляется четырьмя символами К = {А, М, ©, р), где А — алфавит, на который не налагается никаких ограничений, М ={ш\, т2,..., тк...} — счетное множество имен клеток или точек в пространстве, ® ={9Ь...,9Ч} — множество локальных операторов, р — режим функционирования, определяющий распределение применений локальных операторов по пространству и времени. Основным понятием в КА-модели является клетка, представленная парой (а,т), а е А , т е М. Множество клеток О = {(ак,тк)}, не содержащее клеток с одинаковыми именами, называется клеточным массивом.

На множестве М определены именующие функции ф(т) е М, значения которых определяют местоположения клеток, взаимодействующих с клеткой т. Когда в качестве имен используются декартовы координаты М={(г, у')}, именующие функции обычно даются в виде сдвигов фя(т) = {г+а,у+Ь}, где а, Ь — целые числа. Множество именующих функций

Т(т) = {т, ф1(т),..., ф„(т)}, п << |М|, (1)

называется шаблоном. Подмножество клеток Б(т) е О

Б(т)= {(и0,т), (иь ф!(т),...,(и„, ф„(т))}, (2)

с именами из Т(т) называется локальной конфигурацией над шаблоном Т(т). Локальный оператор 9(т )е © выражается в виде подстановки

9(т): Б(т) * Б"(т) ^ Б'(т) Ут е М, (3)

в которой шаблоны Т(т ) и Т(т) локальных конфигураций Б(т) и Б'(т) одинаковы, а Б'(т) может быть любой.

Применение 9(т) к определенной клетке (и,т) е О(?) состоит в следующем: в клеточном массиве О(?) множество клеток Б(т) заменяется на Б'(т). Такое определение локального оператора позволяет моделировать эволюцию живых организмов, которые растут и отмирают. При моделировании явлений физики применяются обычно стационарные локальные операторы [10], которые не меняют множества имен клеток, производя замену состояний клеток из Б(т) на состояния из 5'(т)={(и'0,т),(и'1,ф1(т), ..., (и'п, фп(т) )}, вычисляя их в соответствии с функциями переходов

и к fk (u1, ■ ., un, ип+1, ■ ■ ■ , un+h), (4)

где последние к переменных — состояния клеток из Б"(т). Эти состояния не изменяются оператором 9(т), а только служат дополнительными переменными, играя роль контекста. Применение 9(т) ко всем клеткам массива О(?) переводит массив в следующее глобальное состояние О(?+1) и называется итеративным шагом, или итерацией.

Существуют несколько разных режимов перехода О(?) ^ О(?+1). Синхронный режим требует, чтобы все аргументы функций переходов (4) были состояниями клеток из О(?), т.е. Б(т) с О (?) для всех т е М. Замена состояний во всех клетках происходит одновременно. Асинхронный режим предполагает, что каждая 9(т)

применяется по очереди к случайно выбранным клеткам, сразу же изменяя состояния в S(m). Таким образом, аргументами функции переходов могут быть состояния клеток как из ^(t), так и из Q(t+1). В обоих случаях итерацией называется применение ®(m) ко всем клеткам массива.

Последовательность

называется эволюцией КА.

КА моделирует какое-либо явление, если его эволюция соответствует (аппроксимирует) функцию и(ж,і), описывающую это явление. Соответствие это устанавливается только путем масштабирования пространственных координат (ж), времени (і) и моделируемых величин (и).

Диффузионными КА здесь называются такие КА-модели, в которых основу моделируемого процесса составляет диффузия. В дифференциальных уравнениях с частными производными диффузионный член представлен лапласианом вида С(иХХ+иуу), где С - диффузионный коэффициент (м2с-1). После дискретизации по явной схеме «крест» диффузионный член приобретает следующий вид:

где к = щ и т = р - масштабы длины и времени соответственно. В численных методах выбор к и т ограничен условием устойчивости и производительностью. В КА-моделировании р/ и р определяются исходя из модельного значения коэффициента диффузии СКа, который является характеристикой модели и для каждого КА, моделирующего диффузию, вычисляется аналитически из правил переходов или определяется экспериментально путем сравнения полученной КА-эволюции с известным моделируемым процессом. Поскольку С из (6) и Ска имеют один и тот же смысл, то Ска можно выразить аналогично (6) через масштабы р/ и р следующим образом:

Значение С и область моделирования /Х х /у х \ъ определены в задании. Поэтому достаточно выбрать модельные размеры КА х / х 4 , чтобы определить масштаб длины. Обычно размеры КА выбирают на основе компромисса между требуемым разрешением пространственной функции и имеющимися вычислительными ресурсами. Далее, когда масштаб длины определен, все остальные масштабы могут быть вычислены исходя из их физических отношений , т.е.

Существует несколько КА-моделей диффузии. Две из них наиболее известны. Первая, называемая блочно-синхронной моделью (КА-синх), предложена в [1] и исследована в [4]. Вторая, тоже предложенная в [1], называется наивной диффузией (КА-наив) представляет собой асинхронный вероятностный КА. Для всех диффузионных КА модельный коэффициент диффузии СКА является функцией от вероятности применения локального оператора, причем его максимальное значение Стах для каждой модели является ее инвариантом.

0(0), O(l), ..., □(/), 0(t+l ), ..., Q(4nd)

(5)

2. Масштабирование диффузионных КА

C (м;_1, j + м;, j+1 + мг+1, j + Mj, j+1 - 4м;, j), i = x/h, j = y/h , C = dT/h2, (6)

()

(8)

Таблица 1

Максимальные значения диффузии для блочно-синхронной и наивной КА-моделей

Размерность Стах (КА-синх) Стах (КА-наив)

Ш 1,0 1,0

2D 3/2 1/2

3D 23/18 1/3

Пример 1. Применение диффузионно-конвективного КА для моделирования движения насыщенного пара через пористую мембрану. Цель моделирования состоит в том, чтобы определить влияние свойств стенок пор на скорость движения пара в них. Исследуется фрагмент мембраны в двумерном приближении, диаметр поры принимается равным ее ширине. Фрагмент имеет три вертикально расположенные поры равного диаметра: с гидрофобными, гидрофильными и нейтральными стенками (рис. 1). Процесс движения пара представлен как поток абстрактных частиц, движущихся под действием конвективных и диффузионных сил. Источник частиц находится в верхней части фрагмента, так что конвективная сила заставляет частицы двигаться вниз. Диффузионная составляющая учитывает анизотропные свойства, вносимые стенками пор путем изменения вероятности перемещения частиц вдоль их диаметра. Отношение конвективной составляющей к диффузионной задано параметром Ре =0,5 [11]. Исследуемый образец имеет размеры 4 = 4 = I = 0,03 м, диаметр поры = 8• 10-3 м, коэффициент диффузии С = 10-4 м2с-1, плотность пара р = 0,5• 10-3 кг-м-3.

Моделирование производится с использованием асинхронной наивной диффузии. Алфавит А = {0,1}, М = {(г,у: г,у = 0,1,...Д-1)}. Локальный оператор работает как суперпозиция двух вероятностных подстановок © = бсш^б^). Вероятности их применения находятся в соответствии с заданным значением Ре = ^С0ОТ / рш = 0,5 [11]. Локальный оператор

0СОПУ (»'. 3) = {(1. (»'. 3)), (0,0', у))} Рсош > = {(0,0', у)), (1,0', у))}, (9)

действуя с с вероятностью ^еопу , перемещает частицу вдоль поры. Локальный оператор

всшг (Л 3) = {(ь 3)), ( 0 + а, у' + Ъ)} Рш > = {(0, (г, у')), (( + а, у' + Ъ))} , (10)

действуя с вероятностью р'ш, обменивается состояниями с одним из своих соседей, который, в свою очередь, выбираются с вероятностью 1/д, где д - число пустых соседей клетки (г, у). Таким образом, р'ш = рш • _р„ , где _р„ - величина, зависящая от типа поры и расстояния от клетки (г, у) до стенки. Если пора нейтральная или > в (в - расстояние от стенки, на котором ее влияние незаметно), то = 1. Если gw < в, то для гидрофильной и гидрофобной поры

р^ (рЫ1) = 1 - ехр(-#„ МвО, JPw (рЬоЪ) = - ехр(—/п1 в 1) соответственно, п1 , п2 - характеристики поверхности пор. Этих данных достаточно, чтобы получить все масштабы величин КА-модели.

1) Масштаб длины. В соответствии с практикой моделирования пористых срез [11], минимальный диаметр поры должен быть не менее десяти модельных единиц длины. Приняв р/ = 10-4 м, получим ДКА(поры ) = 80, а линейный размер фрагмента N = 300.

2) Масштаб коэффициента диффузии. Для используемой наивной диффузии СКА= 0,5 (табл. 1), следовательно, ^ = С/СКА = 2-10-4 м2с-1.

3) Масштаб времени, р = р//р^ = 0,5- 10-4с.

4) Масштаб плотности (масса абстрактной частицы) рр = р/(р0-р3) = 1,25-10-9 кг.

5) Масштаб потока (масса пара, проходящего через сечение поры за секунду, относительно числа частиц, выходящих из поры за итерацию). С учетом того, что площадь поры при ^КА(Поры ) = 80 равна £Р =5024 клеток, рд = р р / р4 = 2-10-5 кгс-1.

t = 0 t = 40 t = 100

Рис. 1. Три глобальных состояния фрагмента пористой мембраны при моделировании прохождения пара через гидрофобную, гидрофильную и нейтральную поры

Результатом моделирования являются следующие значения потоков, полученные в виде количеств частиц, выходящих из поры за итерацию:

0ка(Фо6) = 3939, 0ка(Фил) = 1975, ^кА(нейтр) = 2812,

которые пересчитываются в физические значения путем умножения на p,Q, что в результате дает

!0(фоб) = 0,079 кгс 1, |0(фил) = 0.039 кгс 1, ^кА(нейтр) = 0,058 кгс 1.

3. Масштабирование гидродинамических КА-моделей

КА-гидродинамика (Lattice-Gas Models) [5 - 7] объединяет несколько моделей процессов в газах и жидкостях, которые представлены в виде абстрактных частиц, движущихся в дискретном пространстве и сталкивающихся друг с другом. Наиболее изученной является серия вероятностных моделей, называемых FHP-моделями по инициалам их авторов [7]. В этой серии модель FHP-1 имеет двумерную гексогональную структуру клеточного пространства с расстояниями между центрами гексагонов hKA = 1. У каждой клетки есть шесть соседей, составляющих ее контекст T”(m) = (mb..., m6}. Алфавит состояний A состоит их булевых векторов s = (s0, sb..., s6) }, |A| = 27 . Компонент вектора состояния s, =1 (i=1,...,6) означает, что в клетке (s, mk) находится частица, движущаяся со скоростью с, = 1 по направлению к i-му соседу, s0 =1 означает, что в клетке находится частица покоя у которой с0 = 0. Частицы имеют модельную массу, равную 1. Режим функционирования КА синхронный двухэтапный, т.е. каждая итерация задается глобальной суперпозицией двух локальных операторов: 0 = 01(02). На первом этапе оператор движения заставляет каждую частицу, соответствующую s, =1, продвинуться на единицу длины по направлению приданной ей скорости с,:

01(m): ((s, ,m)} * ((si, ф,(m))} ^ ((s/,m )}, i = 1,2,...,6. (11)

На втором этапе оператор 02(т) моделирует столкновение частиц :

02(т): {(л ,т)} ^ {(л',т )}, (12)

причем если л имеет два единичных компонента, то функция переходов л' = / (л) имеет два равновероятных значения (рис. 2). Обычно, / (л) задается в форме таблицы.

Рис. 2. Графическое представление оператора столкновения для модели БНР-1

Трехмерная модель жидкого потока [12], названная автором моделью RD-1, имеет додэкаэдрическую структуру клеточного пространства. Контекстный шаблон состоит из 12 соседей. Хотя такая структура модель удовлетворяет условиям изотропии приближенно, ее экспериментальная проверка показала приемлемую точность моделирования [12] .

Гидродинамические КА-модели характеризуются следующими параметрами: равновесной плотностью рКА, скоростью звука иКА и вязкостью УКА(рКА). Модельная вязкость уКА зависит от средней плотности частиц

Рка (т) = , * I I5;(т), (13)

| Лу(ш) | |Ау(т)| і

где Лу(т) - область осреднения, содержащая подмножество клеток удаленных от т не более чем на Я (Я - радиус осреднения).

Таблица 2

Параметры основных КА-моделей потока жидкости

КА-модель иКЛ рКЛ ^КЛ ка)1

БНР-1 1/72 1,2 0,85

ЯБ-1 >/7/13 3,9 0,89

Постановка задачи моделирования включает в себя обычно следующие величины: 1) линейные размеры резервуара (м); 2) расположение в нем источников и стоков; 3) свойства жидкости: плотность р (кгм-3), кинематическая вязкость

1 Величины модельной вязкости даны с учетом поправки на дискретность, которая в [5,7] обозначается как g(p).

V (м2 с :); 4) давление (кг-см 2) или скорость потока и (м-с '), или коэффициент Рейнольдса

Яе = — . (14)

V

Сопоставляя характеристики моделей с физическими величинами, можно получить сразу масштабы вязкости и плотности

^ =-^> м2с-1; Цр = —^,кг• м3. (14)

УКЛ Ркл

Поскольку число Рейнольдса, являясь инвариантом, должно быть равным модельному Яе-кА, то из соотношения

Яе = ЯеКА = иКА — (15)

VКА

можно легко найти характеристические значения длины /КА или скорости иКА в зависимости от цели моделирования и затем подсчитать следующие масштабы:

1 —1 М/ /1

И/ =-—, м; Ми =—, м • с ; м =—, с- (16)

‘КА И/ М и

Масштаб давления может быть найден двумя путями:

2 —1

М- —2 Р —1 —2

= — , кг• м• с ; =-----, кг• м с . (17)

М/ рКА

Так как алфавит и правила переходов определены самой моделью, то для ее построения необходимо определить только линейные размеры и начальный клеточный массив, включая состояния клеток источников и стоков. Это делается следующим образом: 1) линейные размеры модели N1 , N , N вычисляются путем деления физических заданных величин /х, /у, 4 на р* 2) состояния источников и стоков определяются через их плотность путем деления значений внешнего давления на рр. Остальные масштабы необходимы при интерпретации результатов моделирования.

Пример 2. Моделирование потока моторного масла в трубе, частично перекрытой задвижкой.

Труба имеет диаметр Б = 0,7 м, длину / = 7 м, с одного торца трубы приложено давление Р = 7 000 кг-м-1. Другой торец открыт. Задвижка расположена на расстоянии 3,5 м от источника жидкости и перекрывает трубу до половины ее диаметра. Вязкость моторного масла V = 4,5 • 10-3 м2 с-1. Плотность моторного масла р = 880 кгм-3. Целью моделирования является поле скоростей за задвижкой. Для моделирования применяется 3D-модель RD-1. Масштабы плотности и вязкости определяются сразу:

= —= 13,8-10-3 м2 с-1, =-^ = 225 кгм3. (18)

УКА РкА

Масштаб давления можно определить исходя из модельной плотности рКА (табл. 2). Поскольку давление в источнике равно падению давления в трубе, модельную плотность в клетках источника можно принять равной рКА(0) « 2рКА. Пусть рКА(0) = 8. Тогда из (17) и (18) можно вычислить все остальные масштабы:

р p =

AP

рКА (0)

= 875 кг•м 1с 2,

= = 1,97 м • с-1,

Vi = — = 7 -1G-3

Vu

м, pt = — = 3,55 -10 3 c.

Vu

(l9)

Зная масштаб длины, можно определить размеры КА: диаметр БКА = х/р1 = 100 и длину трубы /КА = //р/ = 1000. Таким образом КА типа КЭ-1 определен.

Моделирование КА RD-1 производилось на кластере МВС-1000/М (МСЦ, Москва). Результаты моделирования получены в виде поля трехмерных векторов (^аХ*,у,к), которые равны осредненным значениям КА-скоростей в клетках (рис. 3). Соответствующие физические величины определяются путем умножения на масштаб скорости:

^х, у, х) = р/ • (^аХ/,У, к), X = /•р/, у = у р/ , 2 = к • р/ .

Максимальная модельная скорость потока над задвижкой =0,9, что со-

ответствует физической скорости итах = 1,77 м-с-1, а средняя скорость потока в трубе, вычисленная как осредненная проекция скорости на расстоянии /КА= 700 от источника, (^д) = 0,67, т.е. (и> = 1,3 м-с-1.

Рис. 3. Фрагмент проекции на продольный разрез трубы поля скоростей за задвижкой, полученный путем применения библиотеки MathCad к результату моделирования потока в трубе с помощью КА-модели КД-1

4. Масштабирование асинхронных вероятностных КА

Кинетические асинхронные КА (АКА) применяются для моделирования явлений, состоящих из элементарных действий, непосредственным образам имитирующих движения и взаимодействия молекул и атомов, называемых частицами. В настоящее время известно использование АКА-моделирования (называемого иногда методом Монте-Карло) для исследования процессов на поверхности катализаторов [9], эпитаксиального роста кристаллов [8]. Метод моделирования состоит в следующем. Область моделирования разделяется на клетки, на которых могут размещаться частицы веществ, участвующих в процессе. Частицы движут-

ся, взаимодействуют и претерпевают трансформации по законам, предписанным моделируемым явлением. АКА самым непосредственным образом отображает эти действия.

Алфавитом в АКА может быть: 1) булево множество А = {0,1}, если достаточно указать наличие частицы в клетке; 2) множество символов, обозначающих вещества (Al, CO, H2, 0, и т.д. ); 3) множество целых чисел, если много частиц могут находиться в одной клетке.

Множество координат M={(i,j): i,j = 0,...,N} часто расширяется добавлением в него контекстного пространства, например одной обобщенной клетки с именем Gas, которая символизирует внешнее пространство, откуда поступают или куда уходят частицы. Множество M имеет структуру той или иной кристаллической решетки. С помощью локальных операторов моделируют следующие наиболее часто используемые элементарные действия:

- Адсорбция. Частица ae A адсорбируется из газа на пустую клетку с вероятностью pa.

0a : {0, m )}* {a, Gas)}—{{a, m)}. (20)

- Десорбция. Частица b десорбируется с поверхности с вероятностью pb.

0а : {(«>m)} Pi >К0,m)}• (21)

- Реакция. Если частицы а и b оказываются в соседних клетках, между ними происходит реакция с вероятностью pab и образуется частица ab, которая, например, уходит в газ.

: {m))Ъ,фО))} РаЬ >{(0,т)(0,фО))}. (22)

- Диффузия. Если частица а оказывается рядом с пустой соседней клеткой, она

передвигается в нее с вероятностью pd .

©d : {a,m)(0,ф(m))} Pd >{(0,m)(a,ф(m))}. (22)

Когда ставится задача моделирования кинетического процесса с помощью АКА, задаются вещества, участвующие в процессе, условия, которые необходимо учесть (температура, давление), и начальные значения всех величин. Что касается размеров исследуемого образца, то они ограничиваются чаще вычислительными ресурсами, чем требованиями исследователя, для которого «чем больше, тем лучше».

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

1) Для химической реакции R e R, i = 1,...,n, где R множество реакций, участвующих в процессе, вероятность применения оператора 9/

к

P =~n----’

I

j=1

где kt , kj - скорости реакций в с-1.

2) Вероятность адсорбции вычисляется для атомов каждого вещества по заданному парциальному давлению в газе и силе связи частицы с поверхностью.

3) Вероятность диффузии зависит от сил связи между атомами и вычисляется исходя из энергетических величин.

Переход от модельных величин к физическим для длины и массы тривиален.

- Масштаб длины щ = I, м, где I - физический линейный размер атома (молекулы).

- Масштаб массы = g, кг , где g - физическая масса одного атома (молекулы).

- Масштаб времени р вычисляется на основе сопоставления одной итерации, исчисляемой как интервал между какими-то последовательными событиями, с известным значением какой-то временной характеристики процесса. Например [13], если адсорбируется один атом за итерацию и известен поток атомов Р, м-2с-1, поступающих на поверхность (количество атомов адсорбируемых на кв. метре площади в секунду ) или скорость осаждения атомов ЫЬ, с-1 (число слоев адсорбированных атомов в секунду), то

= ^ ^ ^, с> (24)

' ■ Р■ I м I ' Ж I

где - общее число адсорбированных частиц, гКА - число итераций в секунду, измеренное при моделировании, |Ы| - количество клеток в массиве.

Поскольку в моделировании участвуют реальные атомы, визуализация эволюции АКА позволяет наблюдать непосредственно физический процесс.

Пример 3. Упрощенная модель эпитаксиального роста силиконового ^) кристалла. Процесс состоит из двух действий: адсорбции Si атомов из внешнего потока и 2) диффузии адсорбированных атомов (адатомов) по поверхности кристалла. Адсорбция происходит с вероятностью ра =1. Скорость адсорбции в реальном эксперименте ЫЬ& 1,4 с-1. При адсорбции на поверхность слой за слоем адатомы образуют столбики и островки разной высоты и размеров. Верхний адатом в клетке (и,(г, у)) может диффундировать на соседнюю клетку Ф*(г,у) = (г+ау+Ъ), а,Ъ е {-1,0,1}, если и(г'+ау'+Ъ) < «(у). Вероятность такого действия зависит от энергии связи между соседними адатомами, характеризуемой константой В = 0,05 следующим образом. Если клетка (и,(г, у)) имеет N соседних адатомов, то вероятность диффузии р’А = (4 - п)-0,05". Выбор среди соседних клеток, куда адатому разрешено двигаться, равновероятен. Отсюда рЖк) = р&/(4 - п). Процесс моделировался АКА с алфавитом А = {0,1,2,...}, множеством имен Ы = {(г,у): г,у = 0,...,199}. Клетка (у,(г,у)) означает, что в точке (г,у) моделируемой области число адатомов (высота столбика) равна V е А. Локальный оператор © = 0<ш (г, у) 0^0', у) , где

ads

{о.0'. У))} * {(1. с™)} ——^ {(■(о +1. (»'. У))}.

: {(о.(г'>}))>Фк (*.І))} Р“(к) >(о -1.(*.І))+1,Фк (*.І)) (25)

к = 1,2, 3,4.

Всякий раз клетка (г, у) выбирается случайным образом.

Моделирование выполнялось на компьютере Репіїиш IV. Среднее время итерации М = 5,2-10-3 с. На рис. 4 видны сформированные в процессе эволюции островки адатомов. Важным параметром, характеризующим свойства эпитаксиального, покрытия является суммарный периметр островков Р(і), который вычисляется в процессе эволюции (рис. 4, Ъ). Масштаб времени р, = (М^|М|)-1 = 0,178-10-4 с.

**к я

♦V.f^V.v;

V. : t . 4*ь у:« м

f*«: >»irf а * v

Л4.»

Nd

10000-

8000-

6000"

4000"

2000-

1

г = 100 000 2,15 4,02 с

Рис. 4. Глобальное состояние и зависимость периметра островков от времени в примере 3 при моделировании эпитаксиального роста кристалла

Заключение

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

ЛИТЕРАТУРА

1. Toffolli T. Cellular automata as an alternative to (rather than approximation of) differential equations in modeling physics // Physica D. 1984. V. 10. P. 117 - 127.

2. Wolfram S. A new kind of science. Champaign, 1ll., USA: Wolfram Media Inc., 2002. 1200 p.

3. Бандман О.Л. Клеточно-автоматное моделирование пространственной динамики // Системы информатики. Новосибирск: Изд-во СО РАН. Вып. 10. С. 59 - 113.

4. Малинецкий Г.Г., Степанцов М.Е. Моделирование диффузионных процессов клеточными автоматами с окрестностью Марголуса // Журнал вычислительной математики и математической физики. 1998. № 6. С. 1 - 17.

5. Rothman D. H., Zalesky S. Lattice-Gas Cellular Automata. Simple Model of Complex Hydrodynamics. UK: Cambridge University press, 1997. 297 р.

6. Succi S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. N.Y.: Oxford University Press, 2001. 240 p.

7. Frish U., d'Humieres D., Hasslacher B., LallemandP., Pomeau Y., Rivet J.P. Lattice-Gas hydrodynamics in two and three dimensions // Complex Systems. 1987. №. 1. Р. 649 - 707.

8. Neizvestny I.G., Shwartz N.L., Yanovitskaya Z.Sh., Zverev A. V. 3D-model of epitaxial growth on porous {111} and {100} Si Surface. // Comp. Phys. Comm. 2002. V. 147. P. 272 - 275.

9. Elokhin V.I., Latkin E.I., Matveev A.V., Gorodetskii V.V. Application of statistical lattice models to the analysis of oscillatory and autowave processes in the reaction of carbon monoxide oxidation over platinum and palladium surfaces // Kinetics and Catalysis. 2003. V. 44. №. 5. С. 672 - 700.

10. Achasova S., Bandman O., Markova V., Piskunov S. Parallel substitution algorithm. Theory and Application. Singapore: World Scientific, 1994. 180 р.

11. Sahimi M. Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata and simulated annealing // Rev. Modern Physics. 1993. V. 65. Nc. 4. P. 1533 -1660.

12. Медведев Ю.Г. Трехмерная клеточно-автоматная модель вязкой жидкости // Автометрия. 2003. Т. 39. № 3. С. 43 - 48.

13. Neizvestny I. G., Shvarts N.L., and Yanovitskaya Z. Sh. Two-dimensional epitaxial nucleation with large critical-nucleus size // Russian Microelectronics. Science Business Media LLC. 2002. V. 31. Na 2. P. 70 - 78.

Статья представлена кафедрой информационных технологий в исследовании дискретных структур радиофизического факультета Томского государственного университета и оргкомитетом 7-й Российской конференции с международным участием «Новые информационные технологии в исследовании сложных структур», поступила в научную редакцию 10 мая 2008 г.