Моделирование движения поршня в газовой среде клеточным автоматом Текст научной статьи по специальности «Математика»

2010 Дискретные модели реальных процессов №4(10)

УДК 004.007.519-7

МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ПОРШНЯ В ГАЗОВОЙ СРЕДЕ КЛЕТОЧНЫМ АВТОМАТОМ1

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

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

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

E-mail: medvedev@ssd.sscc.ru

Обсуждаются способы задания граничных условий в многочастичной клеточноавтоматной модели потока FHP-MP. Особое внимание уделено нестационарным граничным условиям. Изложен метод перемещения стенок в исследуемой модели.

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

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

Введение

Класс клеточных автоматов, используемый для моделирования процессов газовой и гидродинамики, называется Lattice Gases [1], или, как пишут в русскоязычной литературе, «решёточные газы». В рамках настоящей работы этот класс назван газодинамическими клеточными автоматами. Впервые подобный автомат, именуемый по инициалам создателей HPP, был предложен Харди и др. в 70-х годах прошлого века [2]. Его квадратные клетки имеют по четыре соседа, а поведение детерминировано. В силу небольшого числа правил столкновения он обладает очень узкими границами применимости и для моделирования потоков практически не используется. Спустя десять лет Фриш и др. предложили другой клеточный автомат — FHP, ставший основой для многих последующих моделей клеточно-автоматной газодинамики [3]. Этот автомат использует достаточно богатый набор вероятностных правил столкновения. Каждая его клетка является правильным шестиугольником и имеет шесть соседей. Еще через пятнадцать лет был предложен трехмерный газодинамический клеточный автомат [4, 5]. Он был назван RD, так как его элементарный объем имеет форму ромбического додекаэдра. В этом многограннике двенадцать граней, следовательно, каждая клетка имеет двенадцать соседей. Поведение этого автомата вероятностное.

Граничные условия в газодинамических клеточных автоматах задаются путём назначения некоторым (граничным) клеткам специальных правил столкновения. Например, все частицы, попадающие в один из граничных типов клеток — стенки, отскакивают назад. Из таких клеток-стенок выстроены препятствия, которые поток вынужден обтекать. Также их обычно используют для ограждения моделируемого объема по периметру. При осреднении значений оказывается, что скорость потока на границах таких препятствий равна нулю, что совпадает с физическими представлениями о потоке.

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

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

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

Рис. 1. Коллизии при перемещении стенок в модели РИР с булевым алфавитом

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

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

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

1. Основные определения

Многочастичный клеточный автомат ЕНР-МР — это тройка объектов (К, Д, 0), где К = {с^ с2,... , с,,...} — множество клеток, заданное их индексами г и ] в некотором дискретном пространстве (см. рис. 1). Для каждой клетки с] Е К определено упорядоченное множество соседства N(с^) = {пк(с^) € К : к = 0,1,...,6}, определяемое следующим образом:

д( ) = / {с,,] , с,—1,] , Сг-1^+1,Сг^+1,Сг+1^, ст,]—1, с,,]—1} при нечётном ] ,

д (С]) = 1 г т (1)

^{с*,], С*-1,], С*,]+1, С,+1,]+1, сг+1,], Сг+1,] — 1, сг — 1,] — 1} при чётном ].

Каждая клетка с Е К характеризуется состоянием б(с), которое зависит от дискретного времени £ и представляет собой вектор, имеющий целочисленные компоненты в к (с), к = 0,1,... , 6, указывающие на количество модельных частиц газа в клетке с с единичной массой и вектором скорости ек (с), направленным в сторону соседа Пк(с). Вектор скорости ео (с) имеет нулевую длину, поэтому частицы во (с) принято называть частицами покоя. Остальные векторы имеют единичную длину, так что частицы в каждом направлении несут импульс рк (с) = в к (с) ек (с). Множество состояний б (с) всех клеток с Е К в один и тот же момент времени £ называется глобальным состоянием П (£) = {б (с,]) : V*, ]} клеточного автомата.

Для каждой клетки с,] Е К определены координаты её центра х(с]) и у(с]) на декартовой плоскости следующим образом:

при нечётном ],

х (с]) = ^ . 1

г + ^ при чётном (2)

( , ^3.

У(с]) = -у >

Расстояние й (с^], с,2]2) между клетками с^] Е К и с,2]2 Е К обусловлено соотношением

й (сПЛ , с,2]2) (х (сПЛ ) — х (с,2]2)) + (у (с»и'1) — У (с,2]2 )) . (3)

Осредненные значения скорости и являются модельными значениями скорости потока, а осредненные значения плотности р являются модельными значениями давления. Они вычисляются для некоторой окрестности (с) = {с,] Е К : й (с,], с) ^ г}, где

