Моделирование потоков вязких жидкостей с использованием систем клеточных автоматов

Дипломная работа Моделирование потоков вязких жидкостей с использованием систем клеточных автоматов

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

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

Определения

В настоящем отчете применяют следующие термины с соответствующими определениями:

ГИДРОДИНАМИКА — раздел физики сплошных сред, изучающий движение идеальных и реальных жидкости и газа.

ЛАМИНАРНОЕ ТЕЧЕНИЕ — течение, при котором жидкость или газ перемещается слоями без перемешивания и пульсаций (то есть беспорядочных быстрых изменений скорости и давления).

ТУРБУЛЕНТНОЕ ТЕЧЕНИЕ — явление, заключающееся в том, что при увеличении скорости течения жидкости или газа в среде самопроизвольно образуются многочисленные нелинейные фрактальные волны и обычные, линейные различных размеров, без наличия внешних, случайных, возмущающих среду сил и/или при их присутствии.

КЛЕТОЧНЫЙ АВТОМАТ — набор клеток, образующих некоторую периодическую решетку с заданными правилами перехода, определяющими состояние клетки в следующий момент времени через состояние клеток, находящимися от нее на расстоянии не больше некоторого, в текущий момент времени.

Обозначения и сокращения

КА — клеточный автомат

LGCA — клеточный автомат по модели решетчатого газа (LatticeGasCellularAutomata)

HPP — частный случай модели LGCA, предложенный Харди, Пацисом и Помо (Hardy, dePazzis, Pomeau) в 1973 году

FHP — частный случай модели LGCA, предложенный Фришем, Хасслаэром и Помо (Frisch, Hasslacher, Pomeau) в 1986 году

BGK — оператор столкновений, предложенный Бхатнагаром, Гроссом и Круком (Bhatnagar, Gross, Croock)

HSV- (англ. Hue, Saturation, Value — тон, насыщенность, значение) — цветовая модель, в которой координатами цвета являются цветовой тон, насыщенность и яркость.

Моделирование движения жидкостей и газов является очень важной прикладной задачей при исследовании процессов в промышленных аппаратах. Эта задача относится к механике сплошных сред и привлекает пристальной внимание исследователей.

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

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

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

В основе такого моделирование лежит концепция клеточных автоматов. На ее основе было создано несколько различных подходов к моделированию жидкости. Однако наиболее современным подходом является моделирование, основанное на объединении концепции клеточных автоматов и статистической физики, в частности, с моделью решетчатого газа Больцмана.

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

1. Формирование требований

1.1 Роль гидродинамических процессов в современной технике и технологиях

Еще задолго до становления гидродинамики как науки гидродинамические процессы активно исследовались и использовались. Люди строили сложные системы вентиляции, трубопроводы и фонтаны, активно развивалось мореплавание.

В современном мире роль гидродинамики трудно переоценить. Ее приложения затрагивают самые разнообразные сферы человеческой деятельности. Исследование океанических течений и атмосферы применяется при проектировании летательных аппаратов и водных средств, а также при решении практически всех задач метеорологии. В исследования природных ресурсов гидродинамика применяется в частности при изучении фильтрации грунтовых вод и нефти в подземных месторождениях. Гидродинамические процессы необходимо учитывать при расчете трубопроводов, гидротурбин и насосов. Исследование движения малых объемов жидкости (микрогидродинамика) применяется, например, в струйной печати.

Задачи гидродинамики тесно связаны с другими прикладными разделами физики. Одной из таких смежных задач является задача тепломассообмена, которая возникает, например, при проектировании трубопроводов, предназначенных для передачи горячей жидкости. В данном случае на классическую задачу движения жидкости накладывается задача распространения тепла.

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

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

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

1.2 Необходимость использования компьютерных методов при моделировании гидродинамических процессов

Усложнение математических моделей привело к тому, что человек перестал справляться с расчетами. Ответом на это явилось создание компьютера. Современный компьютер позволяет решать сложные задачи из самых различных областей, в том числе и гидродинамики. Численные расчёты востребованы как в фундаментальных науках, так и в прикладных. Численные расчёты постепенно вытесняют из употребления эмпирические методы.

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

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

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

— большие объемы данных;

