Клеточно-автоматная модель формирования порошковой струи Текст научной статьи по специальности «Математика»

2009 Вычислительные методы в дискретной математике №3(5)

УДК 004.007.519-7

КЛЕТОЧНО-АВТОМАТНАЯ МОДЕЛЬ ФОРМИРОВАНИЯ ПОРОШКОВОЙ СТРУИ1

Ю. Г. Медведев

Институт вычислительной математики и математической геофизики СО РАН,

г. Новосибирск, Россия

E-mail: medvedev@ssd.sscc.ru

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

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

Ключевые слова: клеточный автомат, поток газа, порошковая струя.

Введение

В большинстве классических клеточно-автоматных моделей потоков используется один тип частиц — частицы газа [1]. Двумерная модель, ставшая основой для множества новых моделей, называется FHP (Frish, Hasslacher, Pomeau) [2]. Каждая ее клетка имеет форму шестиугольника и, в силу идемпотентности отношения соседства, имеет семь соседей. В качестве внутренних состояний клеток модели используются булевы векторы. Состояние клетки — это набор модельных частиц газа с единичной массой. Частицы дискретно перемещаются из клетки в клетку через равные интервалы времени. Каждая частица либо имеет вектор скорости с единичным модулем, направленный в сторону одной из шести нетождественных соседних клеток, либо имеет нулевую скорость. В клетке одновременно не могут находиться две или более частицы с одинаковыми скоростями. Каждый из семи разрядов вектора состояния клетки указывает на наличие или отсутствие в этой клетке частицы с соответствующей ему скоростью. Следовательно, одновременно в клетке могут находиться от нуля до семи частиц. Для расширения возможностей модели FHP была создана многочастичная модель FHP-MP (multi-particle) [3]. Она допускает одновременное присутствие в клетке более одной частицы с одинаковой скоростью. Эта модель позволяет работать в более широком диапазоне давлений, так как максимальное количество частиц в клетке ограничено только программной реализацией. Кроме того, она более пригодна для создания композиций

хРабота выполнена при поддержке интеграционного проекта СО РАН №32, 2009 г. и проекта 2.6 программы Президиума РАН, 2009 г.

с другими клеточно-автоматными моделями [4]. Построению и исследованию одной из таких композиций и посвящена эта работа. Композиция построена для моделирования газовой струи, несущей в себе порошок, путем добавления в модель ЕИР-МР частиц порошка и задания новых правил поведения. Новая композиционная модель с двумя типами частиц — газ и порошок — названа ЕИР-ОР (gas-powder).

В статье описана новая клеточно-автоматная модель порошковой струи ЕИР-ОР, приведены результаты вычислительных экспериментов по моделированию порожденной взрывом струи газа, несущей порошок, произведено сравнение с результатами моделирования струи газа без порошка, полученными в результате вычислительных экспериментов с использованием модели ЕИР-МР.

1. Описание модели ЕИЕ-СЕ

Под клеточным автоматом модели ЕИР-ОР будем понимать тройку объектов (К, N 0), где К = {сі, с2,... , с^ ... } — множество клеток, заданное их координатами в некотором дискретном пространстве. Каждая клетка с Є К характеризуется состоянием в (с) и координатами х(с) и у(с) на Декартовой плоскости. Состояние в(с) зависит от дискретного времени і, а координаты остаются неизменными. Следовательно, между любыми двумя клетками с1 Є К и с2 Є К с легкостью можно подсчитать расстояние ^(с1, с2). В модели ЕИР-ОР состояние в(с) — это не один вектор, как у моделей-предшественников, а два: для газовых частиц целочисленный, как в ЕИР-МР, и для частиц порошка булев, как в ЕИР. Такая композиция является новой клеточноавтоматной моделью.

