CC BY
56
2008
НАУЧНЫЙ ВЕСТНИК МГТУ ГА серия Аэромеханика и прочность
№ 125
УДК 519.6:517.9
О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТРЕХМЕРНЫХ ВИХРЕВЫХ ТЕЧЕНИЙ ИДЕАЛЬНОЙ ЖИДКОСТИ В БЕЗГРАНИЧНОЙ ОБЛАСТИ ИЗОЛИРОВАННЫМИ ВИХРЕВЫМИ ЭЛЕМЕНТАМИ
Д.В. БОГОМОЛОВ, А.В. СЕТУХА
Сформулирована модификация вихревого численного метода для трехмерного уравнения переноса завихренности в идеальной жидкости, основанная на использовании изолированных вихревых отрезков. Проведены методические исследования по регуляризации интегрального выражения для скорости через завихренность по закону Био-Савара. Осуществлено численное решение задачи о чехарде вихревых колец конечной толщины.
В работах С.М. Белоцерковского и его последователей [1, 2] развиты численные схемы решения трехмерных задач обтекания тел, основанные на аппроксимации поверхностей тел и вихревых областей в жидкости системами вихревых отрезков и замкнутых вихревых рамок. Существенный упор делался на отражение свойства замкнутости вихревых линий в дискретной модели. Последнее требование, однако, вызывает определенные трудности при практических расчетах, вызванные сильным растяжением вихревых элементов.
Альтернативный подход основан на использовании изолированных вихревых элементов [3, 4]. В настоящей работе построена численная схема такого типа с использованием изолированных вихревых отрезков. При построении численной схемы предварительно осуществлена регуляризация интегрального представления для поля скоростей через завихренность путем умножения ядра интеграла на регуляризирующий множитель. Отличия предложенной схемы состоят в том, что вычисление сглаженного поля скоростей осуществляется по сравнительно простым формулам с применением только арифметических операций, а учет изменения завихренности в точках жидкости осуществляется за счет изменения длины и ориентации вихревых отрезков.
Исходные уравнения течения идеальной несжимаемой жидкости
Пусть M = (x, у, z) - точка, в которой расположена жидкая частица, причем M = M(X, t), где X = (Xb Х2, X3) - лагранжевы координаты частицы, t - время. Пусть w = w(M, t) - скорость жидкости, w = rot w - завихренность. Движение точки M и изменение вектора завихренности W в точке M описываются уравнениями:
dM d7 dw Эw dw
-----=w, —=wx—+wy—+wz—, (1)
dt dt Эх Эу 3z
где d/dt - производная по времени в индивидуальной жидкой частице, т.е. при X = const. Пусть X есть начальное положение частицы и пусть Wo есть начальное распределение завихренности, причем вектор w Ф 0 только в некоторой ограниченной области D. Тогда течение жидкости описывается уравнениями (1) с начальными условиями:
W = Wo, M = X при t = 0. (2)
Вектор w в точке M выражается через завихренность с помощью закона Био-Савара:
w=4 Ш 77"dX XdX ydXZ, (3)
4 D M
где w - завихренность в точке N с лагранжевыми координатами (Xx,Xу,Xz), r=NM . Интеграл в правой части выражения (3) является несобственным. Осуществим его регуляризацию
за счет искусственного сглаживания особенности. Введем параметр е > 0 - радиус области сглаживания вокруг особой точки. Будем решать задачу (1), (2), где
* =7-^1 “хуе(г)С<Д;, Уо(г) = (4я)-уг3г, Уе(г) = 6е(|г|)Уо(г), е > 0, 0е(|г|) = в1(|г|/е), 4- 0
в1(г) = С3[0,¥) - некоторая функция, удовлетворяющая условиям в1(0) = в1'(0) =в1"(0) = 0, в1"(г) < Сг3 при всех г е [0, ¥), С - некоторая константа, в1(г) = 1 при г > 1. Указанная функция в1(г) в расчетах бралась в виде в1(г) = (63г5 - 90г7 + 35г9)/8 при г < 1.
Дискретизация задачи
Осуществим разбиение области Б на ячейки Ц|, ] = 1,..., п так, что область Б есть объединение всех ячеек Ц| и пересечение двух разных ячеек есть множество нулевого объема.
Пусть ^ есть объем ячейки Б|. В каждой ячейке Ц| выберем точку М0 . Построим для каждой точки М0 закрепленный вектор (т.е. вектор с фиксированными началом и концом)
10= К0м0 длины 1 и число Г (циркуляцию отрезка) так, чтобы точка М0 являлась серединой этого отрезка и чтобы выполнялось равенство Гj 10= ^ ю0(М^. Построенные вихревые элементы (10, Г 0) аппроксимируют начальное распределение завихренности.
Далее, в каждый момент времени гк = кДг, к = 1, 2,..., где Дг > 0 - шаг дискретизации по времени, будем аппроксимировать поле завихренности вихревыми элементами (1 'к, Г к), где
1 к = ^ , а числа Г считаем не зависящими от времени.
Пусть (1 j,Гj) - некоторая система вихревых элементов и пусть Mj - центр отрезка ^. Определим поле скоростей, индуцируемое этими вихревыми элементами, по формуле:
П >
и(М)=X1 j хУе(г)Гj, где г = MjM. (4)
j=l
Перемещение дискретных вихревых отрезков по времени будем аппроксимировать схемой второго порядка точности типа метода Эйлера с уточнением:
N+1 = N + и^кж Кк+1 = кк + и(кк)Дг, 1=1,...,п,
Кк+1 = Кк + и к Дг, Кк+1 = кк+и к Дг, 1=1,...,п, где и к = 0,5(и(Кк)+и(Кк+1)), векторы и и и определяются формулой (4) по вихревым элементам (1 к,Гj) и (1|к,Гj), соответственно, где 1^к = К^Мк .
Аппроксимация поля скоростей на начальном шаге
Рассмотрим поле скоростей, индуцируемое вихревым кольцом в форме тора (рис. 1). Пусть центр тора совпадает с началом координат, Я - радиус окружности, вокруг которой образован тор, и г0 - радиус круга в сечении тора. Положение точки М = (х, у, ъ) внутри тора определяется координатами г, ф, у: х = (Я + гсо8у)соБф, у = (Я + гсоБу^іпф, ъ = геіпу.
Пусть поле завихренности Юо распределено так, что Юо(М) есть вектор, сонаправленный
с вектором кхОМ , где к - орт оси Оъ, и |ю0| = (г0 - г)/г0.
Построим систему точек Мук с тороидальными координатами гі = (і - 0,5)Аг, і = 1,...,Кг, У = (І - 1)Ау, І = 1,...,Ку, фк = кАф, к = 1,...,Кф, где Аг = г0/Кф, Аф = 2л/Кф, Ау = 2я/КУ. Каждая
точка лежит в ячейке Бук, определяемой неравенствами |г - г;| < Аг/2, |у - ур| < Ау/2, |ф - фк| < Аф/2. Объем ячейки Dijk приближенно определяется по формуле ёук = (Я + ^совур^АгАуАф.
г | М у
\ ^
к
Рис. 1.
Пусть Тр есть орт вектора Юо в точке Мр. Положим 1 = 0,8ЯАф и построим совокупность векторов 10к и Г ijk так, что Гр/10к = dijk•Юо(Mijk).
Для тора с радиусом г0 = 0,25 были получены распределения вектора скорости и по формуле (4) (при сквозной перенумерации вихревых элементов) при различных значениях параметров е, Кг (при условии N = 2,5Кг, N = 5Кг). На рис. 2 показаны зависимости компоненты вектора скорости и2 от координаты х на оси Ох. Анализ этих результатов показывает, что при малом значении параметра е возникают значительные погрешности при вычислении скорости и нефизичные выбросы. Параметр е (радиус вихря) должен быть существенно больше, чем расстояния между вихрями. Уменьшение параметра е можно производить только с одновременным измельчением разбиения.
Рис. 2.
Для оценки точности аппроксимации поля скоростей были рассчитаны значения циркуляции О вектора скорости по окружностям различных радиусов г, лежащим в плоскости 0x2 и имеющим центр в точке (1, 0, 0). На рис. 3 приведены зависимости в(г), полученные в расчетах для наборов параметров е = 110-2, 7,510-2, 510-2, и К- = 4, 8, 16, соответственно, а также аналитическая зависимость, полученная применением теоремы Стокса. Видно, что уже при _2
е = 7,510 и К- = 8 достигается хорошее совпадение численной и аналитической зависимостей.
Численное моделирование движения двух вихревых колец
Для проверки работоспособности полной схемы решения нестационарной задачи было осуществлено численное моделирование движения системы из двух вихревых колец. Одно кольцо имело конфигурацию и распределение завихренности, описанные выше с г0 = 0,1, и для него использовалось разбиение на вихревые элементы, описанное выше. Второе кольцо получалось смещением первого по оси 02 на расстояние И = 0,5. Расчеты проводились при значениях вычислительных параметров е = 0,05, Кг = 8 при условии Ку = 2,5К К = 5К при различных значениях шага дискретизации по времени А1 = 0,1; 0,2; 0,5; 1,0.
На рис. 4 показаны траектории движения центра нижнего вихревого кольца в сечении плоскостью 0x2. При этом центр вихревого кольца определялся как центр масс вихревых элементов, которые на первом шаге располагались в сечении у = ух. Маркерами отмечены положения центров вихревых колец в контрольные моменты времени.
— Д1=1
— Дї=0.5
— М=0.2
— М=0А ♦ 1=100 ■ 1=150 ж 1=200
Рис. 4.
На рис. 5 показаны поля скоростей в плоскости 0x7 в начальный момент времени и для значения безразмерного времени 1 = 380 (т.е. в момент, когда нижнее кольцо сделало 1.5 оборота).
Выводы
В расчетах удалось получить эффект чехарды вихревых колец (прохождения одного кольца сквозь другое). При этом удается получить устойчивые решения на временном интер-
вале, за который кольца делают два оборота, причем для различных значений вычислительных параметров получены решения, хорошо согласующиеся между собой. Дальше начинает наблюдаться эффект диффузии колец и их взаимного разрушения.
t = 0
t = 380
Рис. 5.
ЛИТЕРАТУРА
1. Белоцерковский С.М. Ништ М.И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. - М.: Наука. 1978.
2. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. - М.: Наука. 1995.
3. Басин М.А., Корнев Н.В. Аппроксимация вихревого поля в безграничной среде // ЖТФ. 1994. № 64. С. 179 - 185.
4. Hou T.J., Lowengrib J. Convergence of the Point Vortex Method for the 3-D Euler equations // Comm Pure Appl. Math., V.43, 1990, pp. 965 - 981.
ON NUMERICAL SIMULATION 3-D VORTEX INVISCID FLOW IN UNBOUNDED SPACE USING THE ISOLATED VORTEX ELEMENTS
Bogomolov D.V., Setukha A.V.
Vortex numerical method with isolated vortex elements was formulated for the three-dimensional vorticity transport equation for an ideal incompressible fluid. The regularization of numerical scheme for Biot-Savart law was researched. The numerical solution for the problem of the two vortex rings motion was obtained.
Сведения об авторах
Богомолов Дмитрий Валерьевич, 1977 г.р., окончил ВВИА им. проф. Н.Е. Жуковского (2005), адъюнкт кафедры высшей математики ВВИА им. проф. Н.Е. Жуковского, автор 5 научных работ, область научных интересов - вычислительная аэродинамика, вихревые методы.
Сетуха Алексей Викторович, 1966 г.р., окончил МГУ (1988), доктор физико-математических наук, доцент, заведующий кафедрой высшей математики ВВИА им. проф. Н.Е. Жуковского, автор более 60 научных работ, область научных интересов - вычислительная аэродинамика, вихревые методы, интегральные уравнения, краевые задачи.