— требование достаточно быстрого получения результата;

— недопустимость ошибок при вычислениях.

2. Разработка концепции

2.1 Классические континуальные модели и представления

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

Одно из первых описаний движения сплошной среды было получено в 1755 году Эйлером:

(1)

где — плотность жидкости, — вектор скорости, — время, — давление жидкости.

Это уравнение описывает движение невязкой, несжимаемой жидкости. Такую модель идеальной сплошной среды фон Нейман назвал «моделью сухой воды».

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

Для понимания процессов рождения и диссипации вихрей, а также процессов взаимодействия (трения) жидкости со стенками сосуда, понадобилось перейти к более сложной модели — модели «мокрой воды», учитывающей влияние вязкости жидкости.

Впервые уравнения движения вязкой жидкости выписал французский учёный и инженер Анри Навье. Для этого потребовалось ввести тензор напряжений, то есть учесть не только нормальные силы (давление), но и касательные силы. В правую часть уравнения (1) Навье ввёл дополнительный член, ответственный за проявление вязкости.

(2)

где — коэффициент кинематической вязкости.

Большой вклад в исследование этого уравнения внёс Стокс. Поэтому уравнения (2), а также их обобщения на случай движения жидкостей с другими свойствами называются уравнениями Навье-Стокса. Уравнения Эйлера — частный случай уравнений Навье-Стокса при .

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

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

Общий критерий возникновения турбулентности установлен Осборном Рейнольдсом в 1883 году. Он подкрашивал ламинарную струйку потока во входной части стеклянной трубки и следил, когда течение станет турбулентным, фиксируя при этом критическое значение безразмерного определяющего параметра, названного впоследствии в его честь числом Рейнольдса.

Представив скорость в виде суммы средней и пульсационной составляющих, а давление в виде, Рейнольдс получил уравнения для средних величин, носящие его имя:

. (3)

Здесь подразумевается суммирование по индексу j. По сравнению с уравнением Навье-Стокса (2) это уравнение включает дополнительные напряжения — так называемые напряжения Рейнольдса. Попытки найти их вид из первых принципов физики оказались безуспешными, поэтому уравнение (3) стало базой для развития эмпирических теорий.

2.2 Обзор дискретных моделей решетчатых газов (LGCA)

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

2.3 Модель HPP

Первой методикой, базирующейся на применении решеточного газа, была модель, предложенная Харди, Пацисом и Помо в 1973 году (модель HPP).

Представим квадратную решетку с четырьмя ячейками в каждом узле, где такая ячейка связана по линиям со всеми ближайшими соседними узлами. Эти ячейки могут быть пусты или заняты максимум одной частицей с единичной массой m=1. Таким образом, каждая ячейка имеет только два возможных состояния (рисунок 1). Скорость, а поэтому и импульс, может быть представлена в каждой частице как вектор, связывающий узел с соседним по линии, на которой расположена частица. Эти вектора называются скоростями решетки. Микроскопическое взаимодействие строго локализовано и включает только частицы в одном конкретном узле. Частицы обмениваются импульсом, сохраняя сумму массы и импульса в каждом узле. После такого столкновения каждая частица перемещается по связанной с ней линии к соседнему узлу.

Вся микродинамика описывается такими столкновениями и перемещениями. Макродинамические значения массы и плотности импульса рассчитываются путем усреднения (вычисляются средние значения на больших областях пространства от сотен до тысяч узлов). Однако для описанной выше модели макроскопические параметры не подчиняются уравнениям Навье-Стокса.

Рисунок 1 — Столкновения в модели HPP

2.4 Модель FHP

Прошло более десяти лет до того момента, как Фриш, Хасслаэр и Помо (1986) нашли третье существенное условие помимо сохранения массы и импульса: решетка должна обладать определенной симметрией. В двумерном пространстве, например, четырехсвязная симметрия (квадратная решетка) не достаточна, в отличие от гексагональной симметрии. Так появилась модель FHP.

Классическая модель FHP обладает следующими свойствами:

— частицы перемещаются по однородной гексагональной решетке, в которой каждая частица связана с шестью соседними;

— векторы, соединяющие центры соседних ячеек называются скоростными векторами, и могут быть описаны следующим образом:

(4)

— каждая ячейка содержит шесть скоростных каналов, каждому из которых соответствует скоростной вектор;

— в каждом скоростном канале каждой ячейки не может находиться более одной частицы;

— все частицы обладают единичной массой;

— изменение состояние решетки во времени происходит путем чередования фаз столкновений и перемещений;

— в столкновении участвуют только частицы, находящиеся в одной ячейке.

Во время фазы перемещений все частицы передвигаются по скоростным каналам в соответствующие соседние ячейки. В фазу столкновений частицы в ячейках перераспределяются по скоростным каналам в соответствии с правилами, представленными на рисунке 2. Если распределение частиц по каналам ячейки отличается от показанного на рисунке, состояние ячейки не изменяется.

Существуют расширения модели FHP, известные как FHP-II и FHP-III. В FHP-II вводится дополнительный скоростной канал, соответствующий частицам, находящимся в покое и несколько новых правил. В модели FHP-III количество различных столкновений доведено до максимума.

В целом модели FHP соответствуют равнениям Навье-Стокса, однако они обладают рядом недостатков, о которых будет сказано ниже.

Рисунок 2 — Cтолкновения в модели FHP

2.5 Сравнительный анализ классического и дискретного подходов

Традиционное моделирование потоков жидкости (и других физических процессов) в целом начинается с нелинейных дифференциальных уравнений в частных производных. Затем эти дифференциальные уравнения дискретизируют конечными разностями. Полученное алгебраическое уравнение или система обыкновенных дифференциальных уравнений решается стандартными численными методами. Такой подход кажется понятным и привычным, однако он обладает некоторыми трудностями. Во-первых, при численном решении дифференциальных уравнений в частных производных при разложении в ряд Тэйлора возникают ошибки усечения при переходе к конечным разностям. В то же время возникает вопрос о сохранении определенных величин при переходе к дискретизированным формам уравнений. Последнее является наиболее важным при интегрировании по большим промежуткам времени в таких прикладных областях как, например, моделирование океана или парное моделирование океана и атмосферы. Небольшая «утечка» может за определенное время опустошить океан. Вторая проблема такого типа численных методов — неустойчивость численного решения.

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

Существенным недостатком этих моделей является их ресурсоемкость. Так как происходит симуляция поведения каждой частицы, нетрудно представить, какой объем вычислительных ресурсов требуется для моделирования потока жидкости, например, в трубопроводе, не говоря уже о моделировании океанических течений.

Также простота правил перемещения и столкновения частиц, необходимая для возможности моделирования более-менее больших объемов жидкости, вызывает определенные шумы. И для получения соответствия законам Навье-Стокса требуется проводить значительные усреднения.

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

2.6 Уравнение Больцмана и распределение Максвелла

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

Рассмотрим вероятностное распределение частиц в микроскопическом объеме, определенное таким образом, что обозначает количество частиц в момент времени t, находящихся в параллелепипеде, ограниченном точками r и r+dr, скорости которых лежат в интервале между c и dc. Предположим, что на все частицы действует сила F. Тогда, если частицы вещества не сталкиваются между собой, то изменение распределения f со временем должно происходить на основании уравнения:

(5)

где m — масса частицы, а dt — малый интервал времени.

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

(6)

Оператор столкновений должен быть определен так, чтобы выполнялись законы сохранения массы, импульса и энергии. На практике широко используется оператор BGK (Bhatnagar — Gross — Croock)

(7)

где — время релаксации, которое имеет физический смысл времени столкновения между молекулами или времени перераспределения энергии между внутренними колебаниями степеней свободы молекулы [4], — плотность равновесного распределения частиц на данном отрезке, соответствующая уравнению Максвелла-Больцмана (8).

(8)

где k_B — константа Больцмана, T — температура.

2.7 Модель решетчатого газа Больцмана (LBM)

Продолжением и развитием дискретного подхода стала модель решетчатого газа Больцмана. Однако в основу этой модели положены совершенно иные представления о решетке и ячейках, поэтому LBM считается самостоятельным численным методом моделирования гидродинамики.

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

Чтобы преодолеть эту трудность понятие конкретной частицы заменили плотностью распределения частиц, а вместо итерационных правил столкновения описываются оператором столкновений.