Для каждой клетки с Є К определено упорядоченное множество из семи элементов N (с) = {пДс) : п0(с) = с, Пі (с) Є К & ^(с, пДс)) = 1,і = 1, 2,..., 6}. Входящие в него клетки Пі (с) находятся в отношении соседства с клеткой с и называются ее соседними клетками или соседями. Таким образом, структура множества клеток К клеточного автомата представляется неориентированным графом, в котором вершинами являются клетки, а ребрами — отношение соседства. Этот граф имеет регулярную структуру и степени вершин, равные семи.

Состояние 5 клетки с Є К есть множество 5 (с) = (с) , Бр (с)}, состоящее из двух

векторов. Первый из них (с) имеет целочисленные компоненты в^ (с) , і = 0,1,..., 6, указывающие на количество частиц газа в клетке с с единичной массой и единичным вектором скорости еі (с), направленным в сторону соседа пі(с). Эти частицы несут импульс р^ (с) = в^ (с) еі (с). Второй вектор Бр (с) = (вр (с) , вр (с) ,... , вр (с)) является булевым, его компоненты определяют наличие или отсутствие в клетке с частиц порошка с вектором скорости еі (с) и импульсом рР (с) = МР вр (с) еі (с), где і = 0, 1,..., 6, а Мр — масса одной частицы порошка. В направлении каждого из соседей пі(с) может быть ориентировано не более одного вектора скорости. Множество состояний в (с) всех клеток с Є К в один и тот же момент времени і называется глобальным состоянием П(і) = {в(с1), в(с2),... , в(сі),... } клеточного автомата.

На рис. 1 изображена клетка с, единичные векторы еі (с) находящихся в ней частиц газа и порошка и ее соседи пі(с), і = 0,1,... , 6. Вектор скорости еі (с) либо направлен в сторону одной из соседних клеток пі(с) (при і = 1,... , 6), либо равен нулю (при і = 0). Масса частиц, движущихся в направлении еі (с), равна ті (с) = в^ (с) + Мрвр (с).

Компоненты в^ (с) , в^5 (с) ,... , в^ (с) газовой составляющей вектора состояния в° клетки с принимают целочисленные значения. Суммарная масса частиц газа в клетке с

равна

6

т° (с) = Е 5? (')■ (!)

*=0

где ер — г-й компонент вектора состояний зс. Масса частицы порошка равна Мр, поэтому масса всех частиц порошка в клетке с равна

6

Шр (с) = МР Е ^ (с). (2)

г=0

Рис. 1. Векторы скорости частиц, нумерация соседей Общая масса частиц газа и порошка в клетке с равна

т (с) = тс (с) + тр (с).

(3)

Импульс частиц, направленных к соседу п*(с), р* (с) = рр (с) + рр (с) в клетке с Е К, есть сумма импульсов газа рр (с) и порошка рр (с). Суммарный по всем направлениям импульс равен

6

р = 2^ р*

(4)

*=1

Из (4), учитывая рис. 1, легко подсчитать проекции рх и ру импульса р (с) на Декартовы оси Ох и Оу:

(|р2| + |р3| - |р5| - |р6|) ,

Ру = |р4| - |р1| + 2 (|рЗ| + |р5| - |р2| - |рб|) .

(5)

(6)

Определено разбиение клеток с Е К на типы. Клетками среды сср Е Кср называются клетки, в которых выполняются законы сохранения массы и импульса. Клетки стенок сст Е Кст — это клетки, в которых выполняется закон сохранения массы, но может нарушаться закон сохранения импульса. И наконец, источники сист Е Кист — клетки, в которых могут нарушаться как закон сохранения массы, так и закон сохранения импульса. Множества клеток среды Кср, стенок Кст и источников Кист попарно

не пересекаются (Кср П Кст = 0, Кср П Кист = 0, Кст П Кист = 0). Объединение этих множеств совпадает с множеством всех клеток автомата (Кср иКст и Кист = К). Поведение стенок и источников задает граничные условия клеточного автомата.

В модели ЕИР-ОР используется клеточный автомат с синхронным режимом функционирования. На каждой итерации происходит смена состояний в (¿) всех клеток с Е К на состояния в (£ +1) = в (в (¿)), где в (в (¿)) Е 0 — соответствующая функция переходов клетки с. Клеточный автомат при этом переходит из глобального состояния П(£) в новое глобальное состояние П(£ + 1). Для того чтобы модель адекватно отображала физический процесс, функция в в клетках среды сср Е Кср должна удовлетворять законам сохранения массы газа

ЕЕ«(8? Ы) = ££(сср)

сЄК ¿=0

сЄК ¿=0

массы порошка

ЕЕ 9 (^ (сср)) = ЕЕ ^ (сср)

сЄК ¿=0

сЄК ¿=0

и общего импульса

ЕЕ9 (рі (сср)) = ЕЕ рі (сср)-

(9)

с€К ¿=1 с€К ¿=1

Каждая итерация выполняется в две фазы: сдвиг и столкновение. Функция переходов в состоит, таким образом, из суперпозиции функций в1 (сдвиг) и в2 (столкновение):

9 (в) = 92 (9! (в)).

ПО)