r — радиус осреднения. Осредненная скорость вычисляется как отношение суммарного импульса всех частиц в клетках cj Е Av (с) к суммарной массе этих частиц:

6

S S Pfc (cij)

Cij€Av(c) fc=0

u (c) =----------------------------------------------------6-. (4)

S S sk (c»j)

Cij GAv(c) k=0

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

1 6

Р(c) = 14 ( )I S S Sfc (с^^^ (5)

|Av (c)| CijeAv(c) k=0

где |Av(c)| — мощность окрестности осреднения Av(c).

2. Стационарные граничные условия в модели FHP-MP

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

Клеточный автомат FHP-MP функционирует в синхронном режиме. На каждой итерации происходит смена состояний s (t) всех клеток с Е K на состояния s (t + 1) = = в (s (t)), где в (s (t)) Е 0 — соответствующая функция переходов клетки с. Клеточный автомат при этом переходит из глобального состояния Q(t) в новое глобальное состояние ^(t + 1). Для того чтобы модель адекватно отображала физический процесс, функция в в клетках среды с Е Kc должна удовлетворять законам сохранения массы

Е Ев (sfc(с)) = Е Е sfc(с) (6)

c€Kc k=0 c€Kc k=0

и импульса

S S в (Pk (с)) = S S Pk (с). (7)

c€Kc i=1 c€Kc k=1

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

в (s) = в2 (в1 (s)) . (8)

На стадии сдвига в каждой клетке с Е K каждая частица, учтённая в компоненте sk (с), k = 1,... , 6, вектора состояния s (с) перемещается в соседнюю клетку nk(с),

расположенную на единичном расстоянии в направлении вектора её скорости ek (с).

Частицы покоя, соответствующие компоненте s0 (с), остаются в клетке с. Таким образом, k-е компоненты sk (с) вектора состояния s (с) клетки с после сдвига принимают значения

вх (sfc (с)) = /sk (N((k+2) mod 6) + 1 (с)) для k = 1, 2,..., 6, (д)

I sk (с) для k = 0.

На стадии столкновения 92 в каждой клетке с Е К происходит изменение направления движения частиц согласно некоторым правилам столкновения, зависящим только от состояния и типа клетки с и не зависящим от состояний соседних клеток из N (с).

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

Е #2 (^ (с)) = Е ^ (с), с Е Кс,

к=0 к=0 (10)

Е 92 (Рк(с)) = Е Рк(c), с Е Кс. к=1 к=1

В клетках с Е Кад, являющихся стенками, частицы «отражаются» в обратном на-

правлении, нарушая при этом закон сохранения импульса:

в2 (вк (с)) = (5((к+2) т°а 6)+1 (с™) для к = 1 2,..., 6’ с Е К«,. (11)

^ (сад) для к = 0,

Такое поведение частиц в клетках-стенках моделирует условие нулевой скорости потока на границах препятствий.

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

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

#2 (^ (с)) = рк, с Е К5. (12)

Граница клеточного автомата, выстроенная из таких клеток, моделирует открытый в атмосферу край объекта.

3. Нестационарные граничные условия в модели FHP-MP

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

В случае нестационарных граничных условий функция переходов 9 состоит из композиции трёх функций 91 (сдвиг), 92 (столкновение) и 93 (перемещение границ):

9 (з) = 93 ( 92 (91 (8))м) , (13)

К1

где М =--------отношение абсолютной величины скорости движения поршня к скоро-

сти звука в газе в модельных единицах. Фактически, М является модельной величиной, обратной числу Маха. Итерацией здесь будем называть однократное выполнение

$2 ($і (з)). Стадии сдвига и столкновения проходят так же, как и со стационарными граничными условиями.

На стадии перемещения границ подвижная стенка (поршень) перемещается на одну клетку на каждой М-й итерации. При этом через поршень не должно просочиться ни одной частицы. На рис. 2, а изображен фрагмент клеточного автомата до перемещения поршня, располагающегося в столбце і = 1 с некоторым набором частиц в клетках. В каждом направлении может быть ориентировано несколько (0,1, 2,...) векторов скорости частиц. На рисунке все такие векторы, сколько бы их ни было, изображены одной стрелкой. На рис. 2, б изображен тот же фрагмент после перемещения поршня на одну клетку вправо. Клетки столбца і = 1 стали клетками среды, клетки столбца і = 2 — стенками. Частицы столбца і = 2 переместились в столбец і = 3, добавившись к тем частицам, которые уже там были, а на их место переместились частицы столбца і = 1 . Коллизий не возникло. Таким образом, функция перемещения границ имеет вид

$3 (^ (Сі,з)) = ^ (Сі,з) + ^ (Сі,2) ,

$3 (Зк (Сі,2)) = ^ (Сі,і) ,

$3 (Зк (Сі,і)) = 0, (14)

Сг,2 переходят из Кс в Кад,

Сг,2 переходят из Кад в Кс.

Рис. 2. Перемещение стенки в модели РИР-МР с целочисленным алфавитом

Перемещение частиц из столбца і = 2 в столбец і = 3 влечет за собой то, что в стенке не окажется ни одной частицы с вектором скорости, направленным влево, вверх или вниз, что, в свою очередь, гарантирует то, что ни одна из частиц не проникнет сквозь поршень. Образовавшийся за поршнем в столбце і = 1 вакуум займут частицы, переместившиеся туда слева (при их наличии) за следующие М итераций, пока поршень покоится.

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

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

4. Компьютерное моделирование

При проведении вычислительных экспериментов ставилась задача посмотреть качественную картину процесса, получив только модельные значения количественных характеристик, не вычисляя их в физических единицах. За единицу длины принято расстояние между соседними клетками. Единицей массы служит масса одной частицы. Скорость потока в точке вычисляется по формуле (4), а давление газа в точке — по формуле (5), так как оно пропорционально плотности модельных частиц. Код программ написан на языке Си. Хранение состояния автомата осуществляется в двумерном массиве. Из-за синхронного режима функционирования приходится использовать два таких массива. Функция сдвига #1 использует в качестве аргумента значения из одного из них, а результат помещает в другой. Для записи результата массивы используются поочередно: один на чётных, а другой на нечётных итерациях. Функции столкновения #2 и перемещения границ записывают свои результаты в тот же массив, в котором находятся их аргументы, и применяются к каждому из двух массивов поочередно, после того как функция #1 запишет в него свой результат.

Программным ограничением является максимальное количество частиц в одной клетке, имеющих одинаковый вектор скорости е^ (с), равное 255, вследствие того, что каждый из семи компонентов (с) вектора состояния хранится в одном байте памяти. Еще один байт занимает информация о типе клетки. Таким образом, состояние каждой клетки занимает 16 байт памяти с учетом дублирования во втором массиве. К примеру, клеточный автомат, эксперименты с которым описаны ниже, имеющий размер 100x200 клеток, требует около 0,3 мегабайта памяти.

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

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

Функция #2 последовательно перебирает все клетки автомата по следующим правилам. Если клетка имеет тип «стенка», то векторы скорости всех частиц в ней меняются на антипараллельные (11). Клетки-источники поддерживают плотность частиц газа на заданном уровне (12). Клетки среды производят столкновение по следующему алгоритму. Вначале вычисляются концентрация частиц в клетке и проекции их импульсов на декартовы оси. Затем полным перебором подсчитывается количество

состояний, сохраняющих массу и импульс, обусловленные выражениями (10). После этого равновероятно выбирается одно из этих состояний; его и полагают результатом функции #2. Из-за полного перебора столкновение является самой затратной вычислительной функцией.

Функция #3 производит замену типов клеток по нижеописанному алгоритму, если номер итерации кратен М, или ничего не делает — в противном случае. Во-первых, частицы газа из клеток, куда будет перемещен поршень, перемещаются на одну клетку в направлении движения поршня и с учетом векторов скорости добавляются к тем частицам, которые там уже были. Во-вторых, частицы газа, находящиеся в клетках поршня, перемещаются на одну клетку в направлении движения поршня на освободившееся место. В-третьих, клетки поршня меняют тип с клеток-стенок на клетки среды. И, в-четвертых, клетки, куда переместился поршень, меняют тип с клеток среды на клетки-стенки (14).

После этого происходит следующая итерация и т. д. Значение М для эксперимента 1 выбрано равным 20, а в экспериментах 2 и 3 — равным 10. В каждом эксперименте вначале поршень двигался вправо 100 раз (это соответствует 2000 итерациям в эксперименте 1 и 1000 итерациям в экспериментах 2 и 3). Затем поршень оставался в покое, т. е. следующие 30 перемещений не выполнялись. И, наконец, поршень совершал 100 перемещений влево, возвращаясь в начальное положение.

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

Заключение

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

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

Рис. 3. Результаты компьютерного моделирования. Поля скорости и давления

ЛИТЕРАТУРА

1. RothmanD.H., ZaleskiS. Lattice-gas cellular automata: simple models of complex

hydrodynamics. Cambridge University Press, 1997.

2. Hardy J, Pomeau Y., and de Pazzis O. 2D Lattice-Gas mode // J. Math. Phys. 1973. No. 14. P. 1746.

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

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

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

6. Медведев Ю. Г. Многочастичная клеточно-автоматная модель потока жидкости FHP-MP // Вестник Томского госуниверситета, сер. «Управление, вычислительная техника и информатика». 2009. №1(6). С. 33-40.