В LBM предусмотрено наличие девяти направлений (скоростных каналов). Это сама клетка и восемь ее соседей. Такая решетка получила название D2Q9 (2 измерения, 9 соседей). На первый взгляд создается впечатление, что такая решетка не обладает симметрией и должна приводить к тем же трудностям, что и модель HPP. Однако, как уже было сказано, по решетке рассматривается не непосредственное движение частицы, а их распределение. В этом случае симметрия достигается путем введения определенных коэффициентов в плотность распределения.

Рисунок 3 — Скоростные каналы LBM

Вектора имеют единичную длину и задают направления скоростных каналов:

(9)

Локальная плотность потока для ячейки вычисляется как сумма значений плотности распределения по всем скоростным каналам:

. (10)

Сумма произведений плотности распределения потока на векторы скоростных каналов составляет плотность импульса. Разделив плотность импульса на плотность частиц, мы получим вектор локальной скорости для ячейки потока:

. (11)

Итерацию клеточного автомата, функционирующего по методу решетчатого газа Больцмана, можно записать следующим образом:

(12)

Это выражение является дискретным аналогом уравнения Больцмана, что было доказано Стерлингом и Ченом в 1996 году.

Перемещение функции распределения по скоростным каналам за шаг времени представлено на рисунке 4.

Рисунок 4 — Перемещение функции распределения Для оператора столкновений плотность распределения частиц по ячейкам решетки заменена плотностью распределения по скоростным каналам. Он вычисляется следующим образом:

(13)

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

(14)

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

(15)

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

(16)

2.8 Задание граничных условий для LBM

По словам Чена (1996), поиск граничных условий является задачей ничуть не менее важной и трудоемкой, чем поиск самого алгоритма вычислений. Трудность их определения связана в первую очередь с тем, что даже с точки зрения классической гидродинамики реальные граничные условия еще не были поняты до конца. Модель решетчатого газа Больцмана является известной методикой моделирования жидкости и газа, однако не существует единых правил задания граничных условий. Большинство авторы статей при описании граничных условий указываю только, что при использовании LBM граничные условия задаются очень легко, однако не описывают каким образом. Другие приводят несколько способов заданий граничных условий. По мнению этих авторов, выбор конкретного способа зависит от поставленной задачи и требует определенной интуиции. Поэтому в данной работе будет уделено особое внимание этому вопросу.

Следует сразу разделить граничные условия на два вида. Назовем первую группу граничных условий «тормозящими». Во-первых, к этой группе относятся границы потока и стенок, между которыми движется этот поток. Также при введении в поток препятствия, например щели, взаимодействие с ним будет происходить с учетом тех же граничных условий. Вторая группа — это «нетормозящие» граничные условия. Их рассмотрение связано с тем, что мы рассматриваем не весь поток, а только его часть. Таким образом, на одной границе необходимо прикладывать силу, вызывающую движение потока, а на другой обеспечивать «сток».

2.9 «Тормозящие» граничные условия

Наиболее простым алгоритмом, описывающим столкновение со стенкой, является алгоритм полного отражения. Его идея заключается в том, что сталкиваясь со стенкой, все частицы, задаваемые плотностью распределения, переходят в ту же ячейку, но в канал противоположного от стенки направления. При этом можно производить обратный отскок частицы под углом 180 градусов, либо отражение (рисунок 5). Авторы, использующие такой поход, отмечают отличие результатов моделирования в пограничных слоях от получаемых аналитически для простейших потоков.

Рисунок 5 — Столкновение со стенкой Существует другой алгоритм — отражение «на полпути» от стенки. При таком подходе границы потока и стенки располагается между первым и вторым и, соответственно, предпоследним и последним рядами ячеек. При этом в крайних рядах не происходит никаких столкновений и перераспределений, а попавшие туда частицы принудительно возвращаются на следующем шаге назад в поток. При этом подходе при интеграции ячеек с целью получения макрохарактеристик потока граничный слой не учитывается.

2.10 «Нетормозящие» граничные условия

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

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

На практике подбор граничного условия определяется поставленной задачей.

2.11 Соответствие реальных величин параметрам модели

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

