Многочастичная клеточно-автоматная модель потока жидкости FHP-MP Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2009 Управление, вычислительная техника и информатика № 1(6)

УДК 519.68

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

МНОГОЧАСТИЧНАЯ КЛЕТОЧНО-АВТОМАТНАЯ МОДЕЛЬ ПОТОКА ЖИДКОСТИ FHP-MP

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

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

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

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

1. Описание модели FHP-MP Основные определения

Под клеточным автоматом CA модели FHP-MP будем понимать тройку объектов (W, A, N), где W = {w1, w2, ..., w, ...} - множество клеток, заданное их координатами в некотором дискретном пространстве. Каждой клетке w є W поставлен в соответствие конечный автомат A , называемый элементарным автоматом.

Внутренним состоянием автомата является целочисленный вектор, а не булев вектор как в классической модели БИР. Каждой клетке у є Ж сопоставлены некоторые координаты х(у) и у(у) на декартовой плоскости. Следовательно, между любыми двумя клетками у1 є Ж и у2 є Ж можно подсчитать расстояние сС(м>1, у).

Для каждой клетки V є Ж определено некоторое упорядоченное множество МУ) = {МУ): Мо(у) = у, МУ) є Ж & сС(у, МУ)) = 1, (і = 1, 2, ..., Ь)}, элементы которого находятся в отношении соседства с клеткой у и называются ее соседними клетками, или соседями. Константа Ь определяет количество нетождественных соседей каждой клетки у є Ж. Каждая клетка является соседом сама себе. Входное состояние элементарного автомата А в клетке у є Ж поставлено в соответствие внутренним состояниям соседей этой клетки. Таким образом, структура множества клеток Ж клеточного автомата представляется графом, в котором вершинами являются клетки, а ребра соответствуют отношению соседства. Этот граф

имеет регулярную структуру и степени вершин равные Ь. Состояние клетки у є Ж представлено вектором 5(у) с целочисленными компонентами 50Су), ^(у), ..., sЬ(w). Множество состояний 5(у) всех клеток м> є Ж в один и тот же момент времени ґ называется глобальным состоянием с(ґ) = {^(^1), 5^У2),..., 5(у), ...} клеточного автомата СА. На рис. 1 изображена клетка у, векторы скорости с находящихся в ней частиц и ее соседи МУ), і = 0, 1, ..., 6. Таким образом, количество соседей каждой клетки у модели БИР-МР равно семи, одним из соседей является сама клетка у, т.е. Мо(у) = У.

Рис. 1. Векторы скорости частиц. Нумерация соседей

Состояние клетки

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

В отличие от ЕИР, в модели БИР-МР компоненты 50У), 51(^), ..., &'6^) вектора состояния ж(^) клетки V принимают не булевы, а целочисленные значения. Таким образом, масса частиц в клетке V равна

ь

т(V) = ^ 5г (V), (1)

г=0

где Ь = 6 - количество возможных направлений вектора скорости, - г-й компо-

нент вектора состояний ж. Физическая интерпретация значений компонентов вектора ж(^) следующая: определяет количество частиц с единичной массой в клет-

ке V, векторы скорости с которых направлены в сторону соседа М(у).

Модельный импульс р в клетке V е Ж есть сумма всех импульсов рг = &'гсг, направленных к соседям М^), где г = 0, 1, ., Ь, а Ь = 6:

р=Ё

і=1

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

л/3

Рх =—(( + 5з - - 56); ()

Ру = 54 - 51 +1 (3 + 55 - 52 - 56 ), (4)

где • - сумма частиц в клетке V, с векторами скорости, направленными к соседу М-^).

Поведение клеточного автомата

Определим разбиение клеток V е Ж на типы. Клетками среды wср е Жср назовем клетки, в которых выполняются законы сохранения массы и импульса. Клетками стенок wст е Жст называются клетки, в которых выполняется закон сохране-

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

