An extension of the cellular-automaton FHP-I flow model to the FHP-MP multiparticle model .pdf Одним из перспективных направлений имитационного моделирования физических процессов является моделирование клеточными автоматами. Клеточно-автоматные модели потоков, называемые решеточными газами (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, называемый элементарным автоматом.Внутренним состоянием автомата является целочисленный вектор, а не булев вектор как в классической модели FHP. Каждой клетке w е W сопоставлены некоторые координаты x(w) и y(w) на декартовой плоскости. Следовательно, между любыми двумя клетками wi е W и w2 е W можно подсчитать расстояние d(wi, w2).Для каждой клетки w е W определено некоторое упорядоченное множество N(w) = {N(w): N0(w) = w, N(w) е W & d(w, N(w)) = i, (i = i, 2, b)}, элементы которого находятся в отношении соседства с клеткой w и называются ее соседними клетками, или соседями. Константа b определяет количество нетождественных соседей каждой клетки w е W. Каждая клетка является соседом сама себе. Входное состояние элементарного автомата A в клетке w е W поставлено в соответствие внутренним состояниям соседей этой клетки. Таким образом, структура множества клеток W клеточного автомата представляется графом, в котором вершинами являются клетки, а ребра соответствуют отношению соседства. Этот графимеет регулярную структуру и степени вершин равные b. Состояние клетки w е W представлено вектором s(w) с целочисленными компонентами s0(w), si(w), sb(w). Множество состояний s(w) всех клеток w е W в один и тот же момент времени t называется глобальным состоянием o(t) = {s(wi), s(w2),..., s(wi), ...} клеточного автомата CA. На рис. i изображена клетка w, векторы скорости с находящихся в ней частиц и ее соседи Ni(w), i = 0, i, 6. Таким образом, количество соседей каждой клетки w модели FHP-MP равно семи, одним из соседей является сама клетка w, т. е. N0(w) = w.Состояние клеткиСостояние клетки w в каждый момент дискретного времени t однозначно определяется набором находящихся в ней частиц. Вектор скорости с каждой из них либо направлен в сторону одной из соседних клеток Ni(w) (при i = i, ., 6), либо равен нулю (при i = 0).В отличие от FHP, в модели FHP-MP компоненты s0(w), si(w), s6(w) вектора состояния s(w) клетки w принимают не булевы, а целочисленные значения. Таким образом, масса частиц в клетке w равнаbm( w) = Y s, (w),(iii=0где b = 6 - количество возможных направлений вектора скорости, s\ - i-й компонент вектора состояний s. Физическая интерпретация значений компонентов вектора s(w) следующая: s\ определяет количество частиц с единичной массой в клетке w, векторы скорости с которых направлены в сторону соседа Ni(w).Модельный импульс p в клетке w е W есть сумма всех импульсов pi = sici, направленных к соседям Ni(w), где i = 0, i, ., b, а b = 6:Используя (1) и рис. 1, достаточно просто подсчитать проекции px и py импульса p на декартовы оси Ox и Oy:л/3px =- (s2 + s3 - s5 - s6 );(3)py = s4 - S1 +1 (s3 + S5 - S2 - S6где Sj - сумма частиц в клетке w, с векторами скорости, направленными к соседуN(w).Поведение клеточного автоматаОпределим разбиение клеток w е W на типы. Клетками среды е Цгср назовем клетки, в которых выполняются законы сохранения массы и импульса. Клетками стенок wCT е WCT называются клетки, в которых выполняется закон сохранения массы, но может нарушаться закон сохранения импульса. И, наконец, источники е Wjjct суть клетки, в которых могут нарушаться как закон сохранения массы, так и закон сохранения импульса. Множества клеток среды W^,, стенок WCT и источников попарно не пересекаются (W^ n WCT = 0, n = 0, WCT n Wист = 0). Объединение этих множеств совпадает с множеством всех клеток автомата (W^ u WCT u = W). Поведение стенок и источников задает граничные условия клеточного автомата.В модели FHP-MP используется клеточный автомат с синхронным режимом функционирования. На каждом такте происходит смена состояний s(t) элементарных автоматов A во всех клетках w е W на состояния s(t + 1) = 5(s(t)), где 5(s(t)) -функция переходов элементарного автомата A. Клеточный автомат CA при этом переходит из глобального состояния о(?) в новое глобальное состояние о(? + 1).Каждый такт работы клеточного автомата выполняется в две фазы: сдвиг и столкновение. Функция переходов 5 элементарного автомата A состоит, таким образом, из композиции функций 51 (сдвиг) и 52 (столкновение):5(s) = 52(51(s)). (5)о££5ооКаждая из функций 1 и 2 должна удовлетворять законам сохранения массы:£ £5j (s (w))=£ Xs, (w) , j е {1, 2}, (6)и импульса:££5j (s, (w)Ci (w)) = ££s, (w)Ci (w) , j е {1, 2}. (7)wеW i=1wеW i=1С точки зрения динамики потока жидкости наличие этих двух фаз интерпретируется следующим образом. Столкновения реализуют диффузию в жидкости, а сдвиг - процесс переноса вещества в потоке. Далее подробно описаны эти две фазы.В фазе сдвига в каждой клетке w е W каждая частица, учтенная в компонентах si(w), при i = 1, 6, вектора состояния s(w), перемещается в соседнюю клетку Ni(w), соответствующую ее вектору скорости с;. Частицы, учтенные в компоненте s0, остаются в клетке w. Таким образом, i-й компонент si(w) вектора состояния s(w) клетки w после сдвига принимает значение 8 (s И) = {* (%+2)m0đ6)+1 '' = ''2'",* (8)U (w),для i = 0.Несмотря на то, что при сдвиге масса и импульс частиц в клетке изменяются, в пределах всего клеточного автомата они сохраняются, т.е. условия (6) и (7) выполняются.В фазе столкновения происходит изменение направления движения частиц согласно некоторым правилам столкновения, не зависящим от состояний соседних клеток, т.е. 52 зависит только от внутреннего состояния своего элементарного автомата. В модели FHP-MP функция 52 вероятностная. Ниже описаны правила столкновения для клеток разных типов: среды, стенок и источников.Среда. В клетке среды >еср е W^, функция 52 выбирается такой, чтобы сохранялись масса m(w):X 82 (s, ^ )) = £s ^ ) (9)i=0i=0и импульс p(w):X 82 (Ci (>"ср )) = X Si (Wcp ) Ci Кр ) (10)i=0 l=0частиц в клетке. Значение функции 52 выбирается равновероятно из всевозможных состояний клетки w^, удовлетворяющих условиям (9) и (10). При выполнении этих условий условия (6) и (7) тем более выполняются.Стенки. В клетках >ест е WCT, являющихся стенками, частицы «отражаются» в обратном направлении, нарушая при этом закон сохранения импульса.й ( ()) [S((i+2)mod6)+1 (Wcm ), для i = 1, 2,b,Из-за того, что количество частиц в клетке не меняется, условия (9), а следовательно, и (6) выполняются. Условие (7) может нарушаться, так как меняются направления векторов скорости c частиц, но это допускается граничными условиями. Такое поведение частиц в клетках-стенках моделирует условие нулевой скорости потока на границах препятствий.Источники. Каждая клетка-источник >еист е поддерживает заданную концентрацию частиц w^Wjict). Для этого она генерирует частицы со всевозможными направлениями вектора скорости в случае, если текущая концентрация частиц «(ww) < «0^ист). Количество генерируемых частиц равно разности «0^ист) -и(>еист) заданной и текущей концентраций. Из клеток-источников можно создавать различные объекты. Например, установив их в пространстве в одну линию (как правило, у границы клеточного массива), можно получить источник равномерного потока частиц заданной концентрации. Отдельно установленный источник будет моделировать форсунку. Естественно, при генерации новых частиц ни масса т^ист), ни импульс р^ист) не сохраняются. Граничные условия в клетках-источниках допускают нарушение условий (6) и (7).Осредненные значения. При моделировании потоков практический интерес представляют не столько значение параметров автомата на микроуровне, т.е. масса m(w) и скорость ci(w) частиц в каждой клетке w е W, сколько осредненные значения их скоростей (и) и концентраций («) по некоторой окрестности Av(w), которая включает все клетки Wj е W, удаленные от клетки w не более чем на некоторую величину r, называемую радиусом осреднения.Осредненная скорость вычисляется как сумма всех векторов частиц, попадающих в окрестность осреднения Av(w), деленная на мощность окрестности осреднения:(w)= TA4^ its .(13)lAv (w) wjеAv(w) i=0Осредненные значения скорости и концентрация частиц, являющиеся модельными значениями скорости и давления, соответствуют значениям скорости и давления моделируемой жидкости и являются параметрами макроуровня.Заметим, что осредненные значения модельных скорости и концентрации будут соответствовать их физическим аналогам только в том случае, когда окрестность осреднения Av(w) состоит исключительно из клеток среды w е W. Если последнее не выполняется для окрестности осреднения, то будем полагать значения (и) и (п) неопределенными. Это условие не позволяет определять значения (и) и (п) на расстоянии ближе, чем радиус осреднения r от границ (стенок и источников), в том числе и находящихся внутри моделируемого объекта.2. Экспериментальное исследование модели FHP-MPДля проверки свойств предложенной модели была создана ее программная реализация, позволяющая проводить вычислительные эксперименты как на однопроцессорных компьютерах, так и на многопроцессорных и многомашинных вычислительных комплексах. Код программ написан на языке Си, параллелизм реализован с помощью библиотеки MPI. Ниже описаны численные эксперименты, проведенные с моделью FHP-MP. Получены качественные характеристики моделируемого потока.Двумерная аппроксимация потока жидкости между двумя параллельными плоскостямиЭтот эксперимент, ставший уже классическим, позволяет проверить модель на соответствие физике. Суть его в том, что продольная скорость потока {u)y (потокдвижется вдоль Oy в положительном направлении) при скорости на границах (и) = 0 должна быть распределена вдоль направления Ox по параболическомузакону.Клеточный автомат, использующийся в этом вычислительном эксперименте, имеет размеры 100 х 2000 клеток (вдоль декартовых осей Ox и Oy соответственно). Клетки с координатами в интервале [(2, 1), (99, 1)] являются источниками;клетки с координатами в интервалах [(1, 1), (1, 2000)] и [(100, 1), (100, 2000)] -стенки. Остальные клетки - клетки среды. Такая двумерная конструкция является сечением параллелепипеда бесконечной ширины (вдоль оси Oz) и аппроксимирует трехмерный поток между двумя параллельными плоскостями.На рис. 2 изображена проекция скорости потока (и) на ось Oy в поперечномсечении потока (вдоль Ox). Кривая с маркерами построена исходя из результатов численного моделирования. Проекция аппроксимирована уравнением параболы (на рис. 2 - кривая без маркеров):(и)у (r) = -0.0008r2 + 2,(14) где r - расстояние от середины до рассматриваемой точки поперечного сечения.2 5 -, Осредненная скорость потока (и )y , модельные единицы0,5 --30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 Расстояние вдоль оси Ох, клетки Рис. 2. Парабола Пуазейля и экспериментальная скорость потокаКривая (14) - теоретически обоснованная парабола Пуазейля - одно из немногих аналитических решений уравнения Навье - Стокса, имеющее видИ(r ) =^ ((2 - r 2 ),(155где dP - падение давления на участке трубы длиной /, ц - динамическая вязкостьжидкости, R - радиус трубы (в двумерном случае R - расстояние между плоскостями).Количество итераций, после которого проведено осреднение, T = 20000. Радиус осреднения r = 15 (клеток). Согласно условию, налагаемому на (12) и (13), ос-редненные значения не могут быть получены на расстоянии до стенок ближе, чем r = 15 клеток. Следовательно, теоретическая и экспериментальная кривые в диапазоне, изображенном на рисунке, не опускаются до нуля.Максимальная скорость потока в проведенном эксперименте, как видно на рис. 2, получилась равной двум, что несколько больше, чем в классической модели FHP. Но данная скорость не является предельной; в новой модели можно увеличивать концентрацию, вырабатываемую источниками, и скорость потока будет значительно превосходить полученную. Следовательно, новая модель позволяет моделировать турбулентные потоки.Результаты этого эксперимента показывают, что модель FHP-MP адекватно отражает процессы в потоках и соответствует физике.Поток с задвижкойСледующий вычислительный эксперимент был проведен с целью исследования обтекания препятствий. Для этого в клеточный автомат, использовавшийся в предыдущем эксперименте, были добавлены граничные условия в виде препятствия, имеющего форму задвижки. Задвижка была установлена на расстоянии трети длины трубы от источников перпендикулярно направлению распространения потока и перекрывала трубу наполовину. Получившееся поле скорости похоже на обтекание задвижки потоком. К сожалению, при мелком масштабе изображения, охватывающем весь клеточный автомат, видны только высокоскоростные части потока. Низкие скорости, каковыми являются скорости за задвижкой, приходится рассматривать отдельно. На рис. 3 в более крупном масштабе приведен фрагмент автомата, находящийся непосредственно за задвижкой. Направление каждой стрелки на рисунке совпадает с направлением потока в ее основании, а длина пропорциональна скорости потока в соответствующей точке. Поток направлен слева направо, задвижка изображена в левой части рисунка. На этом фрагменте отчетливо видны завихрения потока, что указывает на его турбулентность.Обтекание круглого препятствияЕще один эксперимент - с круглым препятствием - также показывает турбулентность потока. В исходный клеточный автомат было добавлено препятствие, имеющее форму круга. Поле скорости этого автомата изображено на рис. 4.За круглым препятствием видны оторвавшийся вихрь и «змейка», очевидно являющаяся вариантом дорожки Кармана для случая, когда расстояния до ограничивающих поток плоскостей слишком мало.ЗаключениеВ статье предложена новая клеточно-автоматная модель потока FHP-MP, в которой допускается более одной частицы в клетке с равными векторами скорости. Вычислительные эксперименты, проведенные с этой моделью, показывают ее корреляцию с классической моделью FHP и с физическими законами. Экспериментально полученное распределение скорости в модели FHP-MP согласуется с параболой Пуазейля. Максимальная скорость потока в проведенном эксперименте получилась несколько больше, чем в модели FHP. Кроме того, в новой модели можно увеличивать начальную концентрацию частиц, получая скорость потока, значительно превосходящую скорость в классической модели. Следовательно, модель FHP-MP позволяет моделировать турбулентные потоки, что и подтвердили эксперименты с задвижкой и круглым препятствием. В результате этих экспериментов получены потоки с завихрениями, в том числе наблюдалась дорожка Кармана за препятствием, что является признаком турбулентного потока.

Скачать электронную версию публикации