Представим, что поток обтекает препятствие, например, цилиндр. Очевидно, что поведение потока будет зависеть от скорости входного (невозмущенного) потока и от линейных размеров (диаметр L) преграды. Также жидкость характеризуется кинематической вязкостью. Эти три параметра (, L и) имеют единицы измерения [м/c], [м] и [м²/с]. Легко заметить, что из этих параметров можно составить ни что иное как одно безразмерное число, а именно число Рейнольдса:

. (17)

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

Допустим, мы хотим смоделировать движение потока воды в трубе, длина которого составляет 1 м. На вход вода подается со скоростью 0,001 м/с. Предположим, что вода находится в стандартных условиях, и поэтому ее плотность составляет 1000 кг/м³, а кинематическая вязкость равна 10^-6 м²/c.

Примем, что длина потока, то есть 1 м, соответствует 100 ячейкам, а за 1 секунду происходит 10 итераций. Тогда начальная модельная скорость потока будет равна 0.001*100/10=0,01. Это означает, что за одну итерацию частицы потока проходит расстояние в 0,01 ячейки. Кинематическая вязкость составит 0,001. Чтобы определить плотность как число частиц (молекул воды) в ячейке, необходимо плотность воды разделить на ее молярную массу и умножить на число Авогадро, разделив на число ячеек в 1 м. Мы получим (1000*6,022*10²³)/(0,018*100)=3,34*10²⁶.

Стоит отметить, что для данного потока можно было выбрать и другой набор модельных величин. Но в любом случае при их выборе должно соблюдаться два условия:

— Соотношение должно оставаться постоянным. Только при выполнении этого условия поведение модели будет соответствовать поведению моделируемого потока при заданных условиях.

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

Для перехода между реальными и модельными величинами введем два параметра:

mчисло ячеек, содержащееся в одном метре;

s — число итераций, проходящих за одну секунду.

Формулы перехода представлены в таблице 1, где переменные с индексом «р» обозначают реальные величины, а без индекса — модельные; - молярная масса жидкости; - число Авогадро. Все реальные величины должны быть выражены в СИ.

Таблица 1 — Формулы перехода между реальными и модельными величинами

3. Техническое задание

3.1 Назначение системы

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

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

3.2 Цели создания системы

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

Критерием выполнения поставленной задачи будет служить соответствие получаемых результатов общепринятым методам вычислений и существующим экспериментальным данным.

В случае достижения данной цели приложение в дальнейшем может быть использовано для выполнения расчетов реальных физических процессов.

3.3 Требования к системе

Требования к системе в целом

Требования к структуре и функционированию системы

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

Требования к функциям, выполняемым системой

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

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

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

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

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

Требования к видам обеспечения

Требования к математическому обеспечению

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

Требования к программному обеспечению

Для работы программы на компьютере должна быть установлена операционная система MS WindowsXPSP2 или старше.

Требования к лингвистическому обеспечению

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

1. Требования к техническому обеспечению системы

— Многоядерный процессор;

— не менее 512 Мб оперативной памяти;

— наличие видеокарты;

— наличие 1 Гб свободного места на жестком диске.

4. Технический проект

4.1 Общесистемные решения

Схема функциональной структуры

Главной функцией создаваемого приложения является моделирование потока жидкости. С точки зрения пользователя данная функция включает в себя три задачи: выбор параметром модели, непосредственный расчет по методу LBM и вывод результата в числовом виде, а также в наглядной графической форме.

4.2 Описание автоматизируемых функций

Задание параметров модели

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

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

— количество ячеек в одном метре;

— количество итераций в секунду.

Должны задаваться следующие параметры жидкости (в скобках указаны реальные — модельные единицы измерения):

— длина потока (м — яч.);

— ширина потока (м — яч.);

— продолжительность моделирования (с — итер.);

— кинематическая вязкость (м²/с — яч.²/итер.);

— плотность (кг/м³ — частиц/яч.);

— молярная масса (кг/моль);

— начальная скорость (м/c — яч./итер.)

Пользователь может вводить параметры жидкости как в реальных единицах измерения, так и в модельных единицах.

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

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

4.3 Расчет по методу LBM