В фазе сдвига в каждой клетке с Е К каждая частица, учтенная в компонентах вр (с) и вр (с), г = 1,... , 6, векторов состояния (с) и бр (с) перемещается в соседнюю клетку ^¿(с), соответствующую ее вектору скорости е^ (с). Частицы, соответствующие компонентам вр (с) и вр (с), остаются в клетке с. Таким образом, г-е компоненты вр (с) и вр (с) векторов состояния (с) и бр (с) клетки с после сдвига принимают значения

9і (дс (с)) 9і (¿Р (с))

^ (N((¿+2) шса 6) + 1 (с)) для і = 1, 2, . . . , 6, (с) для і = 0;

(N((¿+2) шса б)+1 (с)) для і = 1, 2,..., 6,

Г11)

П2)

вр (с) для г = 0.

Несмотря на то, что при сдвиге масса и импульс частиц в клетке изменяются, в пределах всего клеточного автомата они сохраняются, т. е. условия (7), (8) и (9) выполняются.

В фазе столкновения в2 происходит изменение направления движения частиц согласно некоторым правилам столкновения, не зависящим от состояний соседних клеток. В модели ЕИР-ОР функция в2 вероятностная. Ниже описаны правила столкновения для клеток разных типов: среды, стенок и источников.

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

6