В модели БИР-МР используется клеточный автомат с синхронным режимом функционирования. На каждом такте происходит смена состояний «(О элементарных автоматов А во всех клетках V е Ж на состояния + 1) = 8(«(ф, где 8(«(ф -функция переходов элементарного автомата А. Клеточный автомат СА при этом переходит из глобального состояния с(0 в новое глобальное состояние + 1).

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

5(«) = 52(51(«)). (5)

Каждая из функций 51 и 52 должна удовлетворять законам сохранения массы:

X X 5У ( ^)) = ХХ^' И , У е {1, 2} (6)

wеW г=1 wеW г=1

и импульса:

XX 5] ( (^с1 (W)) = 'Х—1^г Мс (w) , У е {1, 2}. (7)

wеW г=1 wеW г=1

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

В фазе сдвига в каждой клетке V е Ж каждая частица, учтенная в компонентах

&’г^), при г = 1, ..., 6, вектора состояния «^), перемещается в соседнюю клетку

МДУ), соответствующую ее вектору скорости с. Частицы, учтенные в компоненте •0, остаются в клетке V. Таким образом, г-й компонент 5,^) вектора состояния «^) клетки V после сдвига принимает значение

Яг (%+2}тоё6)+1 М) , для г = I,2, •••, Ь;

sг (V), для г = 0.

(8)

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

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

Среда. В клетке среды wср е Жср функция 52 выбирается такой, чтобы сохранялись масса ш(у):

частиц в клетке. Значение функции 52 выбирается равновероятно из всевозможных состояний клетки wср, удовлетворяющих условиям (9) и (10). При выполнении этих условий условия (6) и (7) тем более выполняются.

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

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

Источники. Каждая клетка-источник wист е Жист поддерживает заданную концентрацию частиц и0^ист). Для этого она генерирует частицы со всевозможными направлениями вектора скорости в случае, если текущая концентрация частиц и^ист) < и0^ист). Количество генерируемых частиц равно разности и0^ист) -п^ист) заданной и текущей концентраций. Из клеток-источников можно создавать различные объекты. Например, установив их в пространстве в одну линию (как правило, у границы клеточного массива), можно получить источник равномерного потока частиц заданной концентрации. Отдельно установленный источник будет моделировать форсунку. Естественно, при генерации новых частиц ни масса ш^ист), ни импульс р^ист) не сохраняются. Граничные условия в клетках-источниках допускают нарушение условий (6) и (7).

Осредненные значения. При моделировании потоков практический интерес представляют не столько значение параметров автомата на микроуровне, т.е. масса ш(у) и скорость с^) частиц в каждой клетке V е Ж, сколько осредненные значения их скоростей (и) и концентраций (п) по некоторой окрестности Ау^), кото-

(9)

г=0

г=0

и импульс р^):

(10)

г=0

I=0

для г = 1,2,...,ь, для г = 0.

(11)

рая включает все клетки Wj е Ж, удаленные от клетки V не более чем на некоторую величину г, называемую радиусом осреднения.

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

И М= А 1 М X Х^ , ()

|Ау М| wjеА^) г =0

где ^-Av(w)| - количество клеток, попадающих в окрестность осреднения Ау^), сг -единичный вектор скорости, соответствующий г-му разряду вектора состояния s{wj), а • - значение г-го разряда вектора состояния клетки Wj е Ay(w).

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

1 Ь

(п)^)=А ( X X •. (3)

|Ау М1 wj еAv(w) г=0

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

Заметим, что осредненные значения модельных скорости и концентрации будут соответствовать их физическим аналогам только в том случае, когда окрестность осреднения Ау^) состоит исключительно из клеток среды V е Ж. Если по-

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

2. Экспериментальное исследование модели FHP-MP

Для проверки свойств предложенной модели была создана ее программная реализация, позволяющая проводить вычислительные эксперименты как на однопроцессорных компьютерах, так и на многопроцессорных и многомашинных вычислительных комплексах. Код программ написан на языке Си, параллелизм реализован с помощью библиотеки МР1. Ниже описаны численные эксперименты, проведенные с моделью БИР-МР. Получены качественные характеристики моделируемого потока.

Двумерная аппроксимация потока жидкости между двумя параллельными плоскостями

Этот эксперимент, ставший уже классическим, позволяет проверить модель на соответствие физике. Суть его в том, что продольная скорость потока (и)у (поток

движется вдоль Оу в положительном направлении) при скорости на границах (и) = 0 должна быть распределена вдоль направления Ох по параболическому

закону.

Клеточный автомат, использующийся в этом вычислительном эксперименте, имеет размеры 100 х 2000 клеток (вдоль декартовых осей Ох и Оу соответственно). Клетки с координатами в интервале [(2, 1), (99, 1)] являются источниками;

