Моделирование возникновения неустойчивости в течении Тейлора — Куэтта методом конечных элементов с частицами PFEM-2

Язык труда и переводы:
УДК:
532.5
Дата публикации:
27 ноября 2022, 15:44
Категория:
Математическое моделирование физических процессов и технических систем
Авторы
Попов Андрей Юрьевич
МГТУ им. Н.Э. Баумана
Аннотация:
Рассмотрено моделирование неустойчивости течения вязкой несжимаемой жидкости на примере задачи о течении Тейлора — Куэтта о течении между двумя вращающимися цилиндрами. Показано отсутствие неустойчивости в двумерном расчете и ее возникновение в трехмерном. При исследовании использована разработанная на кафедре «Прикладная математика» МГТУ им. Н.Э. Баумана реализация метода конечных элементов с частицами PFEM-2. Возникновение неустойчивости продемонстрировано посредством возникающих вихревых структур по всей толщине жидкости между цилиндрами.
Ключевые слова:
метод конечных элементов с частицами, течение Тейлора — Куэтта, неустойчивое течение, свободное программное обеспечение
Основной текст труда

Введение

При моделировании течений несжимаемых жидкостей и газов большой практический интерес представляет переход к их неустойчивости, однако расчет в таких задачах характеризуется высокой вычислительной трудоемкостью по причине трехмерной природы соответствующих явлений, а также необходимостью использования очень мелких сеток (в сеточных методах) либо большого количества частиц (в методах частиц) для корректного разрешения мелкомасштабных эффектов. В настоящей работе применительно к одной из задач о возникновении неустойчивости рассматривается сравнительно новый гибридный эйлерово-лагранжев метод конечных элементов с частицами (Particle finite element method, 2nd generation, PFEM-2), который путем комбинирования расчета на сетке и перемещения частиц позволяет снизить вычислительную нагрузку, особенно в задачах с преобладанием переноса среды. Для расчетов используется программная реализация метода, разработанная на кафедре «Прикладная математика» МГТУ им. Н.Э. Баумана, допускающая возможность расчета в том числе и в параллельном режиме на многопроцессорной системе. В данной работе мы ограничиваемся модельной задачей о течении Тейлора — Куэтта между двумя вращающимися коаксиальными цилиндрами [1, с. 659] и демонстрируем возникновение неустойчивости при трехмерном расчете в отличие от двумерного.

Постановка задачи

 

Дж. Тейлор первым привел пример ламинарного течения, в котором при определенных условиях возможен переход к неустойчивости [2]. Мы будем рассматривать течение между двумя коаксиальными цилиндрами радиусов R_{1} и R_{2} ( R_{1}<R_{2} ), которые вращаются с постоянными угловыми скоростями \omega _{1} и \omega _{2} , соответственно (течение Тейлора — Куэтта). В силу симметрии поле скоростей будет иметь только окружную компоненту, зависящую только от радиальной координаты:

V_{\varphi }=C_{1}r+C_{2}{\dfrac {1}{r}}, где константы C_{1} и C_{2} определяются из граничных условий и равны соответственно

C_{1}={\dfrac {\omega _{1}R_{1}^{2}-\omega _{2}R_{2}^{2}}{R_{1}^{2}-R_{2}^{2}}},\qquad C_{2}={\dfrac {R_{1}^{2}R_{2}^{2}(\omega _{2}-\omega _{1})}{R_{1}^{2}-R_{2}^{2}}}.

Для случая вращения цилиндров в одну сторону ( \omega _{1} и \omega _{2} одного знака) Сайндж привел критерий устойчивости течения [3]

\omega _{1}R_{1}^{2}<\omega _{2}R_{2}^{2}.

Если же цилиндры вращаются в разные стороны (пусть для определенности \omega _{1}<0,\omega _{2}>0 ) течение в целом будет неустойчивым.

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

Используемый метод

В данной работе для моделирования возникновения неустойчивости в течении Тейлора — Куэтта используется метод конечных элементов с частицами PFEM-2 [4]. Производится решение системы из уравнений Навье — Стокса и уравнения несжимаемости