После задания параметров потока по команде пользователя должен начинаться расчет по методу решетчатого газа Больцмана. Расчет должен прекращаться, когда достигнута заданная продолжительность моделирования или при выходе параметров модели за допустимые границы.

4.4 Визуализация результата

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

Для визуализации результата выбрана цветовая модель HSV. Выбор остановлен на этой модели, поскольку она позволяет одновременно визуализировать два аспекта такого сложного параметра, как вектор скорости: ее модуль и угол.

HSV (англ. Hue, Saturation, Value — тон, насыщенность, значение) — цветовая модель, в которой координатами цвета являются:

— Hue — цветовой тон, (например, красный, зелёный или сине-голубой). Варьируется в пределах 0—360°, однако иногда приводится к диапазону 0—100 или 0—1.

— Saturation — насыщенность. Варьируется в пределах 0—100 или 0—1. Чем больше этот параметр, тем «чище» цвет, поэтому этот параметр иногда называют чистотой цвета. А чем ближе этот параметр к нулю, тем ближе цвет к нейтральному серому.

— Value — значение цвета. Также задаётся в пределах 0—100 и 0—1.

Наглядно данную модель можно представить в виде конуса (рисунок 6).

Рисунок 6 — Коническое представление модели HSV

Для визуализации вектора скорости мы будем использовать только два параметра этой модели.

Параметр Hue удобно использовать для обозначения угла поворота вектора скорости, поскольку этот параметр также по своей сути является углом. Соответствие цвета и численного значения угла представлено на рисунке 7.

Рисунок 7 — Шкала оттенков HUE

Параметр Valueпримем равным единице.

Параметр Saturation будем использовать для обозначения модуля скорости, поскольку он, как и модельная скорость, принимает значения от 0 до 1. Скорость, равная 0, обозначается белым цветом. С увеличением скорости, увеличивается и яркость.

Помимо визуального представления численные значения модуля и угла скорости каждой ячейки на каждой итерации должны сохраняться в текстовом файле, которые пользователь может просмотреть по окончании моделирования.

4.5 Программа и методика испытаний

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

В ходе этой проверки необходимо установить следующее:

— какое соотношение реальных параметров и параметров модели дает наиболее точный результат моделирования;

— на каких интервалах значений реальных физических величин модель решетчатого газа Больцмана пригодна для использования в качестве альтернативы реальному физическому эксперименту.

4.6 Сопоставление с аналитическим решением

Как уже упоминалось ранее, аналитического решения уравнения Навье-Стокса для общего случая не существует. Известно всего несколько аналитических решений уравнения Навье-Стокса для простейших случаев. Одно из них — плоский поток Паузейля. Это двумерный поток вязкой несжимаемой жидкости, которая движется с постоянной силой между двумя бесконечными плоскостями при постоянных значениях плотности и давления. В этом случае аналитическим решением уравнения Навье-Стокса является парабола (рисунок 8), которую можно задать формулой (18):

(18)

Рисунок 8 — Поток Паузейля В ходе испытаний необходимо выявить, соответствует ли профиль скоростей параболической форме.

4.7 Сопоставление с физическим экспериментом

С эмпирической точки зрения возникновение турбулентности в потоке характеризуют при помощи числа Рейнольдса. Считается, что потоки с небольшим числом Рейнольдса (Re<<1) являются ламинарными, вихревые дорожки Кармана наблюдаются при средних значениях (Re), и турбулентные потоки наблюдаются при очень больших числах (Re>>1).

Для модели число Рейнольдса рассчитывается по формуле (17). При соответствующих числах Рейнольдса в модели должны наблюдаться такие явления как вихревые дорожки за препятствиями и возникновение турбулентности.

4.8 Решения по программному обеспечению

Структура программного обеспечения

1. Структура расчетного модуля

Структура основной части приложения, выполняющей расчет по методу LBM, представлена в виде диаграммы классов на рисунке 9.

Рисунок 9 — Структура расчетной части приложения Класс «Клеточный автомат» хранит информацию о размерах потока и о таком параметре жидкости как время релаксации. Работа клеточного автомата будет состоять из последовательности операций «сделать шаг». Интерфейс класса «Клеточный автомат» приведен в приложении A.

