CC BY
43
УДК 532.5.013
МОДЕЛИРОВАНИЕ ЭВОЛЮЦИИ ПЕРЕПЛЕТЕННЫХ ВИХРЕВЫХ НИТЕЙ МЕТОДОМ ВИХРЕВЫХ ЭЛЕМЕНТОВ
С.А. ДЕРГАЧЕВ, Г.А. ЩЕГЛОВ
Статья посвящена описанию новой модификации алгоритма расчета движения вихревых нитей. Представлены результаты методических исследований, показывающие, что данный алгоритм позволяет корректно моделировать растяжение, сжатие, изгиб и перезамыкание вихревых нитей, а также известные качественные явления, возникающие при взаимодействии вихревых колец. Особое внимание уделено повторению в расчете недавно исследованной экспериментально эволюции переплетенных вихрей. Показано, что с помощью разработанного алгоритма могут быть устойчиво смоделированы сложные эволюции вихревых петель. Главное достоинство нового алгоритма -устойчивость и высокая скорость вычислений.
Ключевые слова: переплетенные вихри, вихревые нити, метод вихревых элементов, вычислительная гидродинамика.
Введение
Применение вихревых бессеточных лагранжевых методов для математического моделирования нестационарных процессов, протекающих при взаимодействии элементов конструкций с потоком жидкости или газа в условиях, когда справедливо допущение о несжимаемости среды, часто является более эффективным с точки зрения затрат вычислительных ресурсов, чем использование сеточных методов [1; 2]. Наиболее развиты вихревые методы расчета плоскопараллельных и осесимметричных течений, в которых вихревые нити имеют постоянную форму прямых линий или круговых колец. Для описания эволюции вихревого следа за обтекаемым телом в этом случае достаточно рассматривать движение одной характерной точки-маркера для каждой из нитей [3; 4].
Расчет пространственных течений вихревыми методами требует учета изменения формы вихревых нитей: их искривления, растяжения, перезамыкания. Анализ литературы показывает, что имеется значительное количество подходов к моделированию указанных процессов как в классической, так и в сверхтекучей жидкости [5-8]. Параметризация формы пространственной вихревой нити, как правило, вводится путем задания на ней значительного количества маркеров [9-11]. Движение маркеров определяется уравнениями движения жидких частиц. Перезамыкание нитей моделируется при помощи изменения принадлежности маркеров, которое проводится по эмпирическим алгоритмам [10]. Для обеспечения вычислительной устойчивости вводятся дополнительные приемы, такие как сглаживание ядра, добавление маркеров на вихревой нити при ее удлинении, сглаживание нити, например путем удаления петель и пр.
При решении вихревыми методами задач обтекания тел требуется моделировать образующиеся за ними вихревые следы. Эти следы могут содержать большое количество протяженных вихревых нитей, которые перезамыкаются и деформируются в потоке сложным образом, образуя замкнутые петли и узлы. Выбор эмпирических моделей и разработка эффективных алгоритмов моделирования эволюции вихревых нитей при пространственном обтекании тел остаются актуальными задачами, поскольку при большом числе вихревых структур в следе вычислительные затраты на моделирование их эволюции оказываются значительными, а накапливающаяся вычислительная погрешность может приводить к ограничению конечного времени моделирования нестационарного процесса.
Целью данной работы является тестирование новой модификации алгоритма расчета движения вихревых структур в следе за обтекаемым телом. Основная задача методических иссле-
дований состоит в выборе параметров, позволяющих моделировать известные качественные явления, возникающие при взаимодействии изолированных вихревых петель. Особое внимание уделено повторению в расчете недавно исследованной экспериментально эволюции переплетенных вихрей [12].
Постановка задачи и метод решения
Рассматривается безграничный объем несжимаемой среды, с которым связана неподвижная система координат 0ХУ2. Принято допущение, что среда имеет малую вязкость, влияние которой проявляется только в процессах перезамыкания вихревых нитей. При этом движение среды описывается уравнением неразрывности и уравнением Гельмгольца:
divV = 0; — = , (1)
dt
где V(r, t) - нестационарное трехмерное поле скоростей; г - радиус-вектор точки в неподвижной мировой системе координат; О = V х V - завихренность. Задано граничное условие затухания возмущений на бесконечности
1imV(r,t)= 0.
В начальный момент времени в среде находится К замкнутых вихревых нитей, движение которых исследуется в течение промежутка времени 0 < X < Т.
Для численного моделирования используется метод вихревых элементов. Параметрами к-й вихревой нити являются координаты принадлежащих ей N точек-маркеров в неподвижной системе координат г и интенсивность Г, одинаковая для всех нитей. Движение вихревой нити описывается перемещением маркеров по траекториям жидких частиц:
-к-=Ук1; -(0)=-0; /=иN; к=1,..к. (2)
ш
При вычислении скорости маркеров считается, что отрезок вихревой нити между маркерами Л?к; = ?к;+1 - гк1 является прямолинейным. Вектор скорости маркера определяется суммарным влиянием всех отрезков, вычисленным по закону Био-Савара [ 1]
к N
- к N - / \
V* =1),
1=1 ]=1
- / ч I а, Г в, • Лг 5? • ЛГ ^ - - - -- -- -- .- - -
где V(-)=—гА- -^р^; а = 81 X Лгц; 81 = Г - Гц; 81 = Г - Гу+1; Л^ = ^ - ^.
V ^ Р2| )
I _а_
4 л а • а
Для исключения неограниченного роста скоростей (- ) и их производных при приближении к оси вихревого отрезка вводится сглаживающая функция [13]
В(Г) = Г3(1 + Г2)-3/2, где г = г/ в, а в > 0 - заданный радиус вихря.
На шаге интегрирования системы (2) положение маркеров и их принадлежность нитям определяется в соответствии с описанным ниже алгоритмом.
Алгоритм расчета движения вихревых нитей
Расчет проводится до достижения конечного времени счета Т. Алгоритм расчета одного шага по времени состоит из следующих операций:
Шаг 1. Определяются векторы скоростей маркеров V;.
Шаг 2. Определяются новые положения маркеров путем интегрирования системы (2) методом Эйлера 1 -го порядка точности с шагом А1.
Шаг 3. Вычисляется угол излома к-й нити в точке, где расположен маркер, как угол между
Aîki ЧЫ
). Если данный угол оказывается меньше заданно-
соседними отрезками = arccos(-
Ark! АГкЫ
го предельного угла ф, то, аналогично [10], проводится сглаживание нити путем коррекции положения маркера rki = (гкМ + rki+J/2. После коррекции угол излома петли в точке, где расположен маркер, становится равным yki = п. Данная операция продолжается, пока во всех маркерах не выполнится условие > ф.
Шаг 4. Производится перераспределение маркеров вдоль k-й нити для выравнивания длины отрезков на следующем шаге. Перед началом расчета задается номинальная длина отрезка h. Для уточнения формы нити между маркерами rki, rki+1 на текущем шаге производится построение кривых aki. Каждая кривая выражается полиномом третьего порядка, параметры которого определяются по координатам четырех маркеров aki(rkil, rki, rki+1, rki+2]. Длины кривых суммируются, и по величине h находится количество маркеров на новом шаге N (t + At). Полученная кривая а = У oki разбивается на равные по длине отрезки. i=1
Шаги 3-4 повторяются для всех контуров k = 1,...K.
Шаг 5. Осуществляется поиск областей перезамыкания нитей. Для каждого маркера rki осуществляется поиск ближайшего маркера Tj (l = 1,..K;j = 1,..N¡ ), для которого выполняются
условия:
rU " rki
DklDii - - - - - -< ^ ; arccos(-^ ) > ф ; Di. = ri j+1 - ri; Dkt = ik1+1 - ^ s
DkiDiJ
(3)
ротивном слу--ij _ -ki = m1n,
где д > 0 - заданная величина [14].
Если условия (3) не выполнены, то перезамыкания нитей не происходит. В пр
чае из удовлетворяющих условию (3) пар маркеров находится такая, для которой
и производится перезамыкание нити между этими маркерами. При этом либо объединяются две петли, либо одна петля разделяется на две.
Шаг 6. Определяется новый контур нити. Для этого вводится допустимая зона перезамыкания из В маркеров. В этой области определяются индексы четырех граничных маркеров Ь е [0, В] ( т = 1, 2, 3, 4 ) из условия максимального сокращения длины нити
j+b2
I
m=j-bi
Ar,
i m
1+b,
+ I
n=1—b.
Ar
kn
_ Ц
rki—b3 rij+b2 +1
+
rij—b rki+b4 +1
) = max
где X > 1 - заданный коэффициент.
Шаги 5-6 повторяются до прекращения перезамыканий. В случае если перезамыкание нитей имело место, осуществляется переход к шагу 3, иначе производится переход к шагу 1 для следующего момента времени.
При программной реализации алгоритм был распараллелен при помощи технологии MPI.
Тестирование алгоритма
Параметрами задачи, определяемыми из методических исследований, являются: s - радиус вихря; At - шаг интегрирования (2) по времени; h - номинальная длина отрезков вихревых нитей; Г - интенсивность вихревых нитей; ф - предельный угол излома нити; X - коэффициент сокращения длины нити при перезамыкании; д - максимальная дистанция для перезамыкания; B - область перезамыкания. Для определения параметров использовался ряд известных задач о движении вихревых колец.
На первом этапе тестирования воспроизводилось явление "чехарды" вихревых колец, ранее смоделированное при помощи других алгоритмов [11]. Данная задача позволила оценить накопление погрешности при работе алгоритма.
Рассчитывалось движение двух одинаковых круговых вихревых колец (K = 2) с единичной интенсивностью Г = 1 и радиусами R1 = R2 = 1, лежащих в параллельных плоскостях на начальном расстоянии b = 1,2.
Ряд проведенных расчетов позволил показать, что при дискретизации порядка 100 маркеров на кольцо удается длительное время моделировать данный неустойчивый процесс. В качестве примера на рис. 1 приведены фазы чехарды колец, которые параметризованы одинаковым числом маркеров N1 = N2 = 125. Расчет производился с параметрами: s = 0,1; At = 0,0005; h = 0,05; Ф = 115°; ц = 0,05; А, = 1,3; В = 10.
Рис. 1. Фазы процесса чехарды вихревых колец, полученные в расчете: первое кольцо изображено пунктирной линией, второе - сплошной
На рис. 2 приведены зависимость радиусов колец от времени и зависимость радиусов колец от их взаимного расположения. Видно, что до разрушения кольца совершают 6 полных циклов чехарды. При этом процесс с высокой точностью можно считать периодическим, что свидетельствует о высокой точности моделирования, с малыми потерями энергии и импульса из-за накопления погрешностей.
На втором этапе тестирования воспроизводилась эволюция овальных вихревых колец [15; 16]. Данная задача позволила оценить точность моделирования эффектов изгиба и перезамыкания вихревых нитей.
Ряд проведенных расчетов позволил показать, что все эффекты, наблюдаемые в эксперименте [15], воспроизводятся в расчете. В качестве примера на рис. 3 представлены результаты моделирования эволюции овального вихревого кольца большого удлинения единичной интенсивности Г = 1 с размерами: Ь = 11,0; Б = 2,0; Ь/Б = 5,5. Расчет производился с параметрами: К = 1; N1 = 485; в = 0,1; Лх = 0,0005; И = 0,05; ф = 115°; = 0,05; X = 1,3; В = 10.
-2 -1 0 1 2 0 50 100 150 200
а б
Рис. 2. Изменение параметров колец: а - зависимость радиусов от расстояния между кольцами;
б - зависимость радиусов колец от времени
Формы, принимаемые кольцом в характерные моменты времени в начале эволюции, приведены на рис. 3. Как видно из рисунка, овальное кольцо при первом колебании разделяется на два кольца, далее два кольца объединяются в одно, и происходит циклический изгиб кольца, характерный для овальных колец малого удлинения.
Рис. 3. Фазы эволюции овального вихревого кольца большого удлинения, полученные в расчете
На третьем этапе тестирования воспроизводилось явление столкновения двух круговых колец [14]. Данная задача позволила оценить точность моделирования эффектов растяжения и перезамыкания вихревых нитей.
Изучалось движение двух одинаковых круговых вихревых колец (К = 2) единичного радиуса = Я2 = 1 с единичными интенсивностями Г1 = 1, Г2 = -1, лежащих в начальный момент времени в параллельных плоскостях на расстоянии Ь = 4,0. В качестве примера на рис. 4 представлены три основные фазы процесса столкновения колец. Расчет производился с параметрами: N1 = N2 = 125; в = 0,1; М = 0,0005; И = 0,05; ф = 115°; = 0,4; X = 1,3; В = 10. Из рисунка видно, что кольца притягиваются друг к другу с увеличением радиуса. В определенный момент из двух колец возникает большое количество малых колец, что полностью соответствует картине, наблюдаемой в эксперименте [14].
0,0 с 4,0 с 8,0 с
полученные в расчете
Результаты моделирования эволюции переплетенных вихревых нитей
В эксперименте [12] при помощи специального генератора удалось получить несколько типов переплетенных вихревых колец и визуализировать их эволюцию [17]. При помощи разработанного алгоритма было рассчитано движение двух конфигураций колец из эксперимента.
Была исследована эволюция сцепленных вихревых колец единичной интенсивности Г = 1. Угол между плоскостями колец в начале моделирования составлял 30°. Расчет производился с параметрами: в = 0,05; М = 0,025; И = 0,05; ф = 130°; д = 0,085; X = 1,3; В = 10, К = 2, N1 = N = 127. На рис. 5 представлены три основные фазы эволюции колец в двух проекциях. Из рисунка видно, как сцепленные кольца перезамыкаются в два отдельных независимых кольца, что соответствует картине, наблюдаемой в эксперименте [12; 17].
Далее была исследована эволюция вихревого узла сложной формы, представленного на рис. 6 в двух проекциях. Параметры расчета: и = 100; К = 1; N = 898; в = 0,005; А1 = 0,01; И = 0,5; ф = 140°; д = 2,5; X = 1,3; В = 10. На рис. 5 показаны основные фазы эволюции нити: вихревой узел также трансформируется в два независимых кольца, что соответствует картине, наблюдаемой в эксперименте [12; 17].
Рис. 4. Фазы столкновения вихревых колец,
Рис. 5. Фазы эволюции сцепленных вихревых Рис. 6. Фазы эволюции вихревого узла, колец, полученные в расчете полученные в расчете
Заключение
Проведенные расчеты показали высокую эффективность разработанной модификации алгоритма. Тестовые расчеты эволюции вихревых колец показали корректность моделирования растяжения, сжатия, изгиба и перезамыкания вихревых нитей, а также позволили выбрать параметры расчетной схемы. Расчеты движения переплетенных вихревых нитей показали, что с помощью разработанного алгоритма могут быть устойчиво смоделированы сложные эволюции вихревых петель. Получено качественное сходство результатов расчетов с экспериментальными данными. Достигнутые результаты позволяют использовать разработанный алгоритм для моделирования вихревого следа за обтекаемым телом при помощи множества вихревых петель.
ЛИТЕРАТУРА
1. Трехмерное отрывное обтекание тел произвольной формы / под ред. С.М. Белоцерковского. М.: ЦАГИ, 2000. 265 с.
2. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. Cambridge: Cambridge University Press, 2000. 320 p.
3. Белоцерковский С.М., Гиневский А.С. Моделирование турбулентных струй и следов на основе метода дискретных вихрей. М.: Физматлит, 1995. 367 с.
4. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. М.: МГУ, 2006. 184 с.
5. Kida S., Takaoka M. Vortex Reconnection // Annual Review of Fluid Mechanics, 1994. Vol. 26. Pp.169 -177.
6. Sreenivasan K.R. Vortex reconnection classical and superfluid turbulence compared and contrasted // Proc. of Fluid-Gravity Correspondence Arnold Sommerfeld Centre, LMU, Munich 4 September 2009. [Электронный ресурс]. URL: https://www.theorie.physik.uni-muenchen.de/activities/workshops/archive_09/200909_2-7/talks/sreenivasan.pdf.
7. Hussain F., DuraisamyK. Mechanics of viscous vortex reconnection. Physics of Fluids, 2011. Vol. 23(2):021701.
8. Van Rees W.M., Hussain F., Koumoutsakos P. Vortex tube reconnection at Re=104. Physics of Fluids, 2012. Vol. 24(7): 075105.
9. Leonard A. Vortex Methods for Flow Simulation // Jour. Comput. Phys. 1980. № 37. Pp. 289-335.
10. Chorin A.J. Hairpin removal in vortex interactions II // Jour. of Comput. Phys. 1993. № 107. Pp. 1-9.
11. Богомолов Д.В., Марчевский И.К., Сетуха А.В., Щеглов Г.А. Численное моделирование движения пары вихревых колец в идеальной жидкости методами дискретных вихревых элементов // Инженерная физика. 2008. № 4. С. 8-14.
12. Kleckner D., Irvine W.T.M. Creation and dynamics of knotted vortices // Nature Physics. 2013. № 9. Pp. 253-258.
13. Winckelmans G.S., Leonard A. Contributions to Vortex Particle Methods for the Computation of Three-Dimensional Incompressible Unsteady Flows // Journal of Computational Physics. 1993. № 109. Pp. 247-273.
14. Weibmann S., Pinkall U., Filament-based smoke with vortex shedding and variational reconnection // ACM Transactions on Graphics. 2010. Vol. 29. № 4. Article 115.
15. Alkemade. A.J.Q. On Vortex Atoms and Vortons. PhD Thesis. Delft, (The Netherlands), 1994. 209 p.
16. Марчевский И.К., Щеглов Г.А. Численное моделирование эволюции вихревых колец в идеальной жидкости // Необратимые процессы в природе и технике: Труды Шестой Всероссийской конференции. М.: МГТУ им. Н.Э. Баумана, 2011. Ч.П. С. 98-101.
17. Kleckner D., Irvine W.T.M. Supplementary videos from the paper "Creation and dynamics of knotted vortices". [Электронный ресурс]. URL: http://www.youtube.com/watch?v=rcnw8NeJqjU.
VORTEX ELEMENT METHOD SIMULATION OF KNOTTED VORTEX FILAMENTS EVOLUTION
Dergachev S.A., Shcheglov G.A.
Research paper is devoted to the description of a new modification of the algorithm for calculating the vortex filaments motion. The results of methodological studies showing that this algorithm can correctly simulate the elongation, shortening, bending and reconnection of vortex filaments, as well as known qualitative phenomena arising from the vortex rings interaction are given. Particular attention is paid to the repetition of the calculation recently investigated experimentally the evolution of bound vortices. It is shown that with the developed algorithm can be steadily modeled complex vortex loops evolution. The main advantage of the new algorithm is stability and high-speed computing.
Keywords: knotted vortices, vortex filament, vortex element method, computational fluid dynamic.
REFERENCES
1. Trekhmernoe otrihvnoe obtekanie telproizvoljnoyj formih. pod red. S.M. Belocerkovskogo. M.: CAGI. 2000. 265 p. (In Russian).
2. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. Cambridge: Cambridge University Press. 2000. 320 p.
3. Belocerkovskiyj S.M., Ginevskiyj A.S. Modelirovanie turbulentnihkh struyj i sledov na osnove metoda diskret-nihkh vikhreyj. M.: Fizmatlit. 1995. 367 p. (In Russian).
4. Andronov P.R., Guvernyuk S.V., Dihnnikova G.Ya. Vikhrevihe metodih rascheta nestacionarnihkh gidrodina-micheskikh nagruzok. M.: MGU. 2006. 184 p. (In Russian).
5. Kida S., Takaoka M. Vortex Reconnection. Annual Review of Fluid Mechanics. 1994. Vol. 26. Pp. 169 -177.
6. Sreenivasan K.R. Vortex reconnection classical and superfluid turbulence compared and contrasted. Proc. of Fluid-Gravity Correspondence Arnold Sommerfeld Centre, LMU, Munich 4 September 2009. URL: https://www.theorie.physik.uni-muenchen.de/activities/workshops/archive_09/200909_2-7/talks/sreenivasan.pdf.
7. Hussain F., DuraisamyK. Mechanics of viscous vortex reconnection. Physics of Fluids. 2011. Vol. 23(2):021701.
8. Van Rees W.M., Hussain F., Koumoutsakos P. Vortex tube reconnection at Re=104. Physics of Fluids. 2012. Vol. 24(7): 075105.
9. Leonard A. Vortex Methods for Flow Simulation. Jour. Comput. Phys. 1980. № 37. Pp. 289-335.
10. Chorin A.J. Hairpin removal in vortex interactions II. Jour. of Comput. Phys. 1993. № 107. Pp. 1-9.
11. Bogomolov D.V., Marchevskiyj I.K., Setukha A.V., Theglov G.A. Chislennoe modelirovanie dvizhenija pary vihrevyh kolec v ideal'noj zhidkosti metodami diskretnyh vihrevyh jelementov. Inzhenernaia fizika. 2008. № 4. Pp.8-14. (In Russian).
12. Kleckner D., Irvine W.T.M. Creation and dynamics of knotted vortices. Nature Physics. 2013. № 9. Pp. 253-258.
13. Winckelmans G.S., Leonard A. Contributions to Vortex Particle Methods for the Computation of Three-Dimensional Incompressible Unsteady Flows. Journal of Computational Physics. 1993. № 109. Pp. 247-273.
14. Weibmann S., Pinkall U., Filament-based smoke with vortex shedding and variational reconnection. ACM Transactions on Graphics. 2010. Vol. 29. № 4. Article 115.
15. Alkemade. A.J.Q. On Vortex Atoms and Vortons. PhD Thesis. Delft. (The Netherlands). 1994. 209 p.
16. Marchevskiyj I.K., Theglov G.A. Chislennoe modelirovanie jevoljucii vihrevyh kolec v ideal'noj zhidkosti. Ne-obratimihe processih v prirode i tekhnike: Trudih Shestoyj Vserossiyjskoyj konferencii. M.: MGTU im. N.Eh. Baumana. 2011. Ch.II. Pp. 98-101. (In Russian).
17. Kleckner D., Irvine W.T.M. Supplementary videos from the paper "Creation and dynamics of knotted vortices". URL: http://www.youtube.com/watch?v=rcnw8NeJqjU.
Сведения об авторах
Дергачев Сергей Александрович, 1989 г.р., окончил МГТУ им. Н.Э. Баумана (2012), аспирант МГТУ им. Н.Э. Баумана, область научных интересов - математическое моделирование динамики элементов конструкций аэрокосмических систем.
Щеглов Георгий Александрович, 1972 г.р., окончил МГТУ им. Н.Э. Баумана (1996), доцент, доктор технических наук, профессор кафедры аэрокосмических систем МГТУ им. Н.Э. Баумана, автор более 90 научных работ, область научных интересов - вихревые методы вычислительной гидродинамики, математическое моделирование динамики элементов конструкций аэрокосмических систем.
