Содержание к диссертации
Введение
Глава 1. Общие принципы клеточно-автоматного моделирования потоков жидкости 15
1.1. Основные определения 15
1.2. Клеточно-автоматные модели потоков жидкости 17
1.3. Двумерная модель с 4 соседями (НРР модель) 20
1.4. Двумерная модель с 6 соседями (FHP модель) 22
1.5. Четырехмерная модель с 24 соседями (FCHC модель) 28
1.6. Моделирование потоков жидкости с меньшей вязкостью. Модификация FHP модели 29
Выводы к главе 34
Глава 2. Трехмерная КА модель потока жидкости RD-I с 12 соседями 35
2.1. Структура модели 35
2.2. Правила перехода клеточного автомата 40
2.3. Граничные условия модели 42
2.4. Осреднение значений 44
2.5. Модельные параметры. Соотношения модельных и физических величин 46
В ыводы к главе 5 О
Глава 3. Программный комплекс моделирования потоков жидкости 51
3.1. Назначение и блок-схема комплекса 51
3.2. Основные модули комплекса 53
3.3. Симулятор потока жидкости 57
3.4. Параллельная реализация симулятора 60
3.5. Форматы данных 63
3.6. Методика проведения вычислительных экспериментов 65
Выводы к главе 71
Глава 4. Эксперименты с моделью RD-I 72
4.1. Алгоритм моделирования стационарного потока 73
4.2. Поток вязкой жидкости в трубе. Сравнение с законом Пуазейля 76
4.3. Поток вязкой жидкости в трубе с задвижкой 84
4.4. Моделирование пористых сред. Сравнение с законом Дарси 87
4.5. Оценки эффективности распараллеливания 91
Выводы к главе 95
Заключение 96
Список литературы
- Двумерная модель с 4 соседями (НРР модель)
- Правила перехода клеточного автомата
- Параллельная реализация симулятора
- Поток вязкой жидкости в трубе. Сравнение с законом Пуазейля
Введение к работе
Общая схема решения задач моделирования потоков жидкости обычно состоит их трех этапов. На первом этапе производится поиск математической модели, отражающей изучаемый процесс. Второй этап состоит в выборе численных методов решения построенной модели. На третьем этапе проводят вычислительные эксперименты на компьютере или многопроцессорной вычислительной системе. Существует два принципиально различающихся подхода к решению подобных задач.
В традиционном подходе в качестве модели изучаемого процесса используются дифференциальные уравнения в частных производных. Дифференциальные уравнения, используемые в реальных задачах, невозможно решить аналитически, поэтому возникает необходимость применять численные методы их решения, с которыми при больших размерах задачи связан ряд проблем. [1, 2, 3]. Так, например, в большинстве задач при использовании явных схем решение сходится очень медленно из-за малого шага по времени, поэтому приходится использовать неявные схемы. При вычислении неявных схем на многопроцессорных вычислительных системах трудно достигнуть высокой эффективности распараллеливания.
Поиск путей преодоления вышеперечисленных трудностей ведется по двум направлениям: совершенствование численных методов решения дифференциальных уравнений и способов их распараллеливания и поиск новых моделей, отличных от дифференциальных уравнений [4, 5,6]. Одно из направлений этого поиска осуществляется путем создания и исследования новых моделей, обладающих следующими свойствами.
Дискретность. Модели должны иметь булев базис для возможности прямого отображения на структуру вычислителя, легкой аппаратной реализации на базе современной микроэлектроники, а также для возможности эмуляции на универсальной ЭВМ.
Природные свойства. В основу модели должен быть положен принцип близко действия и другие физические законы, чтобы свойства модельных объектов были максимально схожими с физическими свойствами природных объектов, а интерпретация результатов моделирования была простой.
Пространственный параллелизм. Модели должны обладать естественным параллелизмом, чтобы иметь возможность легкой реализации на спецпроцессорах, а также легко распараллеливаться на крупноблочных многопроцессорных вычислительных системах.
Один из новых классов моделей, разрабатываемых для описания широкого круга задач пространственной динамики, в том числе и для моделирования потоков жидкостей, и обладающих вышеперечисленными свойствами составляют модели, основанные на клеточных автоматах [7, 8]. Впервые клеточный автомат был предложен Дж. Фон Нейманом (J. von Neumann) в 1948 году [9]. Затем эта модель развивалась и находила применение в различных областях, в том числе и в задачах моделирования пространственной динамики. Приведем некоторые важные даты, связанные с развитием клеточно-автоматных моделей потоков жидкости.
1973 год. J. Hardy, О. de Pazzis, Y. Pomeau. Предложена двумерная НРР модель потока [10]. Клетки этой модели имеют четыре соседа. В модели выполняются «лишние» законы сохранения, что приводит к ее анизотропии. Модель не получила применения для описания потоков жидкости. год. U. Frisch, В. Hasslacher, Y. Pomeau. Предложена двумерная FHP модель с шестью соседями [11]. Эта модель изотропна, ее элементарные автоматы имеют вероятностные правила переходов. Доказано, что FHP модель аппроксимирует уравнение Навье-Стокса и может быть использована для моделирования потоков жидкости на плоскости. год. Т. Toffoli, N. Margolus. Построен спецпроцессор, аппаратно реализующий клеточный автомат, так называемая машина клеточных автоматов САМ-8 (Cellular Automata Machine) [12]. Создание этого устройства дало толчок к поиску новых моделей и развитию клеточно-автоматного моделирования как самостоятельной отрасли науки.
1987 год. U. Frisch, D. d'Humieres, В. Hasslacher, P. Lallemand, Y. Pomeau, J.-P. Rivet. Предложена четырехмерная FCHC (Face Centered Hyper Cubic) модель, построенная на гранецентрированном гиперкубе [13]. Проекция этой вероятностной модели с 24 соседями на трехмерное пространство адекватно описывает трехмерные потоки жидкости и полностью удовлетворяет условиям изотропии. Несмотря на это, FCHC модель не получила широкого распространения в силу большой вычислительной сложности.
Кроме задач газовой [14 - 21] и гидродинамики [22 - 29] клеточные автоматы применяются во многих других областях для моделирования самых различных физических процессов, таких как пространственная диффузия [30 ~ 33], эволюция спиральных галактик [34], нелинейные химические системы [35, 36, 37], рост ветвящихся кристаллов [38, 39], фазовые переходы [40, 41], пространственная структура турбулентных потоков [42], рост, поведение и функционирование биологических организмов [43, 44], обработка изображений, распознавание образов [45], движение групп людей и толпы [47 - 52], распространение эпидемий [53], сейсмические волны [54], магнитогидродинамические явления [55, 56] и многое другое [57, 58, 59]. В связи с обилием вышеперечисленных применений клеточно-автоматные модели приобретают все большую актуальность.
Важным условием эффективности моделирования является способ реализации модели [60]. Клеточные автоматы могут быть реализованы на универсальных ЭВМ [61 — 64], на многопроцессорных вычислительных системах [65], в виде специализированных машин клеточных автоматов [12, 66] и даже известны проекты реализации клеточных автоматов на молекулярном уровне [3].
Реализация клеточных автоматов на универсальных ЭВМ является самым доступным способом реализации. Программные модели имеют -7-a.^ низкую стоимость, просты в эксплуатации и могут быть легко модифицированы. Но у этого способа реализации есть существенный недостаток — большое время вычислений.
Для сокращения времени вычислений используют параллельную реализацию. Клеточно-автоматные модели обладают высокой эффективностью распараллеливания на многопроцессорных вычислительных системах и системах с массовым параллелизмом [12, 65, 66]. Параллелизм заложен в самой структуре этих моделей, поэтому они реализуются на любых многопроцессорных архитектурах с сотнями, тысячами и даже десятками тысяч процессоров при минимальной деградации эффективности распараллеливания. " Самым эффективным способом реализации является построение спецпроцессоров. Каждая клетка автомата реализуется конкретным элементом схемы. Этот элемент имеет невысокую сложность, поэтому на современном уровне развития микроэлектроники интегральная схема может реализовать клеточный автомат достаточно большого размера [23].
Задание граничных условий в клеточно-автоматных моделях осуществляется заданием типов клеток автомата. Группы клеток, образующие границы, состоят из клеток типа стенки. Хотя поведение клеток стенки и отличается от поведения обычных клеток, вычислительная ф: сложность обработки одной клетки не зависит от ее типа. Это означает, что время вычислений клеточного автомата заданного размера остается неизменным, какими бы сложными ни были его границы. В силу легкости задания граничных условий клеточно-автоматные модели приобрели особенное значение для исследования потоков вязкой жидкости через пористые среды [67 - 74]. Также клеточно-автоматные модели позволяют моделировать движение нефти в трубопроводах сложной конфигурации и в грунте. ф Моделирование трехмерных потоков является чрезвычайно сложной задачей. С одной стороны, модель должна удовлетворять условиям -8-изотропии, а с другой стороны — иметь приемлемую вычислительную сложность, которая зависит от количества соседних клеток. Условия изотропии накладываются на тензоры изотропии до четвертого порядка включительно [75, 76]. Среди трехмерных моделей нет таких, которые обладали бы изотропией четвертого порядка. Были найдены четырехмерные модели, но они оказались слишком сложными для практического использования [13]. В связи с этим возникла задача поиска трехмерной модели, в которой условия изотропии выполнялись бы приближенно (до третьего порядка), и которая была бы приемлемой по сложности. Предлагаемая в диссертации трехмерная клеточно-автоматная модель, названная RD-I, имеет третий порядок точности соответствия с условиями изотропии, но ее вычислительная сложность на несколько порядков ниже, чем у точных моделей.
В связи с приведенными выше тезисами трехмерное клеточно-автоматное моделирование потоков приобретает большую актуальность.
Цель работы.
Целью диссертационной работы является разработка трехмерной клеточно-автоматной модели потока вязкой жидкости с приемлемой сложностью, которая обладала бы следующими свойствами: малая по % сравнению с известными моделями сложность реализации, эффективная распараллеливаемость программной реализации на многопроцессорных вычислительных комплексах, возможность моделирования пористых сред, а также разработка программного комплекса этой модели и экспериментальное исследование ее вычислительных характеристик.
Задачи.
В соответствии с поставленной целью в диссертационной работе решаются следующие задачи: - разработка алгоритма моделирования потока вязкой жидкости на основе клеточно-автоматной модели с трехмерной структурой и вероятностными правилами переходов автомата; - создание программного комплекса, эмулирующего клеточный автомат разработанной модели и позволяющего производить вычислительные эксперименты с моделью как на последовательных компьютерах, так и на многопроцессорных вычислительных комплексах; - экспериментальное исследование разработанной модели, определение ее параметров и соотношений, связывающих параметры модели с параметрами моделируемого потока; исследование производительности программного комплекса, а также эффективности его распараллеливания на вычислительных кластерах; разработка алгоритмов задания краевых условий экспериментов и экспериментальное исследование потоков жидкости через трехмерную пористую среду.
Методы исследования.
При проведении исследований использовался аппарат дискретной математики, теории автоматов, параллельного программирования и теории вероятностей, а также вычислительный эксперимент на одно- и многопроцессорных ЭВМ.
Научная новизна.
В ходе исследования были получены следующие новые результаты: предложена новая трехмерная клеточно-автоматная модель, названная RD-I, отличающаяся от известных моделей много меньшей сложностью; получены соотношения между физическими и модельными параметрами, сопоставляющие результаты моделирования с реальными величинами, характеризующими моделируемый поток; разработан метод задания краевых условий для предложенной модели, позволяющий задавать краевые условия объектов моделирования с границами любой сложности, включая пористые среды; создан программный комплекс, эмулирующий эволюцию клеточного автомата предложенной модели как на последовательных так и на многопроцессорных вычислителях.
Практическая значимость работы.
Предложенная трехмерная клеточно-автоматная модель RD-I может быть использована для моделирования потоков вязкой жидкости в объектах с границами любой сложности. Созданный в рамках работы программный комплекс может применяться для моделирования потоков вязкой жидкости, в том числе и в пористых средах.
Достоверность результатов.
Достоверность полученных результатов обусловлена использованием аппарата дискретной математики и теории автоматов для построения модели, для получения параметров модели и соотношений, связывающих модельные и физические параметры. Кроме того, все результаты исследований подтверждены вычислительными экспериментами с использованием разработанного программного комплекса.
Основные положення, выносимые на защиту.
Предложенная трехмерная клеточно-автоматная модель потока вязкой жидкости RD-I со структурой автомата в виде ромбического додекаэдра и вероятностными правилами перехода является моделью потока вязкой жидкости.
Сложность граничных условий модели RD-I (включая случаи с пористыми средами) не оказывает существенного влияния на вычислительную сложность алгоритма. *' 3. Результаты моделирования с помощью RD-I укладываются в пределы допустимой погрешности относительно аналитических расчетов. Модельное и физическое числа Рейнольдса также различаются на допустимую величину погрешности. Предел допустимой погрешности находится в обратной зависимости от * используемого радиуса осреднения.
4. Созданный программный комплекс эмулирует эволюцию клеточного автомата модели RD-I как на последовательных ЭВМ, так и на многопроцессорных вычислительных комплексах.
Содержание диссертации изложено в четырех главах. * Первая глава посвящена изложению общих принципов построения клеточно-автоматных моделей. В ней вводятся основные понятия, такие как клеточный автомат, элементарный автомат (клетка), частица, вектор состояния элементарного автомата, отношение соседства. В этой главе даны оценки сложности модели в зависимости от количества соседних клеток. Описаны двумерные НРР модель с 4 соседями, FHP модель с 6 соседями и четырехмерная FCHC модель с 24 соседями. Описан способ получения модельной скорости путем осреднения векторов скорости частиц. Также в этой главе описана программная реализация и приведены результаты щ, экспериментальных исследований FHP модели и предложена ее модификация, позволяющая увеличить число Рейнольдса моделируемого потока. Результаты этих экспериментов послужили основой для построения трехмерной модели и для формирования экспериментальных данных для ее исследования.
Во второй главе диссертации предлагается трехмерная КА модель потока жидкости RD-I с 12 соседями и одной частицей покоя в каждой клетке [77]. Структура автомата основана на полярном комплексе ромбического щ додекаэдра. Если заполнить трехмерное пространство ромбическими додекаэдрами с единичным расстоянием между двумя противолежащими -12-гранями, то координаты узлов решетки автомата будут совпадать с координатами центров додекаэдров. При этом соседние узлы будут находиться в центрах додекаэдров, имеющих общую грань. Узлы решетки при данном пространственном распределении совпадают с центрами плотно упакованных шаров в одном из вариантов плотной упаковки. Предложенная модель RD-I достаточно проста в реализации, т.к. производит операции над булевыми векторами. Также в этой главе изложен метод задания граничных условий [78] и выведены формулы для вычисления параметров, таких как модельная вязкость и структурный коэффициент [79].
В третьей главе диссертации представлен программный комплекс, созданный для реализации трехмерной клеточно-автоматной модели RD-I на последовательных компьютерах и параллельных вычислительных системах [80]. В ней сказано о целях создания и о назначении программного комплекса, описаны форматы данных, используемых при моделировании [81], алгоритмы ввода-вывода, алгоритм эмуляции элементарного автомата, алгоритм осреднения, алгоритм распределения клеток автомата по процессорам в параллельной реализации. Также в этой главе приведены рекомендации по использованию комплекса.
Четвертая глава диссертации содержит описание вычислительных экспериментов [82]. Серия экспериментов по моделированию потока вязкой жидкости в трубе круглого сечения была поставлена для того, чтобы подтвердить соответствие модели RD-I реальному потоку, скорость которого в трубе можно рассчитать аналитически из закона Пуазейля [83]. Далее в четвертой главе описан эксперимент по моделированию движения нефти с такими же физическими характеристиками в трубе такого же размера, как и в предыдущем эксперименте, но с наполовину закрытой задвижкой, имеющей форму полукруга в поперечном сечении трубы. Еще одна серия экспериментов, описанная в четвертой главе диссертации, представляет собой моделирование потока вязкой жидкости через пористую среду. В них исследовалась зависимость скорости потока от плотности пор. Результаты -13-этих экспериментов коррелируют с законом Дарси. Также в этой главе даны оценки эффективности распараллеливания, полученные в результате серии экспериментов с использованием параллельной версии программного комплекса на мультикомпьютерах.
Реализация результатов работы.
Работа выполнена в соответствии с темой N 25 «Исследование физических и информационных процессов с использованием моделей мелкозернистого параллелизма» (№ госрегистрации 01.20.00.04430) плана научно-исследовательских работ отдела математического обеспечения высокопроизводительных вычислительных систем института вычислительной математики и математической геофизики СО РАН, а также поддержана грантом РФФИ № 01.00.00026 (1999г.) и Программой фундаментальных исследований РАН № 17 (2004г.)
Результаты диссертационной работы были использованы в учебном процессе на кафедре параллельных вычислительных технологий Новосибирского государственного технического университета, а также для решения ряда типовых задач моделирования потоков.
Личный вклад автора.
Все результаты, представленные в диссертационной работе, получены автором лично.
Апробация работы.
Основные результаты диссертации докладывались и обсуждались на конгрессе по прикладной и индустриальной математике (Новосибирск, 1998), на конференции, посвященной 70-летию Сибирского физико-технического НИИ (Томск, 1998), на конференции с международным участием «Новые информационные технологии в исследовании дискретных структур» (Томск, 2000), на международной конференции по вычислительной математике -14-(Новосибирск, 2001), на конференции «Моделирование неравновесных систем» (Красноярск, 2001), на Сибирской научной школе-семинаре по параллельному программированию (Томск, 2001), на конференциях молодых ученых ИВМиМГ (Новосибирск, 1999; 2002; 2004), на конференциях с международным участием «Новые информационные технологии в исследовании сложных структур» (Томск, 2002; Иркутск, 2004), а также на семинарах отделов Математического обеспечения высокопроизводительных вычислительных систем и Математических задач геофизики Института вычислительной математики и математической геофизики СО РАН.
Публикации.
Результаты исследований по теме диссертационной работы опубликованы в виде 3 статей в научных журналах, 8 докладов, 3 тезисов выступлений на научных конференциях. Всего опубликовано 14 работ.
Структура и объем работы.
Диссертация состоит из введения, четырех глав, заключения и списка цитируемой литературы. Содержание диссертации изложено на 108 страницах, содержит 41 иллюстрацию и 8 таблиц. Библиографический список включает 127 наименований.
Двумерная модель с 4 соседями (НРР модель)
В этой модели структура клеточного автомата имеет вид двумерной квадратной решетки. Количество клеток-соседей, таким образом, равно четырем. Каждый такт работы автомата, как следует из (1.1), состоит из двух фаз: сдвига и столкновения. При сдвиге, как и во всех остальных моделях, происходит перемещение всех частиц в соседние клетки в направлении их векторов скорости (переход фрагмента КА из состояния, изображенного на рис. 1.1, в состояние, изображенное на рис. 1.2). Правила столкновения в этом автомате детерминированные и определяются следующим образом. Если в какой-то момент времени в клетке оказываются ровно две частицы, движущиеся во встречных направлениях, то они меняют направления своего движения на два противоположных ранее свободных направления (переход фрагмента КА из состояния, изображенного на рис. 1.2, в состояние, изображенное на рис. 1.3). Формально это правило можно определить следующим образом. Если состояние клетки s = 0101, то оно переходит в состояние s\ = 1010, если состояние s = 1010, то оно переходит BSI =0101, иначе состояние клетки не меняется. Приведенные правила являются единственными нетривиальными правилами столкновения (такими правилами, в которых выходное состояние не совпадает с входным), при которых сохраняются масса и импульс частиц в клетке. Нетривиальные правила столкновений для НРР модели изображены на рис. 1.4. Более того, в НРР модели сохраняются импульсы частиц вдоль горизонтальной и вдоль вертикальной осей. Это свойство является «лишним» законом сохранения. Из-за него НРР модель дает большую погрешность при моделировании потоков жидкости [75, 98], но это не мешает применять ее для моделирования других процессов, как, например, распространение электромагнитных и акустических волн [99, 100]. Поведение НРР модели инвариантно относительно дискретных преобразований, таких как поворот -на тс/4, зеркальная симметрия и сдвиг по горизонтали и по вертикали на шаг решетки. Более того, поведение инвариантно относительно дополнения, когда 1 и 0 (частицы и дыры) меняются местами.
FHP модель так же, как и НРР, является двумерной моделью потока. Ее клетки образуют гексагональную структуру, т. е. каждая клетка имеет ровно шесть соседей. Известен целый класс FHP моделей с одной, двумя и тремя частицами покоя, и длиной вектора состояния клетки семь, восемь и девять соответственно [25, 101]. Их правила столкновения построены по схожим принципам, поэтому далее рассматривается только модель с одной частицей покоя [11]. Каждая клетка является конечным недетерминированным автоматом, и ее состояние можно представить в виде семиразрядного булева вектора (рис. 1.5). Нетривиальные (изменяющие направления движения частиц) правила столкновений для FHP модели приведены в таблице 1.1. В случаях, когда имеется два варианта состояния после столкновения, выбирается одно из них равновероятно. Всего правил столкновения в FHP модели 27, т.к. вектор состояния клеток в ней семиразрядный.
Графически нетривиальные правила столкновений изображены на рис. 1.6. Частицы покоя выделены на нем черным цветом. Правила, получающиеся из изображенных состояний преобразованиями симметрии, на рисунке не приведены, хотя в таблице 1.1 они присутствуют. На рис. 1.6а показаны вероятностные правила столкновений, а на рис. 1.66 — детерминированные. Эти правила сохраняют суммарные массу и импульс -частиц в клетке (условия (1.4) и (1.5)). В отличие от НРР, в FHP модели нет «лишних» законов сохранения, из рис. 1.56 видно, что импульс вдоль осей не сохраняется. На рис. 1.7, 1,8 и 1.9 показан фрагмент автомата FHP модели и его состояния после сдвига и после столкновения соответственно.
Поведение FHP модели инвариантно относительно дискретных преобразований, таких как поворот на я/6, зеркальная симметрия и сдвиг на шаг решетки по каждому из трех осевых направлений. Поведение также инвариантно относительно дополнения, когда 1 и 0 (частицы и дыры) меняются местами.
Практический интерес для исследователя представляют осредненные значения скоростей и концентраций частиц в некоторой окрестности Av(w). Они вычисляются по формулам (1.8) и (1.9). Окрестность Av(w) для FHP модели имеет форму шестиугольника (рис 1.10). На рисунке жирными линиями выделены клетки, попадающие в окрестность осреднения. Вектор осредненной скорости изображен в центре окрестности.
FHP модель описывает потоки жидкостей и газов также как и уравнение Навье-Стокса. В работе [22], где приведено доказательство соответствия FHP модели уравнению Навье-Стокса, в котором формализм Грина-Кубо, позволяющий представлять коэффициенты переноса макроскопической системы в терминах равновесных корреляционных функций, расширен на КА модели гидродинамики. Другое доказательство этого соответствия приведено в [102], где уравнение Навье-Стокса выведено из уравнений, описывающих поведение клеточного автомата на микроуровне, методом, использующим разложение Чепмена-Энскога.
Одна из характеристик КА моделей — набор тензоров изотропии различного ранга. Они показывают степень инвариантности модели относительно различных вращений и отражений. В КА моделях потоков рассматривают тензоры r-го ранга: где с/сГ/- проекции вектора С{, на Декартову ось а, є {х,у, z}, (і = 1, 2, ..., г), a /= 1, 2, ..., & - индекс суммирования по всем возможным направлениям скоростей вектора С/.
Для того чтобы модель была изотропной, на тензоры от первого до четвертого ранга включительно накладывают следующие условия [3, 23, 102], Тензоры первого ранга должны удовлетворять равенству: Ec/«=0,Vae {x,y,z}. (1.11) В тензорах второго ранга суммируются все попарные произведения ПроеКЦИЙ Ciai he2 Zc c = 5- »Va,y9e {x,y,z}. (1.12)
Малыми греческими буквами обозначают произвольно фиксированные Декартовы координаты а, Д у, Г) є {х, у, z}, D - размерность Декартова пространства, Ь - длина вектора состояния клетки, с - модуль вектора „ „ (0, еслиа р скорости, а о — символ Кронекера, т.е. SaS =\ [\, еслиаг = /ї Тензоры третьего ранга должны удовлетворять условию: ZC W=. (1.13) ХсЛс гс/7 = д(д + 2) + Say6p5 + 6asS - (1 Л4) Для тензоров четвертого ранга должно выполняться: __4Ьс
Считается, что модель изотропна, если ее тензоры изотропии удовлетворяют условиям (1.11) — (1.14). Для FHP модели эти равенства выполняются [11].
Правила перехода клеточного автомата
Одно из преимуществ клеточно-автоматных моделей — это простые граничные условия. Они задаются введением в модель клеток, поведение которых отличается от описанного в первой главе. Определены следующие типы клеток: рабочие клетки, клетки-стенки и клетки-источники частиц.
Элементарные автоматы каждого из этих типов имеют различные таблицы переходов. Фаза сдвига в клетках всех типов происходит одинаково, различия поведения есть только на фазе столкновения. Поведение рабочих клеток описано в п. 2.2, В клетках стенки частицы «отражаются» по различным правилам, определяющим свойства стенки, нарушая при этом закон сохранения импульса. При допущении равенства нулю скорости потока вблизи стенок, правилом отражения клеток-стенок является изменение всех векторов скорости частиц на противоположные (рис. 2.4а). Именно такие стенки были использованы в экспериментах, описанных в четвертой главе диссертации. Таким образом, стенки и источники определяют граничные условия автомата. Из правил поведения автомата следует, что объем вычислений при моделировании зависит только от общего количества клеток в автомате, но не зависит от их типов. Поэтому сколь сложными бы ни были граничные условия, время моделирования от этого не увеличится.
Функционирование клеток стенки определяется способом отражения частиц от препятствий. Один из них, как это было сказано выше, — отражение в обратном направлении (рис. 2.4а). При другом, более сложном способе частицы отражаются по законам геометрической оптики (рис. 2.46). Здесь необходимо учитывать не только состояние рассматриваемой клетки, но также и положение соседних клеток препятствия. Третий способ — частица отражается в любом допустимом направлении (туда, где нет препятствия) равновероятно (рис. 2.4в). Четвертый (интегрированный) способ состоит из этих трех способов отражения с заданной вероятностью срабатывания каждого из них при каждом конкретном столкновении.
С точки зрения практического применения наибольший интерес представляет первый способ отражения. При осреднении скоростей частиц на расстоянии до стенки d = 0 осредненные значения скорости равны нулю, что соответствует условию равенства нулю скорости потока жидкости возле стенки [83, 124]. Простейшим примером отражения частиц от стенки служит отражение одиночной частицы. Правила столкновения для некоторых таких случаев отражения в обратном направлении приведены в таблице 2.4. В левом столбце таблицы даны вектора состояний клеток, содержащих одну частицу, до столкновения, в правом столбце — после столкновения.
При одновременном столкновении двух или более частиц со стенкой результатом столкновения Si(s) является поразрядная дизъюнкция соответствующих строк таблицы 2.4.
Функционирование клеток-источников заключается в генерации частиц со всевозможными направлениями вектора скорости с определенной условиями эксперимента вероятностью/?. Это означает, что вектор состояния клетки источника с вероятностью р изменяется на s = 111111111110. Выстроив клетки-источники в пространстве в одной плоскости, можно получить источник равномерного потока частиц. Изменяя вероятность рождения частиц, можно варьировать концентрацию частиц в потоке.
Как уже отмечалось в п. 1.2, параметрам моделируемой жидкости соответствуют не микроскопические параметры модели, а ее макроскопические (осредненные) значения, причем осредненная скорость модельных частиц соответствует скорости потока, а осредненная концентрация частиц соответствует давлению. Существует два способа осреднения: по пространству и по времени. Можно также применять третий способ, комбинированный по пространству и времени.
Окрестностью осреднения с центром в клетке wo є W является множество клеток автомата, удовлетворяющее условию Av (w0, г) = {w: (wx - w0xf + (wy - wQyf + (w2 - w0z)2 r2}, (2.7) где wx, wy, wz и Wo , woy, w0z - координаты клеток w и wo соответственно. Таким образом, окрестность осреднения составляют все клетки, находящиеся внутри шара радиуса г с центром в клетке WQ. При осреднении по пространству результатом является сумма векторов скорости частиц (и) во всех клетках, попадающих в окрестность осреднения Av(wQ) клетки w0 в определенный момент времени /, деленная на количество клеток в окрестности: 1 Ь" (" = П-Т S 2Х (0, (2-8) где Сц вектор скорости 1-й частицы в клетке Wk є AV{WQ), bm — количество соседних клеток, a \Av\ — количество клеток в окрестности осреднения.
Важным параметром является радиус осреднения г. Чем больше г, тем меньше влияние автоматного шума на результат. Автоматный шум обусловлен дискретностью модели. Результирующий вектор скорости (и) принимает значение из конечного множества сумм дискретных векторов c/jftj находящихся во всех возможных комбинациях. Следует заметить, что при осреднении по пространству при достаточно больших размерах окрестности получить осредненные значения вблизи границ невозможно.
Параллельная реализация симулятора
Симулятор трехмерного потока на большинстве задач по моделированию потоков жидкости требует намного больше вычислительных ресурсов, чем любой другой модуль программного комплекса. Время решения многих задач с использованием последовательной версии симулятора очень велико, поэтому возникает необходимость в его параллельной версии. Такая версия была создана с использованием библиотеки параллельного программирования MPI [125, 126].
Параллелизм симулятора реализован следующим образом. Массив, хранящий конфигурацию КА, разрезается на слои вдоль координаты X на количество частей, равное количеству процессоров (рис. 3.2). Память под них на процессорах выделяется с перекрытием, изображенным клетками с горизонтальной штриховкой, для хранения состояний одного слоя соседних клеток. Количество слоев, хранящееся на каждом процессоре, не может быть меньше трех, т.к. внутренние слои, как будет показано ниже, обрабатываются отдельно от граничных. Поэтому количество процессоров, на котором происходит распараллеливание, не может превышать nJ3, где пх - размер К А вдоль оси Ох. Так как линейные размеры автоматов в большинстве задач составляют тысячи или даже десятки тысяч клеток, то моделирование можно производить на вычислительных системах с сотнями и даже тысячами процессоров. Следует также заметить, что размеры автомата не обязаны быть кратными количеству используемых процессоров.
После выделения памяти ведущий (нулевой) процессор рассылает по остальным процессорам требуемые части массива конфигурации КА. MPI_Send(part,partsize,MPI_BYTE,r,О, MPI_COMM_WORLD) ; (3.9) Остальные процессоры соответственно принимают эти части, чтобы в последствии их обрабатывать. MPI_Recv са,cas і ze,MPI_BYTE,0,0, MPI_COMM_WORLD,&status); (ЗЛО)
Также ведущий процессор рассылает остальным процессорам правила столкновения, количество тактов и другие глобальные параметры.
На каждом такте в фазе столкновения все элементы массива конфигурации КА, распределенные процессору, обрабатываются им последовательно в соответствии с (3.4), (3.5) и (3.6), при этом не требуются значения соседних элементов.
В фазе сдвига вначале инициируется пересылка граничных слоев процессором с номером rank левому (процессору с номером rank-І) и правому (процессору с номером гапк+1) соседям (стрелки на рис. 3.2). Крайние процессоры (с нулевым и максимальным номерами) не обмениваются с левым и правым соседями соответственно. Для того чтобы избежать тупиковой ситуации, когда все процессоры ожидают информацию от соседа, обмен происходит по следующему правилу. Процессоры с четными номерами обмениваются вначале с левым соседом, а затем с правым, а процессоры с нечетными номерами, наоборот, вначале с правым, а затем с левым. Функции обмена приведены ниже. MPI_Sendrecv_replace{boundleft,boundsize,MPI_BYTE, rank-1,0,rank-1,0,MPI_COMM_WORLD, status); MPI_Sendrecv_replace(boundright,boundsize,MPI_BYTE, rank+l,0,rank+l MPI_COm_WORLD,&status) ; (3.11)
После инициализации этих пересылок запускается сдвиг во внутренних слоях массива конфигурации (клетки белого цвета), где не требуется информация от соседних процессоров. Вычисления производятся в соответствии с (3.7) и (3.8). В то же самое время идет обмен граничными элементами массива с соседями (3.11). После того, как сдвиг во внутренних слоях выполнен и пересылки завершены, производится сдвиг в граничных клетках (с наклонной штриховкой), при этом используются принятые от соседних процессоров данные. Таким образом, обмен состояниями клеток в граничных слоях с соседними процессорами происходит на фоне вычисления сдвига во внутренних слоях массива конфигурации каждым процессором независимо, что существенно уменьшает потери времени при обработке.
После выполнения заданного в эксперименте количества итераций, ведущий процессор собирает у остальных процессоров их части конфигурации, аналогично (3.9) и (3.10). Новая конфигурация КА снова оказывается на одном процессоре в одном массиве. Далее она может быть сохранена в файл или использована для получения осредненных значений скорости и концентрации в зависимости от условий эксперимента. Таким образом, с точки зрения пользователя функционирование программного комплекса с параллельной версией симулятора ничем, кроме меньшего времени вычислений, не отличается от функционирования комплекса с последовательной версией.
Поток вязкой жидкости в трубе. Сравнение с законом Пуазейля
При моделировании потоков вязкой жидкости, осуществляемом при помощи программного комплекса, описанного в главе 3, заданными параметрами являются геометрические размеры объектов, плотность и вязкость жидкости, давление на границах, а также радиус осреднения. Осредненные значения скорости и концентрации частиц можно вычислять, как это было описано в п.2.4, используя осреднение по времени или по пространству. Осреднение по времени позволяет работать вблизи границ моделируемых объектов, но имеет существенный недостаток, занимает очень много машинного времени. Осреднение по пространству выполняется достаточно быстро, и поэтому в тех задачах, где можно обойтись без вычисления осредненных значений скорости и концентрации вблизи стенок его использование предпочтительнее. Искомыми величинами являются скорость потока и давление в каждой точке внутри объекта.
При моделировании потоков жидкости с осреднением по пространству можно выделить три условия. 1. Изменение состояния во времени должно быть очень малым. То есть, для любой клетки iv є должно найтись такое количество итераций Т, после которого начинают выполняться условия: \uw(t)-uw(t+l)\ e, (4.1) где t Т — произвольный такт КА, uw - осредненная скорость потока в клетке wt nw - осредненная концентрация частиц в клетке w, а є О -наперед заданная допустимая погрешность вычислений. Это означает, что скорость потока и концентрация частиц в каждой точке должны оставаться в пределах допустимой погрешности на протяжении некоторого заданного количества итераций.
Размер КА должен быть достаточно большим, чтобы использовать радиус осреднения, обеспечивающий уровень автоматного шума, не превышающий заданную величину єает. Шум возникает вследствие дискретности автомата, поэтому он тем меньше, чем больше дискретных частиц участвуют в осреднении, т.е. чем больше радиус осреднения. С другой стороны, для получения результатов осреднения, расстояние d между противоположными стенками объекта, в котором производится осреднение с радиусом г клеток, должно составлять не менее 2г клеток: d 2r. (4.2)
Если в расстояние между противоположными стенками какого-либо объекта не будет укладываться два радиуса осреднения, то не удастся получить осредненные значения скорости потока и концентрации ни в одной из точек этого объекта. Таким образом, выбор размеров автомата подразумевает оценку размеров самых мелких деталей моделируемого объекта и соизмерение ее с радиусом осреднения.
3. Расстояние от клеток-источников и открытых краев автомата должно быть достаточным для исключения влияния краевых эффектов. Вблизи клеток-источников и открытых краев образуется избыток и недостаток частиц соответственно. Как было показано во второй главе в соотношении (2.11), модельная вязкость vMOd остается неизменной только в определенном диапазоне концентраций частиц 3 vMOd 10. Для того чтобы концентрация не выходила из этого диапазона, расстояние 1гр от источника или открытого края до области осреднения должно превышать расстояние d между противоположными стенками объекта, расположенными перпендикулярно линиям тока: hp d. (4.3) Это дополнительное расстояние также необходимо учитывать при выборе размеров автомата.
При использовании программного комплекса, описанного в главе 3, для моделирования стационарных потоков с осреднением по пространству в алгоритме моделирования имеются некоторые особенности выбора размеров КА и задания необходимого числа итераций.
1. Выбор размеров и конструирование начальной конфигурации КА.
При конструировании КА вначале оценивается его размер с учетом размера моделируемого объекта и его сложности, а также имеющейся вычислительной мощности. Размер КА не должен быть слишком большим, чтобы моделирование могло выполниться за приемлемое время, но в то же время он не должен быть слишком маленьким, чтобы при соблюдении условия (4.2) можно было использовать достаточно большой радиус осреднения. Из размеров К А и размеров моделируемого объекта согласно , / (2.17) вычисляется коэффициент // - . Конструирование начальной мод конфигурации автомата производится с использованием конструктора граничных условий в соответствии с рекомендациями, данными в п. 3.6.
2. Задание такого количества итераций Т, при котором процесс станет стационарным.
Заранее точно определить количество итераций, необходимое для достижения условий стационарности (4.1), для каждого конкретного случая невозможно. Оно зависит от многих параметров, основными из которых являются размеры К А, сложность границ и начальное распределения частиц в клетках. Общей рекомендацией служит использование числа итераций Г, на порядок превосходящее длину самой длинной линии тока I (Т 10 I), но это лишь нижняя эмпирическая оценка. Очевидно, что чем больше размер КА, тем больше нужно итераций для того, чтобы частица, выпущенная источником на одном конце автомата и не встретившая препятствий на своем пути, достигла его противоположного конца, т. к. за одну итерацию она перемещается на одну клетку. Если же она сталкивается со стенкой, то она начинает движение в противоположном направлении и движется до тех пор, пока не столкнется за выпущенной за ней вслед другой частицей.