Клеточный автомат включает в себя большое количество ячеек, которые в свою очередь характеризуются плотностью и скоростью. Для осуществления расчета скорость удобно хранить в виде двух векторов. Однако для удобства последующей визуализации также рассчитываются и хранятся угол и модуль скорости. Интерфейс класса «Ячейка» приведен в приложении Б.

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

Класс «Ячейка» является абстрактным. Клеточный автомат состоит из ячеек потока, граничных ячеек, истоков и стоков. Каждый из этих видов ячеек характеризуется своими способами расчета плотности, скорости, а также плотности распределения.

2. Общая структура приложения

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

Рисунок 10 — Общая структура приложения Модуль Mainобеспечивает интерфейс пользователя и отвечает за главную форму приложения. С помощью этой формы пользователь может вводить исходные данные, и на нее же выводится результат моделирования в графической форме.

Модуль Mainнепосредственно взаимодействует с классом CA (Клеточный автомат) расчетного блока, который уже в свою очередь взаимодействует с остальными классами вычислительного блока LBM.

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

Модуль Vis служит для визуализации результатов моделирования. Визуализация происходит посредством создания буферного кадра. Вывод кадра на форму происходит по отсчету системного таймера. Также в число функций данного модуля входит перевод значений скорости, заданных в модели HSV в модель RGB.

Методы и средства разработки программного обеспечения

При разработке программного обеспечения используется объектно-ориентированный подход. Для написания программного кода используется среда разработки Delphi 7. Для вывода графики используется метод попиксельного отображения путем построения буферного кадра.

Решения по математическому обеспечению

Алгоритм работы клеточного автомата

Алгоритм работы клеточного автомата представлен в виде диаграммы деятельности на рисунке 11.

Рисунок 11 — Алгоритм работы клеточного автомата гидродинамический процесс компьютерный моделирование На этапе создания потока происходит создание объекта класса «Клеточный автомат», который в свою очередь запускает конструкторы различных ячеек, тип которых определяется заданной пользователем схемой.

Также на этом этапе для каждой ячейки задаются начальные плотности и скорости.

После того, как клеточный автомат создан, рассчитывается время релаксации на основании заданной пользователем вязкости по формуле, получаемой из формулы (14):

(18)

После этого производится расчет начального равновесного распределения по формуле (15), а текущее распределение принимается равным начальному. Далее циклически выполняется заданное число итераций клеточного автомата.

Каждая итерация начинается с расчета распределения частиц по скоростным каналам для каждой ячейки по формуле, получаемой из формул (12) и (13):

(19)

Данная формула применяется для ячеек потока. В граничных ячейках вместо расчета просто производится смена направления распределения по формуле (20):

(20)

Далее для ячеек потока плотность и скорость рассчитываются по формулам (10) и (11) соответственно. Для остальных типов ячеек перерасчет плотности и скорости не производится. Далее по формуле (15) рассчитывается равновесное распределение. Расчет модуля и угла скорости Поскольку расчет по методу LBMиспользует скорость, заданную в виде суммы векторов, а отображать необходимо модуль скорости и ее направление, то возникает необходимость перехода между этими двумя способами задания вектора скорости.

Расчет модуля скорости производится по формуле (21):

. (21)

Расчет угла производится по следующему правилу:

(22)

Алгоритм перевода из HSVв RGB

Как угол скорости, так и переменная H, изменяются в интервале [0.360).Поэтому приравниваем переменную H углу скорости:

(23)

Модуль скорости и переменная Sизменяются в интервале [0.1]. Однако если приравнять переменную Sк модулю скорости при низких скоростях яркость является недостаточной, чтобы пользователь мог разглядеть направление скорости. Поэтому зависимость яркости от модуля скорости будем рассчитывать по следующей формуле:

(24)

График зависимости яркости от модуля скорости приведен на рисунке 12.

Рисунок 12 — Зависимость яркости от модуля скорости

5. Рабочая документация

5.1 Подготовка к работе

Перед запуском приложения необходимо скопировать на свой персональный компьютер и поместить в одну папку файл CA. exe, а также папки Picи Save.

Для начала работы с приложением нужно запустить файл CA.exe.

5.2 Описание операций

Установка параметров

После запуска приложения перед вами откроется главная форма приложения.