клетки с координатами в интервалах [(1, 1), (1, 2000)] и [(100, 1), (100, 2000)] -стенки. Остальные клетки - клетки среды. Такая двумерная конструкция является сечением параллелепипеда бесконечной ширины (вдоль оси От) и аппроксимирует трехмерный поток между двумя параллельными плоскостями.

На рис. 2 изображена проекция скорости потока (и) на ось Оу в поперечном

сечении потока (вдоль Ох). Кривая с маркерами построена исходя из результатов численного моделирования. Проекция аппроксимирована уравнением параболы (на рис. 2 - кривая без маркеров):

иу (г) = -0.0008г2 + 2, (14)

где г - расстояние от середины до рассматриваемой точки поперечного сечения.

Расстояние вдоль оси Ох, клетки

Рис. 2. Парабола Пуазейля и экспериментальная скорость потока

Кривая (14) - теоретически обоснованная парабола Пуазейля - одно из немногих аналитических решений уравнения Навье - Стокса, имеющее вид

Н (г) “(я2 -г2), ()

где ёР - падение давления на участке трубы длиной I, п - динамическая вязкость жидкости, Я - радиус трубы (в двумерном случае Я - расстояние между плоскостями).

Количество итераций, после которого проведено осреднение, Т = 20000. Радиус осреднения г = 15 (клеток). Согласно условию, налагаемому на (12) и (13), ос-редненные значения не могут быть получены на расстоянии до стенок ближе, чем г = 15 клеток. Следовательно, теоретическая и экспериментальная кривые в диапазоне, изображенном на рисунке, не опускаются до нуля.

Максимальная скорость потока в проведенном эксперименте, как видно на рис. 2, получилась равной двум, что несколько больше, чем в классической модели БИР. Но данная скорость не является предельной; в новой модели можно увеличивать концентрацию, вырабатываемую источниками, и скорость потока будет значительно превосходить полученную. Следовательно, новая модель позволяет моделировать турбулентные потоки.

Результаты этого эксперимента показывают, что модель БИР-МР адекватно отражает процессы в потоках и соответствует физике.

Поток с задвижкой

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

' ' ' * Ч ^ Ч ^ к / / / /7 ^ //

- t ? / / /

V \ \ \ * / / / / / ?

V ч ч , * і \ f У t ? / ( > ,

\ V ч ч V і ї f Ї / і / / 1 \

ч ж. ч \ Ч ч Ч \ 1 ! 1 * і \ ,

. » . . V _ _ ч V ч ч ч ч Ч \ * І Ч \

! і „ ч * - «*. ^ V V ч ч ч \ X Ч Ч \ ч t г - .

^ . V ^ ч V Ч ж. V ч \ ч \ ч ч V ч ч \ _

\ \ . . ч ч х Ч ^ Ч ч ■V У ч ч Ч Ч ч ч 1 ч V у _

\ \ \ N \ Ч V \ ч Ч Ч Ч ч ч ч ч

. * ч . , . . Ч * Н * Ч * ч ч ч * \ 1 Ч ч ч ч Ч ч ч V Ч

* ♦ Ч * ч ч * * * ч ч ч ч ч ч \ ч

, ч . . . . ♦ \ і * * \ * * * * \ ч , ♦ Ч Ч ч V \

, * \ і ч Ч і \ Ч Ч Ч \ \ Ч >ч

4 ■ ч ч * ч ч Ч Ч N 'і ! * \ ч ч ч ч ! ч Ч ч ч ч ■V ч ч

^ ^ ^ ... * і , ! ! \ 1 f ч ч ч X V ч ч ч,

* ь \ і 4 ‘ ' \ \ * ' * ' ' * ■ ’ 1 ' * 1 1 ! ч * ч 1 ч ' 4 - '

Рис. 3. Поле скорости потока с препятствием в виде задвижки.

Фрагмент за задвижкой

Обтекание круглого препятствия

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

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

Рис. 4. Поле скорости потока с круглым препятствием

Заключение

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

ЛИТЕРАТУРА

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

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

3. Rothman D.H., Zaleski S. Lattice-Gas Cellular Cutomata: Simple Models of Complex Hydro-dinamics. Cambridge University Press, 1997.

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