@2 (« (Сср^ ^ ] « (Сср), сср 6 К; (13)

i=0 i=0

6 6

Е @2 («[ (сср)) = Е «Р ^Р^ СсР 6 К; (14)

i=0 i=0

66

Е@2 (Pi (Сср)) = Е ^ (Сср), СсР 6 К. (15)

i=1 i=1

При выполнении условий (13), (14) и (15) условия (7), (8) и (9) тем более выполняются.

В клетках Сст 6 Кст, являющихся стенками, частицы «отражаются» в обратном направлении, нарушая при этом закон сохранения импульса:

@2 («Т (Сст)) = { :^^‘(!СТ0. ^ “ = "’ 2' . . . ' 6' (16)

Из-за того, что количество частиц в клетке не меняется, условия (13), (14), а следо-

вательно, и (7), (8) выполняются. Условие (15) может нарушаться, так как меняются направления векторов скорости частиц, что соответствует изменению импульса частиц в газовой среде при столкновении с твердыми стенками. Такое поведение частиц в клетках-стенках моделирует условие нулевой скорости потока на границах препятствий.

Каждая клетка-источник Сист 6 Кист по алгоритму, определяемому граничными условиями, генерирует частицы со всевозможными направлениями вектора скорости. Из клеток-источников можно создавать разные объекты, задавая различные способы введения газа в моделируемый объем. Например, установив такие клетки в пространстве в одну линию (как правило, у границы клеточного массива), можно получить источник равномерного потока частиц заданной плотности. Отдельно установленный источник будет моделировать форсунку. Естественно, при генерации новых частиц ни масса шс (сист), ни импульс рс (сист) частиц газа не сохраняются.

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

Е Е Pi(с)

/ \/ Ч с£Аь(со)

(и) (с°) = ------6------. (17)

Е Е Ш (с)

сбА'и(со) i=0

Осредненная плотность частиц (рс) и (рр) подсчитывается в той же окрестности Аи(с°) следующим образом:

1 6

р°>(С°)= Е Е«?(с); (I8)

1 ( 0)1 сеА^(со) i=0

1 6

рР> (с0) = 1А^Ы| Е Е «Р (с)' (19)

1 1 0) 1 сеА^(со) i=0

где |Ау(с0)| — количество клеток, попадающих в окрестность осреднения Ау(с0).

Осредненные значения скорости (и) являются модельными значениями скорости газового потока, осредненные значения плотности (рР) являются модельными значениями плотности порошка, а осредненные значения плотности (рс) являются модельными значениями давления газа.

Следует отметить, что осредненные значения модельных скорости, плотности и давления соответствуют их физическим аналогам только в том случае, когда окрестность осреднения Ау(с0) состоит исключительно из клеток среды с 6 К. В противном

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

2. Результаты тестовых исследований модели

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

Хранение состояния автомата осуществляется в двумерном массиве. Из-за синхронного режима функционирования приходится использовать два таких массива. Функция сдвига @1 использует в качестве аргумента значения из одного из них, а результат помещает в другой. Для записи результата массивы используются поочередно: один на четных, а другой на нечетных итерациях. Функция столкновения @2 записывает свой результат в тот же массив, в котором находятся ее аргументы, и применяется к каждому из двух массивов поочередно после того, как функция @1 запишет в него свой результат.

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

интерпретируемый как булев вектор, семь из восьми разрядов которого отвечают за шесть направлений и одну частицу покоя, а восьмой разряд не используется. Еще один байт занимает информация о типе клетки. Таким образом, состояние каждой клетки занимает 18 байт памяти с учетом дублирования во втором массиве. К примеру, клеточный автомат, эксперименты с которым описаны ниже, имеющий размер 200x400 клеток, требует около полутора мегабайт памяти.

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

ция отработала, указатели на массивы обмениваются своими значениями, выполняется функция @2, а затем следующая итерация.

Функция @2 последовательно перебирает все клетки автомата по следующим правилам. Если клетка имеет тип «стенка», то векторы скорости всех частиц в ней меняются на антипараллельные. Клетки-источники генерируют новые частицы газа с заданной вероятностью. Клетки среды производят столкновение по следующему алгоритму. Вначале вычисляются концентрация частиц в клетке и проекции их импульсов на декартовы оси. Затем полным перебором подсчитывается количество состояний, сохраняющих массу газа, массу порошка и общий импульс, обусловленные выражениями (13), (14) и (15). После этого при помощи генератора случайных чисел, описанного в [5], равновероятно выбирается одно из этих состояний; его и полагают результатом функции @2. Из-за полного перебора столкновение является самой затратной вычислительной функцией, несмотря на ряд оптимизаций, дающих выигрыш в скорости на несколько десятичных порядков.

Ниже описаны вычислительные эксперименты, проведенные с моделью ЕИР-ОР. В первом из них используется только газовая составляющая модели и не используются частицы порошка. Такие начальные условия редуцируют модель до ЕИР-МР. Во втором эксперименте используются частицы как газа, так и порошка. Размеры автомата, топология стенок и другие параметры в обоих экспериментах одинаковы, что дает возможность сравнить их результаты.

Клеточный автомат, использующийся в этих вычислительных экспериментах, имеет размеры 200x400 клеток (вдоль декартовых осей Ох и Оу соответственно). Периметр с координатами в интервалах [(1, 1), (1, 200)], [(1, 1), (400, 1)], [(1, 200), (400, 200)] и [(400, 1), (400, 200)] выстроен из клеток-стенок. Также из стенок с координатами [(100, 1), (100, 80)] и [(100, 121), (100, 200)] построено сопло. Остальные — клетки среды. Так как моделируемый объем замкнут, клетки-источники не используются.

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

В начале обоих экспериментов в клетках среды находятся частицы газа со средней плотностью (рг) = 3 частицы на каждое направление, или 21 частица в клетке. Для моделирования взрыва в левой части камеры расположены две области газа с плотностью (рг) = 60 частиц покоя в клетке. Они имеют форму наклонных полос толщиной по 40 клеток и изображены в столбцах «Давление газа» черным цветом. Кроме того, во втором эксперименте внутри полос газа высокого давления в столбце «Концентрация порошка» имеются две наклонные полосы порошка с плотностью (рр) = 1 частица покоя в клетке. Эти полосы имеют толщину по 10 клеток. Масса частицы порошка МР = 20.

Результаты первого эксперимента после 70, 220 и 390 итераций показывают, что в левой части камеры от участков высокого давления навстречу друг другу идут

Эксперимент 1 (газ без порошка) Эксперимент 2 (композиция газа и порошка)

Давление газа Скорость потока Давление газа Концентрация порошка

о К £ £ Ф

° £ К о N N

« к Г-

Я" ей Он С-) К Ш 1 \У,' /,и [ и ► ■ і|і ¡и,

О ІМІ

к і ь ■ ■> і 11! ЧВ >|*| ’ ІА?» у’.' і

С-) Ё Ф I ШШ

О * ¡¡ДРу "У £%, чі І ■

СХІ сча ^ V V І V \ ^

« К =г сЗ Он 1) н к о Ш ■ ::::::::: ' -'.Г' : •. • .

05 со ^:Г::г':: Мі :

Рис. 2. Формирование струи. Результаты вычислительных экспериментов

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

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

Заключение

Проведенные эксперименты показали, что новая клеточно-автоматная модель ЕЫР-ОР способна описать процессы формирования порошковой струи. Конечно же, в статье показан только принцип применения модели. Для ее практического использования следует производить привязку модельных величин к физическим. Это необходимо делать каждый раз при настройке модели на новые размеры камеры, плот-

ность и степень намола порошка, температуру газа, мощность взрывчатки и другие параметры. Кроме того, в модели не учитывается гравитация и ряд других физических явлений, но при необходимости их довольно легко учесть. Также легко перейти к трехмерному варианту предложенной модели, используя клеточный автомат с пространственной структурой в форме ромбического додекаэдра, подобно тому, как был осуществлен переход от двумерной модели FHP к трехмерной RD [7]. Конечно, получившаяся модель RD-GP будет требовать очень много вычислительных ресурсов, но если закон Мура будет продолжать работать, то в скором будущем суперкомпьютеры смогут покрыть эти потребности.

ЛИТЕРАТУРА

1. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика. 2006. №10. С. 59-113.

2. Frisch U., Hasslacher B., and Pomeau Y. Lattice-Gas automata for Navier-Stokes equations // Phys. Rev. Lett. 1986. No. 56. P. 1505.

3. Медведев Ю. Г. Расширение клеточно-автоматной модели потока FHP-I до многочастичного варианта FHP-MP // Новые информационные технологии в исследовании сложных структур. Томск, 2008. С. 73.

4. Bandman O. L. Cellular Automata composition technuques for spatial dynamics simulation // Bull. Nov. Comp. Center, Comp. Science. 2008. V. 27. P. 1-39.

5. Kalgin K. V. Acceleration of linear congruental generators // Bull. Nov. Comp. Center, Comp. Science. 2008. V. 27. P. 51-53.

6. Кабулашвили В. Г., Тришин Ю. А. Струйные течения при взрывном обжатии пористых облицовок и их применение // Взрыв, удар, защита: Сб. научных трудов. Новосибирск, 1988. №18. С. 149.

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