МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное бюджетное образовательное учреждение
высшего профессионального образования
«КУБАНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ»
(ФГБОУ ВПО «КубГУ»)
Кафедра математического моделирования
ДОПУСТИТЬ К ЗАЩИТЕ В ГАК
Заведующий кафедрой академик РАН,
д-р физ.-мат. наук, профессор
_____________ Бабешко В.А.
_________________2013 г.
ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ (ДИПЛОМНАЯ)
РАБОТА
МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ТЕЧЕНИЯ ЖИДКОСТИ
Работу выполнил И.И. Шкурат
Факультет ФКТиПМ
Специальность 010501 — Прикладная математика и информатика
Научный руководитель,
канд. физ.-мат наук, доцент С.Е. Рубцов
Нормоконтролер доцент, канд. физ.-мат наук М.С. Капустин
Краснодар 2013
РЕФЕРАТ
Дипломная работа 45 с., 29 рис., 5 источников, 2 приложения.
ПРОСТРАНСТВЕННАЯ ДИНАМИКА, КЛЕТОЧНЫЙ АВТОМАТ, МОДЕЛЬ ПОТОКА.
Объектом исследования являются клеточные автоматы моделирующие течение потоков жидкости.
Цель работы:
— Исследование клеточно-автоматных моделей газовой динамики с помощью клеточных автоматов;
— разработка программы, реализующей клеточно-автоматную модель потока жидкости FHP-I;
— модифицировать модель FHP-I, добавлением частиц с новыми свойствами;
— экспериментально проверить модель FHP-I на соответствие физике.
В среде программирования Borland C++ Builder 6 разработано приложение, реализующее клеточно-автоматную модель потока жидкости FHP-I, а также модифицированную модель.
СОДЕРЖАНИЕ
ВВЕДЕНИЕ
1 Клеточно-автоматная пространственная динамика
2 Клеточно-автоматная модель потока жидкости FHP-I
2.1 Основные определения
2.2 Поведение клеточного автомата
3 Клеточно-автоматная модель потока жидкости FHP-I с примесью
4 Описание программной реализации КА модели FHP-I
5 Экспериментальное исследование модели FHP-I
5.1 Эксперимент 1. Двумерная аппроксимация потока жидкости между двумя параллельными плоскостями
5.2 Эксперимент 2. Поток с задвижкой
5.3 Эксперимент 3. Распространение частиц примеси
ЗАКЛЮЧЕНИЕ
СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
ПРИЛОЖЕНИЕ А
ПРИЛОЖЕНИЕ Б
ВВЕДЕНИЕ
Клеточный автомат (КА) был предложен фон-Нейманом в середине прошлого века [1]. На основе этой модели подтверждалась мысль о том, что человек может создать устройство, обладающее присущими ему самому свойствами, в частности, способностью воспроизводить себе подобного. Поэтому КА долгое время рассматривался как модель самовоспроизведения и толковался как упрощенная модель некоторого биологического сообщества, состоящего из множества клеток. Каждой клетке ставится в соответствие конечный автомат, называемый элементарным автоматом, который может находиться в одном из двух состояний: 0 и 1 и изменять это состояние или в зависимости от состояний клеток некоторого своего окружения, называемого соседством. Алгоритм вычисления следующего состояния в зависимости от состояний соседей у всех клеток одинаковый. Все клетки выполняют переход в новое состояние одновременно, т.е. синхронно и параллельно. При этом наблюдается изменение глобальной черно-белой картины распределения состояний по пространству КА. Такая картина называется конфигурацией КА.
Итеративная смена конфигурации при переходах всех элементарных состояний в новые состояния называется эволюцией КА. Эволюционируя, КА моделирует пространственную динамику, которая может иметь завершение, повторяться периодически или изменяться хаотично. Исследования этой «модели мира» показали, что, несмотря на простоту каждой клетки, их кооперативная работа моделирует очень сложные и разнообразные процессы, которые иногда невозможно (или, по крайней мере, неизвестно как) описать другим способом. В таком понимании КА теперь называется классическим КА.
В середине 80-х годов прошлого века произошел всплеск интереса к КА, в связи с построением новых моделей вычислений. Были предложены КА, эволюция которых моделирует процессы диффузии, разделения фаз, реакционно-диффузионные процессы, знаменитую реакцию Белоусова-Жаботинского, движения солитонов, образование диссипативных структур. Переход от булевых пространств к привычным непрерывным пространственным функциям (макровеличинам) производится путем осреднения значений состояний клеток по некоторой заданной окрестности.
Настоящим потрясением основ в моделировании пространственной динамики стало появление клеточно-автоматной газовой динамики, названной «Gas-Lattice», что в переводе на русский звучит как «решеточный газ». Очень важен тот факт, что было строго математически доказано соответствие Gas-Lattice моделей уравнению Навье-Стокса [2].
Предполагается, что этот новый подход к моделированию природных явлений станет дополнением к традиционным моделям математической физики, основанным на дифференциальных уравнениях в частных производных и хорошо развитым численным методам их решения на компьютерах, вычислительная мощность и архитектура которых меняется быстрее, чем создается необходимое математическое обеспечение.
Описание использования клеточных автоматов для моделирования различных физических явлений приведены в статье О.Л. Бандман [6].
1. Клеточно-автоматная пространственная динамика
Мотивацией к развитию и распространению КА-моделей послужили два фактора: технический и математический. Технический фактор, в свою очередь имеет два аспекта. Во-первых, было подсчитано что КА-алгоритмы, воплощенные в «железо» могут увеличить производительность на 3 порядка и выше по отношению к решению аналогичных задач традиционными методами, например, решению дифференциальных уравнений в частных производных. Для подтверждения этого тезиса была построена КА-машина (САМ-8), на которой производились испытания КА-моделей, в частности, испытывались алгоритмы решеточной гидро- и газодинамики. Второй технический фактор состоит в том, что появились высокопроизводительные многопроцессорные суперкомпьютеры и сетевые вычислительные технологии, которые позволили моделировать эволюцию КА очень больших размеров достаточно быстро, ведь КА-алгоритмы чрезвычайно легко распараллеливаются на любое число процессоров из-за присущего им естественного параллелизма. И, хотя в этих случаях больших ускорений вычислений не наблюдается, КА-моделирование на многопроцессорных системах выдерживает конкуренцию с традиционным моделированием.
Второй фактор, вызвавший к жизни КА-моделирование, составляют трудности численного моделирования на основе решения дифференциальных уравнений, возникшие в связи с ростом размеров и сложности решаемых задач, а также необходимостью их размещения на многих процессорах. Основная проблема состоит в том, что решение стационарных задач математической физики обычно выполняется неявными методами, которые не накладывают ограничений на вычислительную устойчивость, что делает их более быстрыми, чем явные. Однако параллельная реализация этих методов требует дополнительных усилий, как математиков, так и программистов, и всегда связана с деградацией эффективности распараллеливания при увеличении количества процессоров. При решении нестационарных задач предпочтение отдается явным методам, которые легко и эффективно распараллеливаются, но ограничены условиями вычислительной устойчивости (требуют применения малого шага по времени). Явные методы обладают многими свойствами КА, иногда рассматриваются как «непрерывные клеточные автоматы».
Сравнение КА-моделирования с традиционными методами численного моделирования может быть сделано сколько-нибудь корректно только для конкретных задач [3]. В тех случаях, для которых известны дифференциальные уравнения, но неизвестна их КА-модель или наоборот, КА-модель известна, а дифференциальное уравнение не существует, выбора нет, и сравнение не имеет смысла. Более того, для процессов с сильно выраженной нелинейностью, численные методы могут оказаться слишком сложными, и КА предпочтительнее. И, наоборот, для известных хорошо проработанных численных методов, а иногда и имеющегося программного обеспечения, предпочтение отдается численному методу.
Основные достоинства КА-моделей:
Отсутствие ошибок округления. Поскольку состояния клеток представлены булевыми векторами, и вычисления ведутся над ними, все разряды в них равноправны, округления не производится. При итеративных вычислениях это бывает важно.
Экономия памяти при больших объемах данных, так как булевы данные размещаются во всех разрядах машинных ячеек.
Виртуальная неограниченность параллелизма. Допускается разрезание области моделирования любым образом. Распараллеливание на любое количество процессоров не приводит к деградации эффективности распараллеливания.
Простота задания граничных условий. Границы областей моделирования представлены особыми клетками, которые имеют граничные правила переходов. Даже при очень сложных формах границ (пористые среды, изогнутые трубки) алгоритм моделирования практически не усложняется по сравнению с тем, когда препятствий нет.
Несомненно, КА-модели и основанные на них алгоритмы имеют и недостатки, с которыми приходится бороться. Главными из них являются следующие:
Автоматный шум, связанный с дискретным представлением данных. Избежать его невозможно, но при достаточно большом радиусе осреднения, он становится незаметным.
Отсутствие формальных методов «предсказания» поведения КА, по заданным функциям перехода. Отсюда отсутствие методов синтеза функции перехода КА по заданным свойствам его эволюции.
2 Клеточно-автоматная модель потока жидкости FHP-Й
2.1 Основные определения
Под клеточным автоматом КА модели FHP-Й будем понимать тройку объектов (W, A, N), где — множество клеток, заданное их координатами в некотором дискретном пространстве. Каждой клетке поставлен в соответствие конечный автомат A, называемый элементарным автоматом [7].
Внутренним состоянием автомата является булев вектор. Каждой клетке w W сопоставлены некоторые координаты и на декартовой плоскости. Следовательно, между любыми двумя клетками и можно подсчитать расстояние .
Для каждой клетки w W определено некоторое упорядоченное множество
,
элементы которого находятся в отношении соседства с клеткой w и называются ее соседними клетками, или соседями. Константа b определяет количество нетождественных соседей каждой клетки . Каждая клетка является соседом сама себе. Входное состояние элементарного автомата A в клетке поставлено в соответствие внутренним состояниям соседей этой клетки. Таким образом, структура множества клеток W клеточного автомата представляется графом, в котором вершинами являются клетки, а ребра соответствуют отношению соседства. Этот граф имеет регулярную структуру и степени вершин равные b. Состояние клетки представлено вектором с булевыми компонентами . Множество состояний всех клеток в один и тот же момент времени t называется глобальным состоянием клеточного автомата КА. На рисунке 1 изображена клетка w, векторы скорости находящихся в ней частиц и ее соседи .
Рисунок 1 ? Векторы скорости частиц. Нумерация соседей
Таким образом, количество соседей каждой клетки w модели FHP-Й равно семи, одним из соседей является сама клетка w, т. е..
Состояние клетки w в каждый момент дискретного времени t однозначно определяется набором находящихся в ней частиц. Вектор скорости каждой из них либо направлен в сторону одной из соседних клеток (при ), либо равен нулю (при ).
Таким образом, масса частиц в клетке w равна
, (1)
где b = 6 — количество возможных направлений вектора скорости, — i-й компонент вектора состояний s. Физическая интерпретация значений компонентов вектора следующая: определяет частицу с единичной массой в клетке w, вектор скорости которой направлен в сторону соседа . На рисунке 2 изображен вектор состояния клетки.
Рисунок 2 ? Вектор состояния клетки
а) три двтжущихся частицы, нумерация разрядов
б) две двтжущихся и одна частица покоя
в) две двтжущихся частицы
Модельный импульс p в клетке w W есть сумма всех импульсов , направленных к соседям где , а :
(2)
Используя (1) и рисунок 1, достаточно просто подсчитать проекции и импульса p на декартовы оси Ox и Oy:
(3)
, (4)
где , — частица в клетке w, с вектором скорости, направленным к соседу.
2.2 Поведение клеточного автомата
Определим разбиение клеток на типы. Клетками среды назовем клетки, в которых выполняются законы сохранения массы и импульса. Клетками стенок называются клетки, в которых выполняется закон сохранения массы, но может нарушаться закон сохранения импульса. И, наконец, источники клетки, в которых могут нарушаться как закон сохранения массы, так и закон сохранения импульса. Множества клеток среды , стенок и источников попарно не пересекаются
(=?,=?,=?).
Объединение этих множеств совпадает с множеством всех клеток автомата (=). Поведение стенок и источников задает граничные условия клеточного автомата.
В модели FHP-Й используется клеточный автомат с синхронным режимом функционирования. На каждом такте происходит смена состояний элементарных автоматов A во всех клетках на состояния , где — функция переходов элементарного автомата A. Клеточный автомат КА при этом переходит из глобального состояния в новое глобальное состояние . Каждый такт работы клеточного автомата выполняется в две фазы: сдвиг и столкновение. Функция переходов элементарного автомата A состоит, таким образом, из композиции функций (сдвиг) и (столкновение):
(5)
Каждая из функций и должна удовлетворять законам сохранения массы:
, (6)
и импульса:
(7)
С точки зрения динамики потока жидкости наличие этих двух фаз интерпретируется следующим образом. Столкновения реализуют диффузию в жидкости, а сдвиг — процесс переноса вещества в потоке. Далее подробно описаны эти две фазы.
В фазе сдвига в каждой клетке частица, учтенная в компонентах, при , вектора состояния , перемещается в соседнюю клетку , соответствующую ее вектору скорости . Частица, учтенная в компоненте , остается в клетке w. Таким образом, i-й компонент вектора состоянияклетки w после сдвига принимает значение
(8)
Несмотря на то, что при сдвиге масса и импульс частиц в клетке изменяются, в пределах всего клеточного автомата они сохраняются, т.е. условия (6) и (7) выполняются.
В фазе столкновения происходит изменение направления движения частиц согласно некоторым правилам столкновения, не зависящим от состояний соседних клеток, т.е. зависит только от внутреннего состояния своего элементарного автомата. В модели FHP-Й функция вероятностная. Ниже описаны правила столкновения для клеток разных типов: среды, стенок и источников.
Среда. В клетке среды функция выбирается такой, чтобы сохранялись масса :
(9)
и импульс :
(10)
частиц в клетке. Значение функции выбирается равновероятно из всевозможных состояний клетки , удовлетворяющих условиям (9) и (10). При выполнении этих условий условия (6) и (7) тем более выполняются. На рисунке 3 (а) показаны вероятностные правила столкновений, а на рисунке 3 (б) детерминированные.
Рисунок 3 ? Нетривиальные правила столкновений
а) вероятностные; б)детерминированные
На рисунках 4, 5, 6 показан фрагмент автомата FHP модели и его состояний после сдвига и после столкновений соответственно.
Рисунок 4 ? Фрагмент КА
Рисунок 5 ? Фрагмент КА. Результат применения правил сдвига к КА на предыдущем рисунке
Рисунок 6 ? Фрагмент КА. Результат применения правил столкновения к КА на предыдущем рисунке
Стенки. В клетках, являющихся стенками, частицы «отражаются» в обратном направлении, нарушая при этом закон сохранения импульса.
(11)
Из-за того, что количество частиц в клетке не меняется, условия (9), а следовательно, и (6) выполняются. Условие (7) может нарушаться, так как меняются направления векторов скорости частиц, но это допускается граничными условиями. Такое поведение частиц в клетках-стенках моделирует условие нулевой скорости потока на границах препятствий.
Источники. Каждая клетка-источник генерирует частицы со всевозможными направлениями вектора скорости. Из клеток-источников можно создавать различные объекты. Например, установив их в пространстве в одну линию (как правило, у границы клеточного массива), можно получить источник равномерного потока частиц. Отдельно установленный источник будет моделировать форсунку. Естественно, при генерации новых частиц ни масса , ни импульс не сохраняются. Граничные условия в клетках-источниках допускают нарушение условий (6) и (7).
Осредненные значения.
При моделировании потоков практический интерес представляют не столько значение параметров автомата на микроуровне, т.е. масса и скорость частиц в каждой клетке , сколько осредненные значения их скоростей и концентраций по некоторой окрестности , которая включает все клетки , удаленные от клетки не более чем на некоторую величину , называемую радиусом осреднения.
Осредненная скорость вычисляется как сумма всех векторов частиц, попадающих в окрестность осреднения , деленная на мощность окрестности осреднения:
, (12)
где — количество клеток, попадающих в окрестность осреднения , -единичный вектор скорости, соответствующий i-му разряду вектора состояния, а — значение i-го разряда вектора состояния клетки . На рисунке 7. изображен осредненный вектор скорости, жирными линиями выделены клетки, попадающие в окрестность осреднения. Вектор осредненной скорости изображен в центре окрестности.
Рисунок 7 ? Осредненный вектор скорости ( радиус осреднения r=2 клетки)
Осредненная концентрация частиц подсчитывается в той же окрестности следующим образом:
(13)
Осредненные значения скорости и концентрация частиц, являющиеся модельными значениями скорости и давления, соответствуют значениям скорости и давления моделируемой жидкости и являются параметрами макроуровня.
3. Клеточно-автоматная модель потока жидкости FHP-Й с примесью
Для моделирования потока жидкости с примесью каждая частица должна иметь свойство, по которому можно определить является ли эта частица примесью либо это «чистая» частица. Для этого вводится новое понятие состояние частицы, и объект модели называемый частицей. Тогда под клеточным автоматом КА модели FHP-Й будем понимать четверку объектов (W, A, N, O). — множество клеток, заданное их координатами в некотором дискретном пространстве, A конечный автомат каждой клетки , некоторое упорядоченное множество для всех клеток . множество всех сгенерированных частиц. Состояние частицы представлено вектором с булевыми компонентами . Где компонента определяет, является ли эта частица действующей, а компонента определяет, является ли эта частица примесью. На рисунках 8, 9, 10 изображены клетки с векторами состояния, где красный цвет стрелки означает частицы являющимися примесью, а голубой “чистые” частицы
Рисунок 8 — Клетка с вектором состояния(2 частицы приме и 1”чистая”)
Рисунок 9 — Клетка с вектором состояния (2 частицы).
Рисунок 10 — Клетка с вектором состояния (1 “чистая” частица покоя , 1 ”чистая” частица, 1 частица примеси).
В данной модели масса частицы являющейся примесью равна массе чистой частицы. Поэтому масса частиц в клетке w не изменяется и равна
(14)
Фаза сдвига элементарного автомата A остается прежней. В фазе столкновения происходят изменения. В прежней модели зависела от внутреннего состояния своего элементарного автомата, а именно от количества и расположения частиц в клетке. Так как в данной модели частицы имеют свои состояния, то теперь зависит также и от состояния каждой частицы. За счет этого правила столкновения усложняются, а частицы получают большую свободу движения. Фактически детерминированные правила столкновения модели FHP-I без примеси трансформировались в вероятностные, в тех случаях, когда в клетке присутствуют оба вида частиц. На рисунке 12 показано, как изменились правила столкновения с вводом в модель особых частиц, а на рисунке 11 показаны детерминированные.
Рисунок 11 — Детерминированные правила столкновений.
Рисунок 12 — Вероятностные правила столкновений.
Каждая из функций и по-прежнему удовлетворяет законам сохранения массы:
, (15)
и импульса:
(16)
В данной модели предполагается, что частицы примеси имеют массу отличную от массы чистых частиц, а также двигаются, учитывая закон всемирного тяготения. Поэтому, так как у частиц появляется новое свойство, вследствие которого они притягиваются к Земле, то к элементарному автомату A добавляется еще одна фаза — фаза воздействия силы тяжести.
В данной фазе происходит перемещение частиц примеси от верхних клеток к нижним. Эти перемещения можно рассмотреть на примере трех клеток. Одна клетка располагается вверху, и она может только отдавать частицы нижним клеткам, будем называть ее клетка-источник. Соответственно две другие клетки располагаются ниже, и они могут только принимать частицы, их будем называть клетками-приемниками.
Сначала рассматривается случай, когда в источнике есть частицы примеси, а в хотя бы одном из приемников присутствует частица, с тем же направлением движения, но которая не является частицей примеси.
Все частицы примеси, расположенные в источнике, перемещаются в клетку-приемник только в том случае, если в приемнике есть частица, движущаяся в том же направлении, что и частица в источнике, но эта частица не является частицей примеси. То в какую именно клетку-приемник попадет частица, зависит от направления движения частицы и от вектора состояния клетки-приемника. Из двух клеток-приемников выбирается та, расположение которой наиболее соответствует направлению движения частицы. Если же в данной клетки уже есть частица примеси, которая движется в том же направлении, тогда рассматривается возможность переместить частицу в соседнюю клетку-приемник. Если в соседней клетки существует частица с тем же направлением движения, но не являющаяся частицей примеси, тогда частица перемещается в данную клетку. Направление движения частицы, которая перемещается из источника в приемник остается прежним. Эти правила перемещения изображены на рисунке 13 . Затем рассматривается случай, когда в источнике есть частицы примеси, а в обоих приемниках вообще отсутствует частица, с тем же направлением движения. Критерии выбора клетки-приемника такие же, как и в первом случае. Эти правила изображены на рисунке 14. Для правильного перемещения частиц, эти правила применяются к клеткам, начиная от тех, которые располагаются в самом низу автомата, и затем двигаясь к клеткам находящимся выше. Это позволяет клетке сначала отдать, а затем быть способной принять частицы.
Рисунок 13 — Правила перемещения частиц при наличии у клетки-приемника частицы
Рисунок 14 — Правила перемещения частиц при отсутствии у клетки-приемника частицы
Применяя эти правила, в клетке нарушается закон сохранения импульсa :
,
но в целом в автомате он не нарушается, то есть условие (7):
по-прежнему выполняется.
4. Описание программной реализации КА модели FHP-Й
Для удобства в проведении экспериментов и получения результатов, которые легко анализировать, было написано Windows — приложение СА_FHP-Й Simulator. Программа написана на языке C++ в среде программирования Borland C++ Builder 6 и реализует модель FHP-I.
Из-за большой вычислительной сложности, основной код программы, в котором осуществляется фаза сдвига, фаза столкновения и генерация частиц, является самостоятельной функцией, т.е. внутри нее не вызываются другие процедуры и функции. Основной код программы представлен в приложении Б. Так как клетки в модели FHP-I имеют гексагональную структуру, то элементы двумерного клеточного массива, в которых располагаются клетки, имеют индексы такие, чтобы выполнялось равенство (i+j) mod 2=0. То есть массив можно представить в виде шахматной доски, и только элементы, которые находятся на белой стороне, содержат рабочие клетки. Поэтому все действия с клеточным массивом разбиты на две части, сначала обрабатываются клетки расположенные на нечетных индексах по строкам и столбцам, а затем расположенные на четных индексах. Нетривиальные правила столкновений были вычислены и реализуется в программе в виде макросов.
Для реализации модели потока жидкости FHP-I с примесью в программа СА_FHP-Й Simulator была доработана. Так как частицы теперь представляют собой не наличие нулей и единиц в клетке, а имеют свои состояния, была добавлена структура частицы с элементами, отвечающими за наличие сомой частицы и признака того, что эта частица является примесью. Вследствие этого были переписаны блоки программы выполняющие фазу столкновения и сдвига таким образом, чтобы они могли взаимодействовать со структурой частицы. Так как правила столкновения частиц получили большее количество вариантов код программы реализующий фазу столкновения также был переработан. В приложение была добавлена возможность визуализации потока уже не в виде векторов скорости, а в виде разноцветных участков демонстрирующих количество примеси в определенной области.
В центре приложения расположено главное окно, справа интерфейс пользователя. В главном окне можно в ручную задавать различные граничные условия в виде препятствий для потока жикости, что очень удобно. В этом же окне отображается состояние автомата в виде поля скорости потока. Получаемые векторы скорости можно настраивать для лучшей визуализации. Для того чтобы можно было рассмотреть участки как с высокими скоростями потока так и с низкими, есть возможность изменения масштаба получаемой картинки поля скорости. Эту картинку также можно сохранить на жесткий диск в формате .bmp
Рисунок 15 ? Приложение CA_FHP-I Simulator
газовый клеточный программа жидкость
5. Экспериментальное исследование модели FHP-Й
Клеточный автомат, использующийся в вычислительных экспериментах, имеет размеры 100 2000 клеток (вдоль декартовых осей Ox и Oy соответственно). Клетки с координатами в интервале [(2, 1), (99, 1)] являются источниками; клетки с координатами в интервалах [(1, 1), (1, 2000)] и [(100, 1), (100, 2000)] -стенки. Остальные клетки — клетки среды. Такая двумерная конструкция является сечением параллелепипеда бесконечной ширины (вдоль оси Oz) и аппроксимирует трехмерный поток между двумя параллельными плоскостями.
5.1 Эксперимент 1. Двумерная аппроксимация потока жидкости между двумя параллельными плоскостями
Этот эксперимент позволяет проверить модель на соответствие физике. Суть его в том, что продольная скорость потока (поток движется вдоль Oy в положительном направлении) при скорости на границах = 0 должна быть распределена вдоль направления Ox по параболическому закону.
На рисунке 16 изображена проекция скорости потока на ось Oy в поперечном сечении потока (вдоль Ox). Кривая с маркерами построена исходя из результатов численного моделирования. Проекция аппроксимирована уравнением параболы (на рисунке 16 — кривая без маркеров):
где r — расстояние от середины до рассматриваемой точки поперечного сечения.
Кривая (14) — теоретически обоснованная парабола Пуазейля — одно из немногих аналитических решений уравнения Навье — Стокса [3], имеющее вид
,
где dP — падение давления на участке трубы длиной l, з — динамическая вязкость жидкости, R — радиус трубы (в двумерном случае R — расстояние между плоскостями).
Количество итераций, после которого проведено осреднение, T = 20000. Радиус осреднения r = 15 (клеток). Согласно условию осредненные значения не могут быть получены на расстоянии до стенок ближе, чем r = 15 клеток. Следовательно, теоретическая и экспериментальная кривые в диапазоне, изображенном на рисунке, не опускаются до нуля.
Максимальная скорость потока в проведенном эксперименте, как видно на рисунке 16 получилась равной 1.5. Результаты этого эксперимента показывают, что модель FHP-Й адекватно отражает процессы в потоках и соответствует физическим явлениям.
Рисунок 16 ? Парабола Пуазейля и экспериментальная скорость потока
5.2 Эксперимент 2. Поток с задвижкой
Данный вычислительный эксперимент был проведен с целью исследования обтекания препятствий. Для этого в клеточный автомат были добавлены граничные условия в виде препятствия, имеющего форму задвижки. Задвижка была установлена на расстоянии трети длины трубы от источников перпендикулярно направлению распространения потока и перекрывала трубу наполовину. Получившееся поле скорости похоже на обтекание задвижки потоком. На рисунке 17 изображен весь клеточный автомат. На рисунке 18 в более крупном масштабе приведен фрагмент автомата, находящийся непосредственно за задвижкой. Направление каждой стрелки на рисунке совпадает с направлением потока в ее основании, а длина пропорциональна скорости потока в соответствующей точке. Поток направлен слева направо, задвижка изображена в левой части рисунка. На этом фрагменте отчетливо видны завихрения потока, что указывает на его турбулентность.
Рисунок 17 ? Поле скорости потока с препятствием в виде задвижки
Рисунок 18 — Поле скорости потока с препятствием в виде задвижки. Фрагмент за задвижкой
5.3 Эксперимент 3. Распространение частиц примеси
Клеточный автомат, использующийся в данном эксперименте, имеет те же размеры, что и в предыдущих экспериментах. Начиная с того момента когда автомат полностью заполнен чистыми частицами, в автомат начинают поступать частицы примеси, которые подразумевают собой некоторое загрязняющее вещество. Поступают частицы примеси из клеток источников с координатами [(2, 1), (22, 1)], в течение 50 итераций. Затем поступление частиц примеси прекращается. Данный эксперимент показывает, как ведут себя частицы с новым свойством. Голубым цветом отображаются концентрации чистых частиц. Красным цветом изображаются области автомата с максимальной заданной концентрацией примеси. На рисунках 19 — 28 видно, что частицы примеси стремятся опуститься вниз автомата, будто на них действует сила тяжести. Другие эксперименты представлены в приложении А.
Рисунок 19 — Концентрация примеси через 10 итераций
Рисунок 20 — Концентрация примеси через 55 итераций
Рисунок 21 — Концентрация примеси через 90 итераций
Рисунок 22 — Концентрация примеси через 125 итераций
Рисунок 23 — Концентрация примеси через 160 итераций
Рисунок 24 — Концентрация примеси через 195 итераций
Рисунок 25 — Концентрация примеси через 230 итераций
Рисунок 26 — Концентрация примеси через 255 итераций
Рисунок 27 — Концентрация примеси через 280 итераций
Рисунок 28 — Концентрация примеси через 315 итераций
ЗАКЛЮЧЕНИЕ
В данной дипломной работе проведено исследование клеточно-автоматной модели потока жидкости FHP-I. Экспериментально показано соответствие модели FHP-I потоку реальной жидкости через сравнения результатов моделирования с аналитическими расчетами распределения скорости потока в трубе по закону Пуазейля.
В результате экспериментов с препятствиями получены потоки с завихрениями, что является признаком турбулентного потока. Следовательно, качественно эта модель соответствует физическим явлениям и позволяет моделировать турбулентные потоки жидкости.
В модель был введен новый вид частиц для возможности моделирования потоков с примесью, что потребовало модификации исходной модели. Были получены изображения автомата, в которых наблюдалось распространение частиц примеси под действием силы тяжести.
СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
Фон Нейман Дж. Теория самовоспроизводящихся автоматов / Фон Нейман Дж. — М.: Мир, 1971. — 500 с.
Frish U. Lattice-gas Automata for Navier-Stokes equation / Frish U. San-Francisco: Mathematical Physics, 1986. 1530 p.
Беркович С.Я. Клеточные автоматы как модель реальности: поиски представлений физических и информационных процессов. — М.: МГУ, 1993. — 112 с.
Frish U., Crutch J.P. Lattice gas hydrodynamics in two and three dimensions / Frish U., Crutch J.P. New -York : Mathematical Physics, 1990. 117 p.
Медведев Ю.Г. Параллельная реализация трехмерной клеточно-автоматной модели потока жидкости. — Новосибирск : ИВМиМГ, 2004. — 112с.
Бандман О.Л. Клеточно-автоматное моделирование пространственной динамики / О.Л. Бандман. — Новосибирск.: СО РАН, 2000. — 113 c.
Медведев Ю.Гр Исследование параллельной реализации клеточно-автоматной модели кумулятивной струи / Ю.Г. Медведев. — Томск.: ТНЦСОРАН, 2000. — 120 c.
ПРИЛОЖЕНИЕ А
Рисунок А1 — Обтекание произвольных препятствий
ПРИЛОЖЕНИЕ Б
Программный код программы СА_FHP-Й Simulator реализующий модель FHP-I.
struct chst — структура частицы
{
bool c; — признак что это частица.
bool g; — признак что это примесь.
};
struct kl
{
chst in[7]; массив для хранения состояния клетки, также 0-й элемент указывает на то, что клетка является препятствием.
chst out[7]; -массив для вычисления новых состояний клетки.
};
chst ch; — частица для генерации новых частицы.
Vector <obst> o; — вектор для хранения координат клеток являющимися препятствием.
struct obst
{
int i;
int j;
};
Так как клетки находятся на позиции с индексами такими, чтобы выполнялось условие: (i+j) mod 2=0, то фрагменты кода программы будут представлены для клеток с нечетными индексами по строкам и столбцам.
Генерация частиц:
for ( i = 1; i <= 97; i += 2 )
{
ch.c=rand() % 2;
if(ch.c==1)
ch.g=rand() % 2;
m[i][3].out[3] = ch;
ch.c=rand() % 2;
if(ch.c==1)
ch.g=rand() % 2;
m[i][3].out[4] = ch;
ch.c=rand() % 2;
if(ch.c==1)
ch.g=rand() % 2;
m[i][3].out[5] = ch;
}
Фаза сдвига частиц:
for ( i = 1; i <= 97; i += 2 )
{
for ( j = 3; j <= 3 + Iter; j += 2 )
{
m[i][j — 2].in[1] = m[i][j].out[1];
m[i + 1][j — 1].in[2] = m[i][j].out[2];
m[i + 1][j + 1].in[3] = m[i][j].out[3];
m[i][j + 2].in[4] = m[i][j].out[4];
m[i — 1][j + 1].in[5] = m[i][j].out[5];
m[i — 1][j — 1].in[6] = m[i][j].out[6];
}
}
Фаза столкновения частиц:
for ( i = 2; i <= 98; i += 2 )
{
for ( j = 4; j <= 4 + Iter; j += 2 )
{
if (m[i][j].in[0] != 1)
{
int ss = 0;
for ( int ks = 1; ks <= 6; ks++ )
ss += m[i][j].in[ks];
ss += m[i][j].out[0];
if (ss == 2)
goto ch2;
if (ss == 3)
goto ch3;
if (ss == 4)
goto ch4;
if (ss == 5)
goto ch5;
goto dru;
ch3: //3 частицы
c2(i,j,0,0,0,1,0,1,1, 1,0,0,0,1,0,1, 0,1,0,0,1,1,0)
c2(i,j,0,0,1,0,0,1,1, 1,1,0,0,0,1,0, 0,1,0,0,1,0,1)
c2(i,j,0,0,0,1,1,0,1, 1,0,0,1,0,1,0, 0,0,1,0,1,1,0)
c4(i,j,0,0,1,0,1,0,1, 1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 1,1,0,0,1,0,0, 0,1,0,1,0,1,0)
c4_2(i,j,0,0,1,0,1,0,1, 1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 1,1,0,0,1,0,0, 0,1,0,1,0,1,0)
c2(i,j,0,1,0,0,1,0,1, 1,1,0,0,0,1,0, 0,0,1,0,0,1,1)
c2(i,j,1,0,0,0,1,0,1, 0,0,0,1,0,1,1, 0,1,0,0,1,1,0)
c2(i,j,0,0,1,1,0,0,1, 0,1,1,0,1,0,0, 1,1,0,1,0,0,0)
c2(i,j,0,1,0,1,0,0,1, 1,0,1,0,0,0,1, 0,1,1,0,0,1,0)
c4(i,j,1,0,0,1,0,0,1, 1,0,1,0,0,1,0, 1,1,0,0,1,0,0, 0,1,0,1,0,1,0, 0,0,1,0,1,0,1)
c4_2(i,j,1,0,0,1,0,0,1, 1,0,1,0,0,1,0, 1,1,0,0,1,0,0, 0,1,0,1,0,1,0, 0,0,1,0,1,0,1)
c2(i,j,1,0,1,0,0,0,1, 0,1,1,0,0,1,0, 0,1,0,1,0,0,1)
c2(i,j,0,0,1,0,1,1,0, 0,0,0,1,1,0,1, 1,0,0,1,0,1,0)
c2(i,j,0,1,0,0,1,1,0, 0,0,0,1,0,1,1, 1,0,0,0,1,0,1)
c2(i,j,0,0,1,1,0,1,0, 1,0,1,0,1,0,0, 0,1,0,1,1,0,0)
c4(i,j,0,1,0,1,0,1,0, 1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 1,1,0,0,1,0,0, 0,0,1,0,1,0,1)
c4_2(i,j,0,1,0,1,0,1,0, 1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 1,1,0,0,1,0,0, 0,0,1,0,1,0,1)
c2(i,j,1,0,0,1,0,1,0, 0,0,1,0,1,1,0, 0,0,0,1,1,0,1)
c2(i,j,0,1,1,0,0,1,0, 1,0,1,0,0,0,1, 0,1,0,1,0,0,1)
c4(i,j,1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 1,1,0,0,1,0,0, 0,1,0,1,0,1,0, 0,0,1,0,1,0,1)
c4_2(i,j,1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 1,1,0,0,1,0,0, 0,1,0,1,0,1,0, 0,0,1,0,1,0,1)
c2(i,j,1,1,0,0,0,1,0, 0,1,0,0,1,0,1, 0,0,1,0,0,1,1)
c2(i,j,0,1,0,1,1,0,0, 0,0,1,1,0,1,0, 1,0,1,0,1,0,0)
c2(i,j,0,1,1,0,1,0,0, 1,1,0,1,0,0,0, 0,0,1,1,0,0,1)
c2(i,j,1,0,1,0,1,0,0, 0,1,0,1,1,0,0, 0,0,1,1,0,1,0)
c4(i,j,1,1,0,0,1,0,0, 1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 0,1,0,1,0,1,0, 0,0,1,0,1,0,1)
c4_2(i,j,1,1,0,0,1,0,0, 1,0,1,0,0,1,0, 1,0,0,1,0,0,1, 0,1,0,1,0,1,0, 0,0,1,0,1,0,1)
c2(i,j,1,1,0,1,0,0,0, 0,0,1,1,0,0,1, 0,1,1,0,1,0,0)
goto dru;
ch2: //2 частицы
//аналогично
ch4: //4 частицы
//—//
ch5: //5 частиц
//—//
dru:
for ( int ks = 1; ks <= 6; ks++ )
m[i][j].out[ks] = m[i][j].in[ks];
goto endend;
end:
for(ri=0;ri<=6;ri++)
{
if(m[(i)][(j)].out[ri].c==1)
{
if(out0.c==1)
{
m[(i)][(j)].out[ri].g=out0.g;
out0.c=0;
}
else
for(rj=1;rj<=6;rj++)
{
if(m[(i)][(j)].in[rj].c==1)
{
m[(i)][(j)].out[ri].g=m[(i)][(j)].in[rj].g;
m[(i)][(j)].in[rj].c=0;
break;
}
}
}
}
endend:
}
}
}
ci-макрос в теле которого происходит изменение состояния клетки. Макрос c4_2 продолжение макроса c4.
#define c4(i,j,i0,i1,i2,i3,i4,i5,i6,o10,o11,o12,o13,o14,o15,o16,o20,o21,o22,o23,o24,o25,o26,o30,o31,o32,o33,o34,o35,o36,o40,o41,o42,o43,o44,o45,o46)
if((m[(i)][(j)].out[0]==(i0))
&(m[(i)][(j)].in[1]==(i1))
&(m[(i)][(j)].in[2]==(i2))
&(m[(i)][(j)].in[3]==(i3))
&(m[(i)][(j)].in[4]==(i4))
&(m[(i)][(j)].in[5]==(i5))
&(m[(i)][(j)].in[6]==(i6)))
{ t=rand()%4;
if(t==0)
{
m[(i)][(j)].out[0]=(o10);
m[(i)][(j)].out[1]=(o11);
m[(i)][(j)].out[2]=(o12);
m[(i)][(j)].out[3]=(o13);
m[(i)][(j)].out[4]=(o14);
m[(i)][(j)].out[5]=(o15);
m[(i)][(j)].out[6]=(o16);
goto end;
}
if(t==1)
{
m[(i)][(j)].out[0]=(o20);
m[(i)][(j)].out[1]=(o21);
m[(i)][(j)].out[2]=(o22);
m[(i)][(j)].out[3]=(o23);
m[(i)][(j)].out[4]=(o24);
m[(i)][(j)].out[5]=(o25);
m[(i)][(j)].out[6]=(o26);
goto end;
}
if(t==2)
{
m[(i)][(j)].out[0]=(o30);
m[(i)][(j)].out[1]=(o31);
m[(i)][(j)].out[2]=(o32);
m[(i)][(j)].out[3]=(o33);
m[(i)][(j)].out[4]=(o34);
m[(i)][(j)].out[5]=(o35);
m[(i)][(j)].out[6]=(o36);
goto end;
}
if(t==3)
{
m[(i)][(j)].out[0]=(o40);
m[(i)][(j)].out[1]=(o41);
m[(i)][(j)].out[2]=(o42);
m[(i)][(j)].out[3]=(o43);
m[(i)][(j)].out[4]=(o44);
m[(i)][(j)].out[5]=(o45);
m[(i)][(j)].out[6]=(o46);
goto end;
}
}
Фаза столкновения частиц с препятствием:
o_sz = size();
for (int io = 0; io <= o_sz — 1; io++ )
{
m[o[io].i][o[io].j — 2].in[1] = m[o[io].i][o[io].j].in[4];
m[o[io].i + 1][o[io].j — 1].in[2] = m[o[io].i][o[io].j].in[5];
m[o[io].i + 1][o[io].j + 1].in[3] = m[o[io].i][o[io].j].in[6];
m[o[io].i][o[io].j + 2].in[4] = m[o[io].i][o[io].j].in[1];
m[o[io].i — 1][o[io].j + 1].in[5] = m[o[io].i][o[io].j].in[2];
m[o[io].i — 1][o[io].j — 1].in[6] = m[o[io].i][o[io].j].in[3];
}
Фаза воздействия силы тяжести:
for ( i = 1; i <= 97; i += 2 )
{
for ( j = 3; j <= 3 + k; j += 4 )
{
for ( p2=1;p2<=6;p2++)
{
if ((m[i+1][j+1].out[p2].c==1) && (m[i+1][j+1].out[p2].g==1) &&
(m[i+1][j+1].in[0].c==0) )
{
if((m[i][j+2].out[p2].c==1) && ( m[i][j+2].out[p2].g==0) &&
(m[i][j+2].in[0].c==0) )
{
m[i][j+2].out[p2].c=1;
m[i][j+2].out[p2].g=1;
m[i+1][j+1].out[p2].c=1;
m[i+1][j+1].out[p2].g=0;
}
else if ((m[i][j].out[p2].c==1) && (m[i][j].out[p2].g==0) &&
(m[i][j].in[0].c==0) )
{
m[i][j].out[p2].c=1;
m[i][j].out[p2].g=1;
m[i+1][j+1].out[p2].c=1;
m[i+1][j+1].out[p2].g=0;
}
else if ((m[i][j+2].out[p2].c==0) && ( m[i][j+2].out[p2].g==0) && (m[i][j+2].in[0].c==0) )
{
m[i][j+2].out[p2].c=1;
m[i][j+2].out[p2].g=1;
m[i+1][j+1].out[p2].c=0;
m[i+1][j+1].out[p2].g=0;
}
else if ((m[i][j].out[p2].c==0) && (m[i][j].out[p2].g==0) &&
(m[i][j].in[0].c==0) )
{
m[i][j].out[p2].c=1;
m[i][j].out[p2].g=1;
m[i+1][j+1].out[p2].c=0;
m[i+1][j+1].out[p2].g=0;
}
