CC BY
77
ПРИКЛАДНАЯ МАТЕМАТИКА И МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
УДК 532.527
И. К. Марчевский, Г. А. Щеглов
МОДЕЛЬ СИММЕТРИЧНОГО ВОРТОНА-ОТРЕЗКА ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ПРОСТРАНСТВЕННЫХ ТЕЧЕНИЙ ИДЕАЛЬНОЙ НЕСЖИМАЕМОЙ СРЕДЫ
Рассмотрен вихревой метод моделирования пространственных течений идеальной несжимаемой среды. В качестве вихревого элемента используется вортон-отрезок, который описывается координатами центра, направляющим вектором и интенсивностью. Построена математическая модель движения такого вортона и рассмотрена тестовая задача о "чехарде" вихревых колец. Показано, что предложенная модель вортона-отрезка позволяет получать более точные результаты по сравнению с другими моделями вор-тонов.
Для эффективного решения инженерных задач, связанных с исследованием динамики конструкций в дозвуковом потоке жидкости или газа, в ряде случаев используются бессеточные вихревые методы расчета течений с использованием упрощенных моделей среды [1, 2].
Вихревые методы [3, 4] основаны на рассмотрении эволюции поля завихренности, что позволяет сосредоточить вычислительные ресурсы в ограниченной области вихревого следа. Это дает возможность получать приемлемые для инженерных расчетов значения нестационарных нагрузок при меньших по сравнению с сеточными методами затратах машинного времени. Для исследования пространственных течений вихревыми методами используются различные трехмерные вихревые элементы — вортоны: вихревая точка [5], вихревой сгусток (vortex blob) [4], вихревой отрезок [6], сферический вортон [7].
Вортон представляет собой элемент вихревой трубки и характеризуется положением в пространстве и вектором интенсивности. Каждый тип вортонов задает некоторое распределение завихренности в пространстве, а скорость, индуцируемая вортоном, вычисляется по закону Био-Савара с учетом этого распределения. Область течения, занятая завихренностью, моделируется большим числом вортонов [6].
При моделировании эволюции завихренности в пространственном течении векторы интенсивностей вортонов изменяются в соответствии с уравнением движения среды. Аппроксимация непрерывного распределения завихренности конечным числом вихревых элементов
может приводить к накоплению значительной погрешности, которая проявляется в получении неверных результатов либо вычислительной неустойчивости [4]. Можно выделить три основных способа уменьшения погрешности вычислений: увеличение числа вортонов, моделирующих течение; уменьшение шага расчета по времени; совершенствование модели вортона. В настоящей работе предложена новая модель вортона-отрезка.
Вортон-отрезок представляет собой прямолинейный участок бесконечно тонкой вихревой нити. Вортоны-отрезки характеризуются интенсивностью и положением в пространстве двух характерных точек, которые движутся по траекториям жидких частиц. Обычно в качестве этих точек выбирают концы отрезка [6]. Однако определение скоростей вблизи краев вортона-отрезка приводит к большим погрешностям при вычислении скорости изменения вектора интенсивности вортона.
Этого недостатка лишен сферический вортон [4], который задает сферически симметричное распределение завихренности и характеризуется центром, вектором интенсивности и параметрами распределения завихренности. Однако сферический вортон не является подобным цилиндрическому элементу вихревой трубки, поэтому для моделирования вихревых структур с достаточной точностью требуется большое число вортонов и специальные приемы моделирования, например метод core spreading [7]. В то же время вортон-отрезок представляется более продуктивным элементом для моделирования вихревых трубок, так как он может естественным образом удлиняться и укорачиваться.
Цель настоящей работы — создание модели вортона-отрезка, которая предполагает вычисление скорости на максимально возможном удалении от краев, т.е. в центре отрезка.
Модель симметричного вортона-отрезка. Рассматривается течение идеальной несжимаемой среды, которое описывается уравнением неразрывности и уравнением Эйлера, записанным в форме Гельмголь-ца для завихренности Q = rot V [8]:
При этом скорость V в точке с радиус-вектором г вычисляется по известному распределению завихренности П(£):
divV = 0; — = rot(V х V).
(1)
где интеграл берется по всей области течения D.
Рис. 1. Вортон-отрезок
Завихренность движется в пространстве по траекториям жидких частиц:
dir -t
dt = V ^ (2)
Вихревой отрезок (вортон-отрезок) [9] характеризуется радиус-векторами концов r\ и r2 и интенсивностью Г (рис. 1).
Скорость, индуцируемая вортоном-отрезком в точке r, равна
Г
V (f) = 4ncä, (3)
где
c =
а = (r — fi) x (r2 — fi) ; 1 / (f —fi ) • (f2 — fi ) (f — f2 ) • (f2 —fi )\ (4)
V |r - ri | |r - r21
Обычно правая часть уравнения (2) вычисляется для концов каждого вортона непосредственно по формулам (3)-(4).
Рассмотрим иной подход к моделированию эволюции вортонов. Каждый вортон будем характеризовать положением его центра г0, вектором Аг (см. рис. 1) и интенсивностью Г. Такой вортон является симметричным относительно точки г0. Предполагается, что в процессе эволюции вихревого элемента точка г0 движется по траектории жидких частиц и всегда остается центром вортона. Ниже выведены уравнения, описывающие изменение Аг в процессе движения.
Уравнения движения симметричного вортона. Запишем уравнение (2) для концов вортона:
d (го + Аг) г
--Л-= V (го + Аг) ;
, dt (5)
d (г0 - Аг) V (5) ---= V (го - Аг).
Скорость V (г) может быть разложена в ряд Тейлора в окрестности
точки г0 с точностью до членов второго порядка:
V (ГО ± Ar) = V (ГО) ± grad V (Го) • Дг + O (|ДГ]2) .
Вычитая второе уравнение (5) из первого и пренебрегая членами высшего порядка малости, получаем
d (ДГ) = grad V (r0) • Ar. (6)
Движение радиус-вектора центра вортона описывается уравнением
§ = * (Го) • (7)
Таким образом, движение симметричного вортона описывается системой обыкновенных дифференциальных уравнений (ОДУ) 6-го порядка (6)-(7). Интенсивность вортона Г остается неизменной [9].
Моделирование пространственного течения. Если распределение завихренности Q в пространственном течении аппроксимируется множеством из N вортонов-отрезков, интеграл в формуле (1) заменяется суммой, слагаемые которой учитывают влияние отдельных ворто-нов и вычисляются по формулам (3)-(4). Для симметричных вортонов требуется решать систему ОДУ, полученную из уравнений (6)-(7):
Е V (ад
(8)
dt
j=i
d(Arj) X N
dt . \j=i
(rciH • A ri, i = 1,N,
где г0г = (хо,уо,го)Т — центр ¿-го вортона; V}- (гог) — скорость, индуцированная ]-м вортоном в точке гог; элементы матрицы В;(гог) = = gradV} (гог) вычисляются путем дифференцирования выражения (3):
( ÔVjdxo ÖVX/дуо ÖVx/ÖZo \
В3 (Гог) =
(9)
ЗУу/дхо дУу /дуо дУу /дго
\ дУ*/дхо дУг/дуо дУ*/дго )
Ниже приведены необходимые соотношения для вычисления компо нент скорости У; (гог) и их производных, входящих в матрицу (9):
Зо = Го; - Гог; ¿1 = Зо + Дг}; в2 = во - Дг}; в1 = |в1|; ¿2 = 1 ¿21; Дг2 = |Дг}|2;
Vi = рт; V2 = р^г; а = 2 (s0 х Arj); h = |a|2 ; |s1| |s2|
c = 2((a - f • ^1 ; K = (si + s2) Ar2 - (si - s2) (A-rj • SO) dh
= -8 ((s0)n Ar? - (Arj)n • (Ar' • SO)) ;
д (rOi)
dc 2 / dh .
scTOx: = ^ (S1 - S2) (Ar'- äcTO-ynK1 +
2
+ 777^2 ((ArJ ^ Sl) S2 - (ArJ' ^ S2) Ъ) . h si S?
Тогда
Г
^ (r0i) = ca;
è (Г_0г)т: = ôCrôf = ГП (dCrOö: • (a)m ^ C ^^ (Af' )k
d (ro)n (гог
" ' (10) Здесь каждый из индексов m, n, k принимает значения x, y, z; emnfe — символ Леви-Чивиты.
Расстояние от точки rO до оси j-го вортона (прямой, проходящей через вортон-отрезок) вычисляется по формуле
1 | a |
' 2 | Arj Г
' 3 I
Чтобы исключить неограниченный рост скоростей V■ (гог) и их производных при приближении к осям вортонов, вводится радиус дискретности е по аналогии с вихрем Рэнкина [8]. Считается, что внутри цилиндра радиуса е вблизи каждого вортона, т.е. при < е, индуцированные им скорости и их производные убывают по линейному закону до нуля на оси вортона.
В этом случае выражения (10) модифицируются:
V; (гог) = а д (г0г); В* (гог)тп = а Вз (г0г)тп , (11)
где
г0г = Гог + <Т (а1 - 1) ; а = ^Т; Г = ДГз ^ - 3).
Таким образом, выражения (10) и (11) полностью определяют правые части системы (8), которую условно запишем в виде
Г= / (Г) . (12)
Начальными условиями для системы (12) являются положения вортонов в момент времени £ = 0, т.е. Г (0) = Го.
:
В качестве метода интегрирования задачи Коши выбран метод Рунге-Кутты второго порядка:
Qi = Qi + f (qi) -у;
(fi+i = qi + / (qi*) At.
Рис. 2. Вихревые кольца
Тестирование симметричного ворто-
на. Вопросы аппроксимации поля завихренности вортонами и индуцированного ими поля скоростей рассмотрены в работах [6, 10]. В данной работе проведено численное исследование эволюции завихренности. Для тестирования предлагаемой расчетной схемы выбран известный качественный эффект "чехарды" двух вихревых колец. Известно [11], что в идеальной жидкости такое движение будет периодическим. Нарушение периодичности в результате численного моделирования позволит судить об адекватности применяемой модели вихревого элемента.
Рассматриваются два тонких вихревых кольца радиусом Ri = R2 = = R = 1 (рис. 2). В начальный момент времени центры колец лежат на одной оси, а расстояние между плоскостями колец h = 1,2. Циркуляции колец одинаковы и равны Г = 1.
Дискретизация вихревых колец была выполнена, как показано на рис. 2: каждое вихревое кольцо моделировалось при помощи N = 60 симметричных вортонов-отрезков, расставленных равномерно по окружности кольца. Радиус дискретности вортонов был выбран равным £ = 0,1, шаг интегрирования по времени At = 0,01.
В результате моделирования получены зависимости радиусов колец (рис. 3, а) и расстояния между кольцами (рис. 3, б) от времени. На рис. 4 приведены зависимости радиусов обоих колец от расстояния между ними.
Движение колец, как следует из графиков, является периодическим, что согласуется с теоретическими результатами исследования "чехарды" вихревых колец в идеальной жидкости [11]. При движении кольца не выходят из плоскости, сохраняют свою форму и через равные промежутки времени восстанавливают взаимное расположение в пространстве. На рис. 5 показаны отдельные фазы половины периода этого процесса.
Для сравнения проведено аналогичное исследование с использованием сферических вортонов методом core spreading [7]. В этом случае каждое вихревое кольцо моделировалось при помощи N = 360 сферических вортонов. Для обеспечения устойчивости вычислительного
1,4 i—
0 4-,-,-,-,-,-
0 20 40 60 80 100 t
а
Рис. 3. Зависимость от времени радиусов колец (а) и расстояния между кольцами (б) при использовании симметричных вортонов-отрезков
1,2 -1 f\.
1,0 -Л Ö
U,о -Л
■0,6-Г\ Л .
U,4 -Л О
0,2 . А .
-1,2 -0,8 -0,4 0 0,4 0,8 1,2 Ь
Рис. 4. Зависимость радиусов колец от расстояния между ними при использовании симметричных вортонов-отрезков
¿ = 0 ¿ = 2,5 ...... ¿=5,0 .■■■■"■■■■.
¿ = 7,5 ¿=10,0
Рис. 5. Фазы процесса "чехарды" вихревых колец
1,4 1,2 1,0 0,8 0,6
0 10 20 30 40 50 60 t
а
h
1,5 |-
б
Рис. 6. Зависимость от времени радиусов колец (а) и расстояния между кольцами (б) при использовании сферических вортонов
RhR2 ■1,6-1—
1,4—
-------0,6--------
-------0,4--------
-------0,2--------
------ 01-------
-1,2 -0,8 -0,4 0 0,4 0,8 1,2 к
Рис. 7. Зависимость радиусов колец от расстояния между ними при использовании сферических вортонов:
-— ;-—
процесса шаг интегрирования по времени потребовалось уменьшить до АЬ = 0,001. Тем не менее в этом случае при более тщательной дискретизации колец и на порядок меньшем шаге по времени полученные результаты далеки от теоретических. Из графиков, изображенных на рис. 6-7, видно, что движение колец не является периодическим, хотя качественный эффект "чехарды" имеет место.
Причиной непериодического движения колец, по-видимому, является деградация завихренности — перераспределение завихренности в процессе ее эволюции: на каждом шаге расчета деформирующиеся вортоны заменяются эквивалентными сферами [7].
Отметим, что вихревые кольца, составленные из сферических вор-тонов, достаточно быстро выходят из плоскости и теряют форму, что отмечено в работе [4]. На рис. 8 показаны положения колец в моменты
¿=8,6
¿=16,2
¿ = 24,8
GO
00
00
¿ = 32,9
¿ = 43,0
¿=53,8
Рис. 8. Разрушение вихревых колец при их моделировании сферическими вор-тонами
времени, когда их взаиморасположение наиболее близко к первоначальному.
Выводы. Предложен вихревой элемент — симметричный вортон-отрезок и получены уравнения его движения. Сравнение со сферическим вортоном показало, что новый вихревой элемент позволяет более эффективно моделировать эволюцию завихренности в пространстве. Предполагается, что симметричный вортон-отрезок можно будет использовать для моделирования пространственного обтекания тел при сравнительно небольших вычислительных затратах.
Авторы благодарят Межведомственный суперкомпьютерный центр РАН за предоставленную возможность использования высокопроизводительных кластеров МВС-60001М, МВС-100.К.
1. Аэрогидроупругость конструкций / А.Г. Горшков, В.И. Морозов, А.Т. Пономарев и др. - М.: Физматлит, 2000. - 590 с.
2. Трехмерное отрывное обтекание тел произвольной формы / Под ред. С.М. Белоцерковского. - М.: ЦАГИ, 2000. - 265 с.
3. Корнев Н. В. Метод вихревых частиц и его приложение к задачам гидроаэродинамики корабля / Дисс____д-ра техн. наук. - СПб, 1998. - 184 с.
4. Winckelmans G. S., Leonard A. Contributions to vortex particle methods for the computation of three-dimensional incompressible unsteady flows // Journal of Comp. Physics, - 1993. - Vol. 109. - P. 247-273.
5. Новиков Е. А. Обобщенная динамика трехмерных вихревых особенностей (вортонов) // Журнал эксп. и теор. физики. — 1983. - T. 84, вып. 3. - С. 975-981.
6. Кирякин В. Ю., Сетуха А. В. О сходимости вихревого численного метода решения трехмерных уравнений Эйлера в лагранжевых координатах // Дифф. уравнения. - 2007. - Т. 43, № 9. - С. 1263-1276.
7. Ojim a A., Kamemot o K. Numerical simulation of unsteady flow around three dimensional bluff bodies by an advanced vortex method // JSME International Journal, Series B. - 2000. - Vol. 43, no. 2. - P. 127-135.
8. Лойцянский Л. Г. Механика жидкости и газа. - М.: Дрофа, 2003. - 840 с.
СПИСОК ЛИТЕРАТУРЫ
9. Нелинейная теория крыла и ее приложения / Т.О. Аубакиров, С.М. Бе-лоцерковский, А.И. Желанников и др. - Алматы: Гылым, 1997. - 448 с.
10. B e a l e J. T., M a j d a A. Vortex methods: Convergence in three dimensions // Mathematics of Computation. - 1982. - Vol. 39, no. 159. P. 1-27.
11. Гуржий А. А., Константинов М. Ю., Мелешко В. В. Взаимодействие коаксиальных вихревых колец в идеальной жидкости // Известия АН СССР. МЖГ. - 1988. - № 2. - С. 78-84.
Статья поступила в редакцию 28.01.2008
Илья Константинович Марчевский родился в 1983 г., окончил МГТУ им. Н.Э. Баумана в 2005 г., ассистент кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Автор 18 научных и учебно-методических работ в области исследования движения и устойчивости конструкций в потоке среды, вычислительной гидродинамики и элементарной математики.
I.K. Marchevskii (b. 1983) graduated from the Bauman Moscow State Technical University in 2007. Assistant of "Applied Mathematics" department of the Bauman Moscow State Technical University. Author of 20 publications in the field of study of motion and stability of constructions in flow of medium, computational hydrodynamics and elementary mathematics.
Георгий Александрович Щеглов родился в 1972 г., окончил МГТУ им. Н.Э. Баумана в 1996 г.. Канд. физ.-мат. наук, доцент кафедры "Аэрокосмические системы" МГТУ им. Н.Э. Баумана. Автор более 30 научных работ в области динамики и прочности конструкций аэрокосмических систем, вычислительных методов в гидроупругости.
G.A. Shcheglov (b. 1972) graduated from the Bauman Moscow State Technical University in 1996. Ph. D. (Phys.-Math.), assoc. professor of "Aerospace Systems" department of the Bauman Moscow State Technical University. Author of more than 30 publications in the field of dynamics and strength of spacecrafts, computational methods in hydroelasticity.
