CC BY
72
2012 Дискретные модели реальных процессов №4(18)
УДК 519.68
РАЗВИТИЕ МЕТОДА НЕПРЕРЫВНЫХ АСИНХРОННЫХ КЛЕТОЧНЫХ АВТОМАТОВ ДЛЯ МОДЕЛИРОВАНИЯ ТУРБУЛЕНТНЫХ ПОТОКОВ
Г. Д. Тимчук, В. В. Жихаревич
Национальный технический университет «Харьковский политехнический институт», Черновицкий факультет, г. Черновцы, Украина
E-mail: [email protected], [email protected]
Представлен метод непрерывных асинхронных клеточных автоматов, моделирующих турбулентные потоки вещества. Предлагается его модификация, которая состоит в реализации вероятностного механизма передвижения клеточных автоматов вдоль компонент векторов скоростей. При этом моделирование процесса переноса вещества осуществляется путём направленного дублирования содержания соседних ячеек. Приведены результаты вычислительного эксперимента по наблюдению турбулентной динамики на примере потока вещества в трубе, которое имеет препятствия.
Ключевые слова: клеточный автомат, компьютерное моделирование, турбулентность, уравнения Навье — Стокса.
Введение
Проблема моделирования турбулентных потоков является актуальной уже более века. Общепринятым подходом в этом отношении считается численное решение уравнения Навье — Стокса и его различных модификаций. При этом возникают трудности, связанные с вычислительной сложностью, что, в свою очередь, стимулирует поиск альтернативных методов, адаптированных для компьютерной реализации.
Последнее время довольно распространёнными становятся клеточно-автоматные методы моделирования [1], которые привлекают своей простотой, универсальностью и естественным параллелизмом. В частности, в работе [2] продемонстрированы возможности метода непрерывных асинхронных клеточных автоматов по моделированию различных процессов, в том числе и волновой динамики, которая, как известно, является одним из решений уравнения Навье — Стокса [3]. При этом модель не предусматривает перемещения вещества, а наблюдение явления турбулентности невозможно.
Для устранения указанного недостатка и построения модели, которая даст возможность исследовать турбулентность, предложена модификация, суть которой состоит в реализации вероятностного механизма передвижения клеточных автоматов вдоль компонент векторов скоростей. Под передвижением в нашем случае следует понимать дублирование содержания соседних ячеек при осуществлении перемещения. Введение такого подхода определяется необходимостью моделирования сплошного непрерывного потока, а не потока отдельных дискретных частиц.
1. Постановка задачи
Как известно, турбулентное течение отличается от ламинарного вихревой динамикой, при которой имеет место вращательно-колебательный характер движения потоков вещества (рис. 1) [4]. Турбулентность может возникать в газах, жидкостях или сыпучих
средах за счёт существования сил внутреннего трения — динамической вязкости (п). Кроме вязкости, для образования турбулентного течения существенными параметрами являются скорость (и), плотность вещества (р) и характерные размеры потока (Ь). Соотношение рХи/п или же Ьи/и, где V = п/р — кинематическая вязкость, называют числом Рейнольдса (Ие) [5]. Превышение числом Рейнольдса некоторого критического значения Иекр определяет возможность возникновения турбулентной динамики. Величина Иекр зависит от конкретного вида течения.
Рис. 1. Пример турбулентного потока (вихревая дорожка Кармана)
В общем случае как ламинарные, так и турбулентные потоки вещества можно описать уравнениями Навье — Стокса, т. е. системой нелинейных дифференциальных уравнений в частных производных [6]:
^ + (7 -У)7 = -1 Ур + VД7 + 7(7,Ъ); (1)
дъ р
др + У (р1 )=0, (2)
где 7 (7,Ъ) —вектор скорости вещества; р(7,Ъ) —давление; V > 0 — кинематическая вязкость; р — удельная плотность; 7 (-,о — внешняя сила; У — оператор Гамильтона; Д — оператор Лапласа. Уравнение (2) отображает закон сохранения массы. Неизвестными величинами являются скорость 7 (7,ъ) и давление р(7, Ъ).
В рамках этой работы предлагается клеточно-автоматная аппроксимация уравнения Навье — Стокса. Таким образом, динамическая картина, полученная в результате клеточно-автоматных взаимодействий, будет альтернативой численному решению уравнения Навье — Стокса. Основная задача данной работы — демонстрация принтти-пиальной возможности построения модели турбулентных потоков на основе непрерывных вероятностных асинхронных клеточных автоматов.
2. Описание клеточно-автоматной модели
Опишем структуру клеточно-автоматного поля в случае моделирования двумерных потоков жидкости (рис. 2). Поле содержит три слоя: 1 — слой давления Р,2 — слой X -компоненты скорости Ух, 3 — слой У -компоненты скорости Уу. Ячейки поля принимают вещественные непрерывные значения.
Площадь ячейки определяется размерностью клеточно-автоматного поля в соотношении к геометрическим размерам фрагмента моделируемой системы.
Процесс моделирования представляет собой итерационный цикл клеточно-автоматных взаимодействий с использованием асинхронной схемы. Данная схема предусматривает циклическое выполнение трех типичных шагов:
Рис. 2. Структура клеточно-автоматного поля двумерной модели потока
1. На клеточно-автоматном поле случайным образом выбирается некоторая клетка (г = 1) с целочисленными координатами ж1, у1. При этом все клетки выбираются равновероятно.
2. Случайным равновероятным образом выбирается некоторая соседняя клетка (г = 2) с целочисленными координатами ж2, у2. В качестве схемы соседства принята окрестность Неймана, т. е. у клетки есть только четыре соседа (рис. 3).
х2=х* /=У+і
х2=хЧ у2=у1 X1,/ х2=х‘+1 /=/
•< Г Ч, ^ II ) Ч-
Рис. 3. Окрестность Неймана и координаты клеток
3. Происходит клеточно-автоматное взаимодействие между двумя клетками, суть которого заключается в модификации непрерывных значений соответствующих слоёв клеток согласно следующей системе уравнений:
( с*' = 4 + ВДж1 - 1Г2)(с2 - 4) + (у1 - у2)(с2 - 4^
< С2' = с2 + ^(ж1 - х2)(с* - с*) + кз(3 - 2і)(4 - с2), (3)
I с3' = с3 + к2(У* - у2)(с1 - с1) + к3(3 - 2і)(с3 - с1).
Здесь і = 1, 2 — значение индекса, задающее выбранную и соседнюю клетки с координатами (ж1, у1) и (ж2, у2) соответственно; с* и с*' — значения ячеек в моменты времени і и і +1 соответственно: с1 = Р, с2 = Ух, с3 = Уу; к1 ,к2,к3 > 0 — коэффициенты пропорциональности, которые отражают пространственно-временные параметры клеточно-автоматной модели, а также удельную плотность и кинематическую вязкость. В частности, смысл коэффициента к1 состоит в определении интенсивности изменения давления при различных значениях скорости потока вещества во взаимодействующих клетках. Коэффициент к2 определяет интенсивность изменения компонентов векторов скоростей пропорционально соответствующему градиенту давления в двух соседних
клетках. Коэффициент к3 определяет интенсивность «диффузии скорости», то есть является клеточно-автоматным отображением кинематической вязкости. Если выразить все величины в системе Си, то размерность коэффициента к1 — Па-с/м, коэффициента к2 —обратная величина м/(Па-с), а коэффициент к3 является безразмерным параметром. Тут следует подчеркнуть, что разность координат двух соседних клеток (ж1 — ж2) или (у1 — у2) в системе (3) необходимо воспринимать как одно из значений — 1, 0, +1, а не как некоторое расстояние.
Система уравнений (3) описывает изменение значений соответствующих ячеек во время одного элементарного взаимодействия клеточных автоматов и является своеобразной аппроксимацией численного решения уравнения Навье — Стокса (1), (2). Систему (3) можно представить в виде с* = с* + Ас*. Совершенно очевидно, что для обеспечения адекватности процесса моделирования прирост или уменьшение некоторого значения ячейки Ас* должны быть как можно меньшей величиной, регулируемой коэффициентами к1, к2, к3, при этом размер клеточно-автоматного поля должен быть как можно больше. С другой стороны, это неизбежно приведёт к слишком долгому процессу моделирования. Здесь по аналогии с различными многочисленными схемами решения возникает проблема поиска компромисса между точностью решения и временем, необходимым для его получения. Количественный анализ точности описываемого подхода не рассматривается и является предметом дальнейших научных исследований. Основной целью данной работы является построение качественной клеточно-автоматной модели турбулентной динамики.
Объясним вид системы уравнений (3), рассмотрев для простоты и наглядности одномерный случай (у1 = у2) без учета вязкости (к3 = 0). При этом система (3) примет вид
Г 4' = 4 + к1(ж1 — ж2)(с2— 4^ (4)
\ с*2' = с*2 + к2(ж1 — ж2)(с1 — с1).
Первое уравнение системы отображает модификацию давления Р для двух взаимодействующих соседних клеток. Если рассматривать две клетки как некий микрофрагмент системы, то в первом приближении можно считать, что изменение давления АР пропорционально разности векторов скоростей. Для наглядности обратимся к рис. 4, а. Пусть выбранная клетка (г = 1) с координатой ж1 находится слева, а соседняя клетка (г = 2) с координатой ж2 = ж1 + 1 — справа. Пусть векторы скоростей направлены вправо, причем V > V2 > 0. Тогда как для избранной, так и для соседней клеток прирост давления составляет АР = Ас1 = к1(ж1 — ж2)(с2 — С,) = к(—1)("2 — V!).
Так как V > V2 > 0, то АР > 0, что и видно на рис. 4,а. Если поменять места-
ми выбранную клетку (г = 1) и соседнюю (г = 2), картина не изменится, поскольку АР = к1(+1)("2 — У1), но в этом случае V > "1, значит, АР > 0. Можно проанализировать и другие случаи (рис. 4, б-г), причём все они должны отражать физику моделируемого явления, то есть давление должно увеличиваться, когда приток вещества превышает отток, и наоборот.
В частности, на рис. 4,б показан случай, когда V > V > 0, следовательно, при этом АР = к1 (—1)("2 — "1) < 0. На рис. 4,в,г показаны случаи, когда векторы скоростей направлены влево, то есть имеют отрицательные значения. При этом рис. 4,в соответствует случаю V < "2 < 0, и здесь результат АР < 0 получается аналогично изображенному на рис. 4,б. Рис. 4,г соответствует случаю "2 < V < 0, и здесь результат АР > 0 получается аналогично изображенному на рис. 4,а.
Второе уравнение системы (4) отражает изменение вектора скорости микрофрагмента вещества под действием градиента давления. Как можно видеть, оно полностью
V V
* >
|-ДР<0
( ] Ч Г2>0 1 ]
Гг >0 я У1 >0 О Г2>0
А А
;§ || }др<0 <—?
1 г. ч^
\\ < 0 У2 < 0 Ух< 0 о V |ЬГ
е
-АР > О
Рис. 4. Демонстрация правил клеточно-автоматных взаимодействий в случае изменения давления
аналогично первому уравнению системы. Для объяснения обратимся к рис. 5. Здесь в упрощенном виде показана модификация векторов скоростей двух взаимодействующих клеток под действием градиента давления. Если сравнить рис. 4, а,б с рис. 5 и отметить, что давление всегда является положительной величиной, то можно увидеть прямую аналогию между приращением или уменьшением скорости за счет разности давлений (рис. 5) и соответствующим приращением или уменьшением давления, обусловленным разностью скоростей в двух соседних клетках (рис. 4, а, б).
Рис. 5. Демонстрация правил клеточно-автоматных взаимодействий в случае модификации векторов скорости
Здесь следует заметить, что изменение вектора скорости за единицу времени (ускорение) пропорционально не только градиенту давления, но и массе вещества в элементе объёма, которая, в свою очередь, является функцией плотности, зависящей от давления сжимаемой среды. Но мы, в первом приближении, пренебрегаем этой зависимостью.
Рассмотрим слагаемые в системе (3), содержащие коэффициент пропорциональности к3. Они отражают вязкость (трение соседних слоев вещества) и по сути пред-
ставляют собой диффузионную составляющую компонентов векторов скорости. Эти слагаемые уравновешивают значения векторов скорости в соседних клетках.
Проиллюстрируем этот процесс с помощью рис. 6. Пусть в момент времени £ взаимодействуют две клетки: выбранная (г =1) с вектором скорости У и соседняя (г = 2) с вектором скорости У2, причем VI > У2. Проанализируем динамику изменения скорости. Для выбранной клетки А У = кз(3 — 2)(У2 — У) = к3(У2 — VI) < 0, а для соседней АУ2 = к3(3 — 4)(У2 — У!) = к3(—1)(V2 — V) > 0, что видно из рис. 6.
Если поменять местами выбранную клетку с соседней, результат остаётся без изменений — компоненты векторов скоростей в следующий момент времени £ + 1 стремятся к усреднённому значению.
Рис. 6. Иллюстрация влияния составляющей системы (3), которая учитывает вязкость вещества — векторы скорости двух соседних взаимодействующих клеток стремятся к усреднённому значению
До сих пор мы рассматривали «микровзаимодействия» на уровне отдельных клеток. Если запустить асинхронный итерационный цикл взаимодействия на множестве клеточных автоматов (клеточно-автоматном поле), который описывается системой (3), то при отклонении от состояния равновесия можно наблюдать колебательную динамику, как это схематически показано на рис. 7 для одномерного случая.
Рис. 7. Демонстрация примера клеточно-автоматной колебательной динамики
В то же время такое поведение системы клеточных автоматов никак не проявляет турбулентную динамику. Это связано с тем, что модель не предусматривает действительного перемещения вещества вместе с моментом импульса, а описывает лишь некую «упругую» среду. Для устранения этого недостатка дополним описанную выше систему итерационных функций (3) функцией перемещения, которая представляет собой реализацию вероятностного механизма дублирования соседней клетки. Для объяснения сути дополнения обратимся к иллюстрации перемещения отдельной частицы на
клеточно-автоматном поле. Пусть имеем элементарное поле размером 5 х 5 клеток, как это показано на рис. 8. Предположим наличие некоторой дискретной частицы, которая в начальный момент времени находится в нижнем левом углу поля. Заставим частицу передвигаться по полю вдоль вектора V, который приведен на рис. 8,а.
об в г д
Рис. 8. Демонстрация движения отдельной частицы на клеточно-автоматном поле: а — вектор скорости; б—д — возможные варианты траектории движения частицы
Для этого организуем вероятностный механизм передвижения частицы. Используя метод Монте-Карло, будем определять вероятность перехода частицы в соседнюю клетку. При этом свяжем вероятность перехода с длиной соответствующей компоненты вектора скорости (VX и Уу). Таким образом, если длина компоненты вектора скорости Уу в 2 раза больше длины VX, то в среднем частица будет переходить в верхние соседние клетки в 2 раза чаще, чем в клетки, расположенные справа от частицы. Следовательно, глобальная динамика поведения частицы — движение вдоль вектора скорости V, как показано на рис. 8,б-д. Кроме того, привязав максимальную величину векторов скорости к единичной вероятности осуществления перехода, можно организовать движение частиц с различными скоростями.
Подобным же образом введём описанный подход в нашей модели; существенным отличием здесь является не перемещение частиц в смысле обмена местами в клетках поля, а дублирование содержания соседних ячеек при осуществлении перемещения, как показано на рис. 9. Введение дублирования вместо обмена содержанием соседних клеток при перемещении определяется, в нашем случае, моделированием сплошного непрерывного потока, а не потока отдельных дискретных частиц.
Рис. 9. Правила клеточно-автоматных взаимодействий в случае перемещения вещества в пределах сплошного потока (дублирование при перемещении)
Несмотря на нарушение законов сохранения на микроуровне (во время локальных микровзаимодействий клеточных автоматов), глобальная макродинамика такой системы вполне адекватна, что можно проиллюстрировать с помощью рис. 10, где показан пример асинхронного вероятностного перемещения фрагмента сплошного вещества в виде некоторой совокупности частиц.
Рис. 10. Демонстрация перемещения вещества
3. Результаты вычислительного эксперимента
Для апробации предложенной в данной работе модели реализован вычислительный эксперимент по наблюдению турбулентной динамики на примере потока вещества в трубе, которая имеет препятствие (рис. 11, 12).
Рис. 11. Векторное поле скоростей
На рис. 11 и 12 приведён пример результата моделирования явления турбулентности на поле клеточного автомата размера 100 х 200. В частности, приведены векторное поле скоростей (рис. 11) и распределение давления (рис. 12).
Значения коэффициентов системы (3): к\ = 100; к2 = 0,00001; к3 = 0,1.
Граничные условия:
1) левая граница: Р = 100; Ух = 0,01;
2) верхняя и нижняя границы: Уу = 0;
Рис. 12. Распределение давления (черный цвет — максимальное давление, белый — минимальное)
3) препятствие (с обеих сторон): V; = 0;
4) правая граница—«прозрачная стенка», то есть тут Р(хп,у) = (хп — 1,у),
VX(xп,y) = VX(xп — 1,у), Vy(хп,у) = Vy(хп — 1,у), где хп — координата клетки, находящейся на правой границе.
Вероятность перемещения клетки (дублирования содержимого соседней клетки) определялась методом Монте-Карло. Если неравенство £ ^ '^^,у выполнялось, где ^ € [0,1] —некоторая случайная величина, а VX,y — X- или У-компонента вектора скорости, то перемещение клетки происходило.
Заключение
Предложена клеточно-автоматная аппроксимация уравнения Навье — Стокса, которое описывает динамику вещества в сплошных средах. Реализован вычислительный эксперимент, позволяющий исследовать турбулентную динамику на примере потока вещества в трубе с препятствием. Приведены векторное поле скорости и распределение давления для двумерного случая. Количественный анализ предложенной клеточно-автоматной модели — предмет дальнейших научных исследований.
ЛИТЕРАТУРА
1. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика. 2005. №10. С. 57-113.
2. Жихаревич В. В., Остапов С. Э. Моделирование процессов самоорганизации и эволюции систем методом непрерывных асинхронных клеточных автоматов // Компьютинг. 2009. Т. 8. №3. С. 61-71.
3. Темам Р. Уравнения Навье — Стокса. Теория и численный анализ. 2-е изд. М.: Мир, 1981. 408 с.
4. Ь^р://ги.и1к1ре^а.о^/и1к1/Турбулентность
5. Ь^р://ги.и1к1ре^а.о^/и1к1/Число_Рейнольдса
6. ЬМр://ги.и1к1ре^а.о^/и1к1/Существование_и_гладкость_решений_уравнений_ Навье_-_Стокса