Сначала необходимо установить параметры жидкости и клеточного автомата. Для этого выберем пункт меню «Параметры».

Выберем последний подпункт, «Соотношение реальных и модельных величин…».

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

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

Выберем пункт меню «Размеры потока…».

Здесь пользователь может задать размеры моделируемого потока в метрах или размеры клеточного автомата в ячейках.

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

Выберем пункт меню «Продолжительность моделирования…». На этой форме можно моделируемое время жизни потока в итерациях или секундах.

Отметим, что задается отрезок времени реального потока, а не время, в течение которого пользователь будет наблюдать процесс моделирования на экране.

Выберем пункт меню «Параметры жидкости…».

Здесь можно задать кинематическую вязкость и плотность моделируемой жидкости в реальных или модельных единицах измерения. Также можно задать молярную массу жидкости. Этот параметр используется для перехода от реальных к модельным величинам. Однако в модели молярная масса жидкости всегда принимается раной единице.

Последний пункт меню — «Начальная скорость…». На данной форме можно задать начальную скорость потока в реальных или модельных величинах. Если флаг «Во всем потоке установлен», то данная скорость устанавливается в каждой ячейке потока, а также ячейках типа «Исток» и «Сток». В противном случае данная скорость устанавливается только ячейках типа «Исток». Остальные ячейки потока считаются неподвижными.

5.3 Задание формы потока

По умолчанию поток имеет прямоугольную форму. На левой границе располагаются ячейки типа «Исток». Направой — ячейки типа «Сток». Верхняя и нижняя стенки — «Граничные ячейки». Препятствий в потоке нет.

Пользователь может поместить препятствие внутрь потока. Для этого нужно установить курсор внутри формы, нажать левую кнопку мыши, переместить курсор и отпустить левую кнопку мыши. Таким образом, задаются два противолежащих угла прямоугольника, изображающего препятствие. Тип ячеек меняются на «Граничные ячейки». И они закрашиваются черным цветом.

Данный рисунок можно сохранить в файл с целью последующего многократного использования. Для этого выберем пункт меню «Шаблоны>>Сохранить». Откроется стандартное меню сохранения файла.

По умолчанию рисуется граничная ячейка. Если мы хотим нарисовать ячейку потока, выберем пункт меню «Рисование >>Простая ячейка». Аналогичным описанному способом нарисуем в уже созданном препятствии щель.

Чтобы снова рисовать граничные ячейки, нужно выбрать пункт меню «Рисование >> Граничная ячейка».

При нажатии пункта «Обновить» все препятствия в потоке исчезают и поток принимает исходный вид (рисунок 14). Однако все параметры потока при этом сохраняются.

Выбрав пункт меню «Шаблоны>>Загрузить», можно загрузить созданный ранее рисунок. При этом размеры потока примут значения, которые они имели при создании рисунка. Поэтому один и тот же рисунок можно использовать только для одного и того же размера потока в ячейках.

5.4 Моделирование

Для запуска моделирования нужно нажать кнопку «Старт». При этом в левом верхнем углу экрана откроется шкала HUE (рисунок 8), а на главной форме приложения можно будет наблюдать течение жидкости.

При нажатии «Стоп» моделирование остановится, но все параметры, в том числе и форма препятствий сохранятся. Если после этого нажать «Старт», моделирование начнется сначала.

Результаты моделирования можно посмотреть в численном виде. Для этого необходимо открыть файл Iteration. txt, находящийся в папке Save. В этом файле для каждой итерации в табличной форме записаны значения модуля и направления скорости. Фрагмент этого файла приведен на рисунке 24.

Рисунок 24 — Фрагмент файла Iteration. txt

5.5 Аварийные ситуации

В ходе моделирования может возникнуть ошибка, представленная на рисунке 25.

Рисунок 25 — Ошибка превышения скорости Данная ошибка означает, что в одной из ячеек скорость превысила 1. Модель LBM предназначена только для моделирования потоков, в которых скорость не превышает одной ячейки за итерацию. Если это условие не соблюдается, то дальнейшее моделирование не имеет смысла.

В случае возникновения такой ошибки необходимо изменить параметры моделирования.

6. Экспериментальная часть