При моделировании течений несжимаемых жидкостей и газов большой практический интерес представляет переход к их неустойчивости, однако расчет в таких задачах характеризуется высокой вычислительной трудоемкостью по причине трехмерной природы соответствующих явлений, а также необходимостью использования очень мелких сеток (в сеточных методах) либо большого количества частиц (в методах частиц) для корректного разрешения мелкомасштабных эффектов. В настоящей работе применительно к одной из задач о возникновении неустойчивости рассматривается сравнительно новый гибридный эйлерово-лагранжев метод конечных элементов с частицами (Particle finite element method, 2nd generation, PFEM-2), который путем комбинирования расчета на сетке и перемещения частиц позволяет снизить вычислительную нагрузку, особенно в задачах с преобладанием переноса среды. Для расчетов используется программная реализация метода, разработанная на кафедре «Прикладная математика» МГТУ им. Н.Э. Баумана, допускающая возможность расчета в том числе и в параллельном режиме на многопроцессорной системе. В данной работе мы ограничиваемся модельной задачей о течении Тейлора — Куэтта между двумя вращающимися коаксиальными цилиндрами [1, с. 659] и демонстрируем возникновение неустойчивости при трехмерном расчете в отличие от двумерного.
Дж. Тейлор первым привел пример ламинарного течения, в котором при определенных условиях возможен переход к неустойчивости [2]. Мы будем рассматривать течение между двумя коаксиальными цилиндрами радиусов и (), которые вращаются с постоянными угловыми скоростями и , соответственно (течение Тейлора — Куэтта). В силу симметрии поле скоростей будет иметь только окружную компоненту, зависящую только от радиальной координаты:
где константы и определяются из граничных условий и равны соответственно
Для случая вращения цилиндров в одну сторону ( и одного знака) Сайндж привел критерий устойчивости течения [3]
Если же цилиндры вращаются в разные стороны (пусть для определенности ) течение в целом будет неустойчивым.
Стоит отметить, что при моделировании течения Тейлора — Куэтта в двумерном случае неустойчивости не возникает в силу несжимаемости жидкости, а истинную картину позволяет передать только трехмерное моделирование. Дополнительно отметим, что хотя вязкость жидкости численным образом никак не определяет решение (поле скоростей зависит только от радиусов цилиндров и скорости и направления их вращения), но определяющей причиной, вызывающей течение Тейлора — Куэтта, является именно влияние вязкости.
В данной работе для моделирования возникновения неустойчивости в течении Тейлора — Куэтта используется метод конечных элементов с частицами PFEM-2 [4]. Производится решение системы из уравнений Навье — Стокса и уравнения несжимаемости
Здесь — поле скоростей; — поле давления; — плотность; — векторное поле массовых сил; — девиатор тензора напряжений с компонентами ( — динамическая вязкость);
где (предполагается суммирование по повторяющимся индексам); — размерность задачи.
Будучи эйлерово-лагранжевым методом, метод PFEM-2 предполагает разделение исходной задачи на две по физическим процессам:
Особенно эффективным PFEM-2 оказывается в задачах с преобладанием переноса среды (convection-dominated problems), поскольку за счет исключения моделирования переноса из задачи МКЭ на сетке позволяет моделировать остальные процессы на более грубых сетках и с более крупным шагом по времени без потери устойчивости. PFEM-2 также применим и в задачах, где большую роль играет вязкость жидкости, что мы и демонстрируем на примере течения Тейлора — Куэтта в настоящей работе.
Лагранжевы частицы являются важной частью метода PFEM-2. Перемещение частиц производится по известному полю скоростей (оно считается в этот момент «замороженным») и выполняется в несколько шагов, так чтобы число Куранта было не более 0,10...0,15 (для глобальной задачи допускается число Куранта ). Технически эта операция может выполняться с помощью явного метода Эйлера либо Рунге — Кутты, скорости при этом вычисляются путем интерполяции скоростей в узлах. Отдельно производится решение уравнений Навье — Стокса (без конвективного слагаемого) и уравнение несжимаемости на неподвижной эйлеровой сетке. Для их конечноэлементного решения могут применяться совместные (monolithic) подходы и подходы «с расщеплением» (segregated). В данной работе мы используем метод дробных шагов [5], который относится к последним и включает в себя три шага ( — шаг задачи по времени):
В то время как поле прогноза скорости не обязано быть бездивергентным, результирующее поле скоростей является таковым и удовлетворяет уравнению несжимаемости.
Обратным эффектом от разделения задачи на «лагранжеву» (частицы) и «эйлерову» (сеточное решение) части является наличие двух наборов переменных (полей скоростей) — в частицах и в узлах сетки, — которые требуют согласования между собой. Оно обеспечивается за счет проекционных процедур в обоих направлениях. После перемещения частиц скорости в узлах сетки обнуляются и вычисляются заново путем проецирования скоростей частиц. На скорость -го узла влияют скорости частиц в ячейках сетки, которым он принадлежит:
Здесь — значение функции формы -го узла, вычисленное в координатах -й частицы.
После нахождения решения на сетке методом конечных элементов оно также используется для корректирования поля скоростей в частицах, однако без полного обнуления скоростей в них, а только путем коррекции на величину изменения поля скоростей в узлах содержащей ее ячейки:
Для численного моделирования течения Тейлора — Куэтта использовалась программная реализация метода PFEM-2, разработанная на кафедре «Прикладная математика» МГТУ им. Н.Э. Баумана [6]. Эта программный код является частью фреймворка с открытым исходным кодом, который покрывает полный цикл решения задачи. Ядром является собственный (in-house) решатель, основанный на конечноэлементной библиотеке deal.II на языке C++. В качестве препроцессора для построения сетки и выделения граничных участков может использоваться пакет SALOME, данные из которого затем импортируются решателем. В процессе расчета решатель формирует данные в формате VTK, которые могут визуализироваться, например, широко распространенным постпроцессором Paraview.
В рамках данной работы рассматривались двумерная и трехмерная задачи. В обоих случаях использовались следующие параметры: кинематическая вязкость , плотность , радиусы цилиндров . Моделирование производилось путем счета на установление. В двумерной задаче использовалась сетка из 2560 ячеек и шаг расчета по времени . Рассматривалось вращение цилиндров в одну сторону с угловыми скоростями что при указанных радиусах цилиндров должно приводить к неустойчивости согласно критерию выше. Произведенный расчет показал, что неустойчивость в этом случае не возникает. Стационарное поле скоростей приведено на рис. 1 (а — сеточное решение, б — поле в частицах).
В трехмерном случае использовалась расчетная область высотой , сечение которой аналогично двумерному случаю. Построенная сетка при этом состояла из 381440 ячеек, был принят шаг расчета . Расчет производился в параллельном режиме (технология MPI) на кластере кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана с использованием в общей сложности 108 вычислительных ядер. Отметим, что в этом случае уже возникает движение жидкости в направлении координаты с образованием структур, показанных на рис. 2 (изображено сечение в различные моменты времени). Возникновение неустойчивости при этом выражается в образовании вихревых ячеек, которые видны на линиях токах поля скоростей на рис. 3. Эти вихри имеют противоположные вращения и занимают всю область между цилиндрами.
В настоящей работе на примере задачи для течения Тейлора — Куэтта показано, что при численном моделировании возникновение неустойчивости происходит в трехмерном случае, в отличие от двумерного. Расчет проведен с использованием ранее разработанной программной реализации метода конечных элементов с частицами PFEM-2. Дополнительный интерес представляет рассмотрение различных случаев соотношения угловых скоростей внутреннего и внешнего цилиндров, а также различных сочетаний направления их вращения.