\left\{{\begin{aligned}\rho \left({\dfrac {\partial {\boldsymbol {V}}}{\partial t}}+({\boldsymbol {V}}\cdot \nabla ){\boldsymbol {V}}\right)=\nabla p+\nabla \cdot {\boldsymbol {\tau }}+{\boldsymbol {g}},\\\nabla \cdot {\boldsymbol {V}}=0.\end{aligned}}\right. Здесь {\boldsymbol {V}} — поле скоростей;  p — поле давления;  \rho — плотность; {\boldsymbol {g}} — векторное поле массовых сил; {\boldsymbol {\tau }} — девиатор тензора напряжений с компонентами ( \mu — динамическая вязкость);

\tau _{ij}=\mu \left({\dfrac {\partial V_{i}}{\partial x_{j}}}+{\dfrac {\partial V_{j}}{\partial x_{i}}}-{\dfrac {2}{3}}\delta _{ij}{\dfrac {\partial V_{k}}{\partial x_{k}}}\right), где i,j,k=1,\ldots ,d (предполагается суммирование по повторяющимся индексам);  d — размерность задачи.

Будучи эйлерово-лагранжевым методом, метод PFEM-2 предполагает разделение исходной задачи на две по физическим процессам:

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

Особенно эффективным PFEM-2 оказывается в задачах с преобладанием переноса среды (convection-dominated problems), поскольку за счет исключения моделирования переноса из задачи МКЭ на сетке позволяет моделировать остальные процессы на более грубых сетках и с более крупным шагом по времени без потери устойчивости. PFEM-2 также применим и в задачах, где большую роль играет вязкость жидкости, что мы и демонстрируем на примере течения Тейлора — Куэтта в настоящей работе.

Лагранжевы частицы являются важной частью метода PFEM-2. Перемещение частиц производится по известному полю скоростей (оно считается в этот момент «замороженным») и выполняется в несколько шагов, так чтобы число Куранта было не более 0,10...0,15 (для глобальной задачи допускается число Куранта \geqslant 1 ). Технически эта операция может выполняться с помощью явного метода Эйлера либо Рунге — Кутты, скорости при этом вычисляются путем интерполяции скоростей в узлах. Отдельно производится решение уравнений Навье — Стокса (без конвективного слагаемого) и уравнение несжимаемости на неподвижной эйлеровой сетке. Для их конечноэлементного решения могут применяться совместные (monolithic) подходы и подходы «с расщеплением» (segregated). В данной работе мы используем метод дробных шагов [5], который относится к последним и включает в себя три шага ( \Delta t — шаг задачи по времени):

  1. {\dfrac {{\boldsymbol {v}}^{n+1/2}-{\boldsymbol {v}}^{n}}{\Delta t}}=-{\dfrac {1}{\rho }}\nabla p^{n}+\nabla \cdot {\boldsymbol {\tau }}^{n+1/2}+{\boldsymbol {g}} — прогноз скорости;
  2. \Delta p^{n+1}=\Delta p^{n}+{\dfrac {\rho }{\Delta t}}\nabla \cdot {\boldsymbol {v}}^{n+1/2} — поиск поля давления;
  3. {\dfrac {{\boldsymbol {v}}^{n+1}-{\boldsymbol {v}}^{n+1/2}}{\Delta t}}=-{\dfrac {1}{\rho }}\left(\nabla p^{n+1}-\nabla p^{n}\right) — коррекция скорости.

В то время как поле прогноза скорости {\boldsymbol {v}}^{n+1/2} не обязано быть бездивергентным, результирующее поле скоростей {\boldsymbol {v}}^{n+1} является таковым и удовлетворяет уравнению несжимаемости.

Обратным эффектом от разделения задачи на «лагранжеву» (частицы) и «эйлерову» (сеточное решение) части является наличие двух наборов переменных (полей скоростей) — в частицах и в узлах сетки, — которые требуют согласования между собой. Оно обеспечивается за счет проекционных процедур в обоих направлениях. После перемещения частиц скорости в узлах сетки обнуляются и вычисляются заново путем проецирования скоростей частиц. На скорость j -го узла влияют скорости {\boldsymbol {u}}_{p} частиц в ячейках сетки, которым он принадлежит:
{\boldsymbol {V}}_{j}={\dfrac {\sum \limits _{p=1}^{P}{\boldsymbol {u}}_{p}\varphi _{j}({\boldsymbol {x}}_{p})}{\sum \limits _{p=1}^{P}\varphi _{j}({\boldsymbol {x}}_{p})}}. Здесь \varphi _{j}({\boldsymbol {x}}_{p}) — значение функции формы j -го узла, вычисленное в координатах p -й частицы.

После нахождения решения на сетке методом конечных элементов оно также используется для корректирования поля скоростей в частицах, однако без полного обнуления скоростей в них, а только путем коррекции на величину изменения поля скоростей в узлах i содержащей ее ячейки:
\delta {\boldsymbol {u}}_{p}^{n}({\boldsymbol {x}}_{p}^{k})=\sum \limits _{i=1}^{M}\varphi _{i}({\boldsymbol {x}}_{p}^{k})({\boldsymbol {V}}_{i}^{n+1}-{\boldsymbol {V}}_{i}^{n}).

Результаты численного моделирования

Для численного моделирования течения Тейлора — Куэтта использовалась программная реализация метода PFEM-2, разработанная на кафедре «Прикладная математика» МГТУ им. Н.Э. Баумана [6]. Эта программный код является частью фреймворка с открытым исходным кодом, который покрывает полный цикл решения задачи. Ядром является собственный (in-house) решатель,   основанный на конечноэлементной библиотеке deal.II на языке C++. В качестве препроцессора для построения сетки и выделения граничных участков может использоваться пакет SALOME, данные из которого затем импортируются решателем. В процессе расчета решатель формирует данные в формате VTK, которые могут визуализироваться, например, широко распространенным постпроцессором Paraview.

В рамках данной работы рассматривались двумерная и трехмерная задачи. В обоих случаях использовались следующие параметры: кинематическая вязкость \nu =0,001 , плотность \rho =1 , радиусы цилиндров R_{1}=0,5,R_{2}=1 . Моделирование производилось путем счета на установление. В двумерной задаче использовалась сетка из 2560 ячеек и шаг расчета по времени \Delta t=0,01 . Рассматривалось вращение цилиндров в одну сторону с угловыми скоростями \omega _{1}=8,\omega _{2}=1, что при указанных радиусах цилиндров должно приводить к неустойчивости согласно критерию выше. Произведенный расчет показал, что неустойчивость в этом случае не возникает. Стационарное поле скоростей приведено на рис. 1 (а — сеточное решение, б — поле в частицах).

Рис. 1. Поле скоростей в течении Тейлора — Куэтта в двумерном случае

В трехмерном случае использовалась расчетная область высотой h=5 , сечение которой аналогично двумерному случаю. Построенная сетка при этом состояла из 381440 ячеек, был принят шаг расчета \Delta t=0,0025 . Расчет производился в параллельном режиме (технология MPI) на кластере кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана с использованием в общей сложности 108 вычислительных ядер. Отметим, что в этом случае уже возникает движение жидкости в направлении координаты z с образованием структур, показанных на рис. 2 (изображено сечение x=0 в различные моменты времени). Возникновение неустойчивости при этом выражается в образовании вихревых ячеек, которые видны на линиях токах поля скоростей на рис. 3. Эти вихри имеют противоположные вращения и занимают всю область между цилиндрами.

Рис. 2. Возникновение неустойчивости в течении Тейлора — Куэтта в трехмерном случае (компонента V_z)
Рис. 3. Вихри в сечении x = 0 при вращении цилиндров в одну сторону

Заключение

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

Литература
  1. Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. 2. Москва, Физматгиз, 1963, 728 с.
  2. Taylor G. Stability of a viscous liquid contained between two rotating cylinders. Philosophical Transactions of the Royal Society of London, 1923, vol. 223, no. 605–615, pp. 289–343. DOI: https://doi.org/10.1098/rsta.1923.0008
  3. Synge J. On the stability of a viscous liquid between rotating coaxial cylinders. Proceedings of the Royal Society, London, 1938, vol. 167, no. 929, pp. 250–256. DOI: https://doi.org/10.1098/rspa.1938.0130
  4. Idelsohn S.R., Nigro N., Gimenez J., Rossi R., Marti J. A fast and accurate method to solve the incompressible navier-stokes equations. Engineering Computations, 2013, vol. 30 (2), pp. 197–222. DOI: https://doi.org/10.1108/02644401311304854
  5. Ferziger J.H., Peri'c M., Street R.L. Computational Methods for Fluid Dynamics. Springer, 2002, 596 p. DOI: https://doi.org/10.1007/978-3-319-99693-6
  6. Popov A., Marchevsky I. MPI-based PFEM-2 method solver for convection-dominated CFD problems. Communications in Computer and Information Science, Springer, 2022, vol. 1618, pp. 261–275. DOI: https://doi.org/10.1007/978-3-031-11623-0\_18
Ваш браузер устарел и не обеспечивает полноценную и безопасную работу с сайтом.
Установите актуальную версию вашего браузера или одну из современных альтернатив.