Научная статья на тему 'Т-СХЕМЫ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ГЕНЕРАЦИИ ЗАВИХРЕННОСТИ НА ГЛАДКИХ ПРОФИЛЯХ В ВИХРЕВЫХ МЕТОДАХ'

Т-СХЕМЫ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ГЕНЕРАЦИИ ЗАВИХРЕННОСТИ НА ГЛАДКИХ ПРОФИЛЯХ В ВИХРЕВЫХ МЕТОДАХ Текст научной статьи по специальности «Физика»

CC BY
122
12
Читать
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ВИХРЕВЫЕ МЕТОДЫ / ГЕНЕРАЦИЯ ЗАВИХРЕННОСТИ / ГРАНИЧНОЕ ИНТЕГРАЛЬНОЕ УРАВНЕНИЕ / МЕТОД ГАЛЕРКИНА / ВТОРОЙ ПОРЯДОК ТОЧНОСТИ

Аннотация научной статьи по физике, автор научной работы — Марчевский И. К., Сокол К. С., Измайлова Ю. А.

Предложены новые численные схемы решения граничного интегрального уравнения, решаемого в вихревых методах вычислительной гидродинамики при моделировании плоского обтекания гладких профилей. В их основе переход от традиционно рассматриваемых сингулярных интегральных уравнений первого рода с ядром Гильберта к уравнениям второго рода с ограниченным или абсолютно интегрируемым ядром относительно интенсивности вихревого слоя. Для их решения применяется метод Галеркина. Показано, что при аппроксимации границы профиля многоугольником можно построить схемы первого и второго порядка точности, рассматривая кусочно-постоянное или кусочно-линейное (разрывное или непрерывное) распределение решения по панелям. Приведены необходимые формулы для вычисления компонент матрицы и правой части системы линейных алгебраических уравнений --- дискретного аналога интегрального уравнения. Они пригодны для описания генерации завихренности при моделировании обтекания как одиночного профиля, так и системы профилей, в том числе подвижных и/или деформируемых. Построенные схемы могут быть использованы в рамках метода вязких вихревых доменов или иных модификаций вихревых методов, поскольку касаются лишь вопроса вычисления конвективных скоростей среды вблизи профиля и не связаны с методами моделирования диффузии завихренности

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по физике , автор научной работы — Марчевский И. К., Сокол К. С., Измайлова Ю. А.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
Предварительный просмотр
i Надоели баннеры? Вы всегда можете отключить рекламу.

T-SCHEMES FOR MATHEMATICAL MODELLING OF VORTICITY GENERATION ON SMOOTHS AIRFOILS IN VORTEX PARTICLE METHODS

New numerical schemes are proposed for solving the boundary integral equation that arises in CFD vortex particle methods of when simulating a plane flow around smooth airfoils. They are based on considering the 2-nd kind integral equation with respect to vortex sheet intensity with a bounded or absolutely integrable kernel instead of traditionally solved singular integral equations of the 1-st kind with Hilbert-type singularity. To solve it, the Galerkin approach is used. It is shown that when approximating the airfoil boundary with a polygon, it is possible to develop schemes of the 1-st and 2-nd order of accuracy, considering a piecewise-constant or piecewise-linear (discontinuous or continuous) distribution of the solution along the panels. The necessary formulae are presented for calculating the components of the matrix and the right-hand side of the system of linear algebraic equations, that is a discrete analogue of the integral equation. They are suitable for modelling of the vorticity generation when simulating the flow around either single airfoil or system of airfoils, including moving and/or deformable ones. The developed schemes can be used in the framework of the viscous vortex domains method as well as other modifications of vortex particle methods, since they only concern the convective velocities of the flow near the airfoil and are not related to methods for modeling viscous diffusion of vorticity

Текст научной работы на тему «Т-СХЕМЫ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ГЕНЕРАЦИИ ЗАВИХРЕННОСТИ НА ГЛАДКИХ ПРОФИЛЯХ В ВИХРЕВЫХ МЕТОДАХ»

УДК 532.5+519.64

DOI: 10.18698/1812-3368-2022-6-33-59

Г-СХЕМЫ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ГЕНЕРАЦИИ ЗАВИХРЕННОСТИ НА ГЛАДКИХ ПРОФИЛЯХ В ВИХРЕВЫХ МЕТОДАХ

И.К. Марчевский1, 2 К.С. Сокол1, 2 Ю.А. Измайлова1, 2

iliamarchevsky@bmstu.ru kuz-ksen-serg@yandex.ru yulia.izmailova@mail.ru

1 МГТУ им. Н.Э. Баумана, Москва, Российская Федерация

2 ИСП РАН, Москва, Российская Федерация

Аннотация

Предложены новые численные схемы решения граничного интегрального уравнения, решаемого в вихревых методах вычислительной гидродинамики при моделировании плоского обтекания гладких профилей. В их основе переход от традиционно рассматриваемых сингулярных интегральных уравнений первого рода с ядром Гильберта к уравнениям второго рода с ограниченным или абсолютно интегрируемым ядром относительно интенсивности вихревого слоя. Для их решения применяется метод Галеркина. Показано, что при аппроксимации границы профиля многоугольником можно построить схемы первого и второго порядка точности, рассматривая кусочно-постоянное или кусочно-линейное (разрывное или непрерывное) распределение решения по панелям. Приведены необходимые формулы для вычисления компонент матрицы и правой части системы линейных алгебраических уравнений — дискретного аналога интегрального уравнения. Они пригодны для описания генерации завихренности при моделировании обтекания как одиночного профиля, так и системы профилей, в том числе подвижных и/или деформируемых. Построенные схемы могут быть использованы в рамках метода вязких вихревых доменов или иных модификаций вихревых методов, поскольку касаются лишь вопроса вычисления конвективных скоростей среды вблизи профиля и не связаны с методами моделирования диффузии завихренности

Ключевые слова

Вихревые методы, генерация завихренности, граничное интегральное уравнение, метод Галеркина, второй порядок точности

Поступила 19.07.2022 Принята 05.08.2022 © Автор(ы), 2022

Работа выполнена при финансовой поддержке Минобрнауки России (проект № 075-15-2020-808)

Введение. Вихревые методы вычислительной гидродинамики имеют узкую область применения, которая ограничена рассмотрением несжимаемых изотермических течений, однако они широко применяются в приложениях [1-4], особенно при расчете гидродинамических нагрузок, действующих на обтекаемые поверхности в потоке среды. Они эффективны для моделирования нестационарных течений при интенсивном вихреобразовании, когда упрощенные модели нагрузок (модель потенциального обтекания, использование стационарных коэффициентов, теория присоединенных масс и т. п.) неприменимы или приводят к значительным расхождениям с данными экспериментов. Обзор современного состояния вихревых методов и истории их развития выполнен в [5-7].

К наиболее важной особенности вихревых методов следует отнести их бессеточный характер. Нематериальные лагранжевы частицы выступают в качестве носителей завихренности. Зная ее распределение, по закону Био — Савара может быть восстановлено поле скоростей, а с помощью аналога интеграла Коши — Лагранжа [8] — поле давления. При моделировании внешних течений вычислительные ресурсы «сосредоточиваются» в той части области течения, где завихренность отлична от нуля. Как правило, ее площадь сравнительно невелика, при этом обеспечено точное удовлетворение граничных условий затухания возмущений на бесконечности. Вихревые методы эффективны при моделировании обтекания подвижных и/или деформируемых тел при произвольных перемещениях и поворотах; для решения сопряженных задач взаимодействия тела и среды разработаны эффективные схемы связывания, пригодные и для расчета движения тел малой массы [9].

Вопросам моделирования эволюции завихренности в области течения посвящено множество исследований, к настоящему времени созданы эффективные алгоритмы расчета ее конвективного переноса и вязкой диффузии. Вопросы моделирования течений в пристеночных областях, где дополнительно имеет место генерация завихренности, исследованы значительно слабее, даже для плоских течений. В большей части публикаций приведены лишь общие рассуждения: для описания процесса генерации завихренности необходимо решать граничное интегральное уравнение (ГИУ) относительно потока завихренности, интенсивности вихревого слоя или плотности потенциала двойного слоя.

Сингулярным и гиперсингулярным ГИУ посвящена работа [10], в которой описаны и обоснованы схемы численного решения указанных задач, составляющие «ядро» метода дискретных вихрей (МДВ). Этот метод эффективен при расчете безотрывного обтекания гладких профилей или профи-

лей с острыми кромками и угловыми точками, в которых моделируется отрыв. Среда полагается идеальной, рассматривается граничное условие непротекания, а завихренность, генерируемая на профиле (кроме окрестностей точек отрыва), считается присоединенной, т. е. связанной с профилем и не являющейся «частью» области течения.

Скачок в развитии вихревых методов связан с появлением алгоритмов моделирования вязких течений, среди которых можно выделить метод вязких вихревых доменов (МВВД) [11], дающий эффективный способ моделирования диффузии завихренности во внутренней части области течения и вблизи обтекаемых профилей. На профилях, в общем случае подвижных и/или деформируемых, ставится условие прилипания, обеспечиваемое генерацией на всей их поверхности свободной завихренности, которая далее эволюционирует в области течения. Присоединенные вихревой слой и слой источников в МВВД также присутствуют, но их интенсивности определяются только скоростью движения границы профиля.

При описании генерации завихренности наиболее естественным представляется определение ее потока на границе профиля (скорости ее генерации); для него можно записать ГИУ, аналогичное по структуре ГИУ из МДВ. Однако алгоритмов, предполагающих его решение, не известно в силу сложности обеспечения балансовых соотношений. В связи с этим вместо потока завихренности рассматривается, по сути, интеграл от него — завихренность, генерируемая в течение малого промежутка времени (шага расчета). С учетом классических оценок можно пренебречь ее диффузией в область течения и рассматривать тонкий вихревой слой на границе профиля. Его интенсивность определяется решением ГИУ, также по форме аналогичного уравнению из МДВ. Тем не менее использование численных схем [10], созданных для МДВ, при моделировании вязких течений, сопровождаемых интенсивным вихреобразованием вблизи профиля, обеспечивает сравнительно невысокую точность, что снижает эффективность современных модификаций вихревых методов в целом.

Указанную проблему можно решить, переходя к эквивалентной формулировке задачи для ГИУ с абсолютно интегрируемым (или даже ограниченным) ядром. Идея такого подхода описана в [12], упоминания о нем имеются в [1, 10] и других работах, однако конкретных численных схем, его реализующих, до последнего времени не было известно.

В настоящей работе для моделирования генерации завихренности в вихревых методах предложены новые численные схемы первого и второго порядка точности, основанные на методе Галеркина и пригодные для расчета обтекания гладких профилей, для которых интенсивность вихре-

вого слоя является ограниченной и непрерывной вдоль профиля. Эти схемы могут быть применены к профилям с острыми кромками и угловыми точками, но в силу невозможности воспроизведения особенности решения их порядок точности будет ниже. Вопросы разработки специальных схем, позволяющих корректно разрешать особенности решения, нетривиальны, требуют отдельного рассмотрения и выходят за рамки настоящей работы.

Математическая модель. В соответствии с обобщенным разложением Гельмгольца [12-15] векторное поле V, заданное в области 5, может быть восстановлено по ротору й(4) = УхУ(4), т. е. завихренности, и дивергенции Я(4) = У^У(4), скорости Uк(4) на границе области 5 и скорости на бесконечности V», если область 5 не ограничена:

а(р)У(р) = УЮ+|(П (4) х Q (р-4) )+{ Я(4) Q (р-4) ^ +

5 5

+ ф( п(4) х ик (4) )х Q(р-4) ^ + ф( п(4) • ик (4)) Q (р-4) ^, (1)

к к

где д(4) = 4 / (2 л 1412) — градиент фундаментального решения двумерного уравнения Лапласа с обратным знаком; п(4) — орт нормали к границе, направленный внутрь области 5. Коэффициент а(р) равен единице в области 5, нулю вне этой области и внешнему углу границы в точке р, отнесенному к 2л. Таким образом, формула (1) формально верна во всей области М2; внутри области течения У(р) определяет «физическое» поле скоростей, вне области течения У(р) = 0. Под к в общем случае будем понимать объединение границ всех обтекаемых профилей и внешней границы области течения 5, если она присутствует. Поскольку У — это поле скоростей несжимаемой среды, то 0(4) = 0.

Процессы переноса и диффузии завихренности в области течения можно описать как результат ее переноса по полю скоростей (У + Ш), где Ш — диффузионная скорость [11, 16, 17]. Условие прилипания обеспечивается генерацией завихренности на границе к , для математического описания которой вместо работы с потоком завихренности [9] удобно оперировать с завихренностью, генерируемой за шаг расчета по времени. Пренебрегая диффузией [18], представим ее в виде тонкого вихревого слоя на границе профиля, который будем называть свободным, а его интенсивность у(4), 4 е к, подлежит определению.

Используя формулу Сохоцкого — Племеля, можно записать выражение для предельного значения поля скоростей на границе профиля (со стороны течения):

а(г)и К (Г) = Уоо+|(П(^) X д(г-4) + ф(у(4) X д(г-4))

$ К

- а(г) ( у (г) X п(г)) + ф ( уаП (4) X д(г -4) ) +

К

+ф qаtt(4)Q(г-4) dk>, г е К, (2)

К

где уаП(4) = п(4) х иК(4); ЦаП(4) = п(4) • иК(4) можно рассматривать как интенсивности присоединенного вихревого слоя и слоя источников. Интегралы по границе профиля в (2) — сингулярные, их следует понимать в смысле главного значения по Коши [10].

Выражение (2) — ГИУ относительно интенсивности свободного вихревого слоя. В плоских течениях вектор завихренности имеет единственную ненулевую компоненту, поэтому 0(4) = 0(4)к, у(4) = у(4)к, уаП(4) = _уап(4)к, где к — орт нормали к плоскости течения, таким образом, (2) — векторное ГИУ относительно скалярной неизвестной у(4):

ф (к х Q(r -4)) у(4)d^4 -а(г)(к х п(г))у(г) = {(г), г е К. (3)

К

Правая часть уравнения (3) зависит от известных величин:

{(г) = а(г)ик (г) - V» -1 (к х Q(r - 4)) 0(4^ -

$

- ф (к х Q(r -4)) ум(4)d^4 - ф qatt(4)Q(г-4) . (4)

К К

Для решения векторного ГИУ (3) с правой частью (4) достаточно решить любое из двух скалярных уравнений, получаемых проецированием его на орт нормали п(г) или орт касательной т(г) = к х п(г) [12]. В первом случае приходим к сингулярному ГИУ первого рода с ядром Гильберта:

ф Рп(г, 4)у(4)dk = /п(г), г е К, (5)

К

где

Рп (г, 4) = п(г) • ( к х Q(r - 4)) = - Т(г) ^(г -4), / (г) = п(г) • ((г).

2л | г -4 |2

Назовем ГИУ (5) ^-моделью; оно совпадает с ГИУ из МДВ [10], обычно получаемым из соотношений теории потенциала; схемы его решения будем называть ^-схемами.

В случае проецирования (3) на направление касательной получается ГИУ второго рода, которое назовем Т-моделью:

ф Рг(г, 4М4)^ - a(r)y(r) = f(r), r e K, (6)

K

где

Рг (r, 4) = x(r) • ( k X Q(r -4)) = n(r)'(r -4) , f (r) = x(r) • f (r).

2л | r - 4 г

Ядро уравнения (6) является ограниченным для профилей, заданных кривыми класса C2:

\nf ем к(г) lim | Рг (г, 4)|=——,

|r-41 ^ 0 4л

где к(г) — кривизна границы профиля в соответствующей точке; если r — утловая точка, то ядро Рг(г, 4) в ее окрестности не ограничено, но абсолютно интегрируемо.

Граничное интегральное уравнение (5) имеет бесконечное множество решений, для выделения единственного следует задать интеграл от решения

$y(4)dZ4 = Г. (7)

K

Если моделируется обтекание системы профилей, то граница области течения K состоит из нескольких компонент и условия вида (7) необходимо задавать для каждой компоненты. Для Т-модели (6) дополнительные условия необходимо задавать для профилей, которые обтекаются внешним образом, поскольку для профилей, ограничивающих область течения снаружи, такое условие получается как следствие уравнения (6). Отметим, что в [12] ошибочно указано, что существенным преимуществом Т-модели является единственность ее решения.

Расчетные схемы на основе метода Галеркина. Рассмотрим процедуру построения расчетных схем для численного решения уравнения (6), дополненного условиями выделения единственного решения (7). При наличии внешней границы области течения условие (7) для общности будем ставить и на ней, что корректно с учетом применяемой регуляризации.

Для численного решения ГИУ необходимо выбрать способ аппроксимации формы обтекаемых профилей; класс, в котором разыскивается приближенное решение; способ построения дискретного аналога ГИУ. Способы решения указанных задач, обеспечивающие различные порядки точности схем, перечислены в [1], в частности, отмечено, что для построения схемы второго порядка точности требуется явный учет кривизны границ профилей.

Ограничимся аппроксимацией границы профиля прямолинейными отрезками, далее называемыми панелями и образующими многоугольник, вписанный в профиль или описанный около него. Тем не менее такой подход позволяет построить схему второго порядка точности. Численное решение будем искать в классе кусочно-постоянных или кусочно-линейных функций, заданных на панелях. Общее число панелей на всех профилях обозначим за N, обозначение К теперь будет относиться к соответствующим многоугольникам.

Кусочно-постоянное численное решение можно представить в виде комбинации

N

у(г)= Еу0Ф0(г), г е К, (8)

I = 1

где уО — интенсивность вихревого слоя на 1-й панели ф0(г), I = 1,..., N, —

базисные функции — индикаторы панелей: ф0(г) = 1 при г е К/ и ф0(г ) = О на остальных панелях.

При кусочно-линейном представлении решения в дополнение к функциям {ф0(г)}^=1 введем кусочно-линейные базисные функции {ф/ (г)}^=1,

которые также финитны, т. е. ф1 (г) = О вне 1-й панели, и ортогональны функциям первого семейства:

Ф1(г)_(г ~ С) 'Ъ r G K

Фг (r)_-;-> r G Ki

где С/, , Ь/ — радиус-вектор центра, орт и длина 1-й панели. Приближенное решение ГИУ (6) будем искать в виде комбинации с неизвестными коэффициентами {у О}^! и {у1}^=х:

N , ч

у(г) = Х(у0ф0(г) + у1ф1(г)), Г е К, (9)

I = 1

где у0, у1 — среднее значение и вариация интенсивности вихревого слоя на 1-й панели. При таком представлении решение в общем случае разрывно на концах панелей.

Повышение порядка в представлении численного решения (введение квадратичных базисных функций) требует перехода к аппроксимации границы криволинейными панелями. Точное решение на профиле с «искусственно» образованными при дискретизации угловыми точками имеет в них слабые особенности. Квадратичные базисные функции позволяют их воспроизводить, и сходимости численного решения к точному на глад-

ком профиле в нормах ¿1 и ¿2 не наблюдается, тогда как при кусочно-постоянном или кусочно-линейном представлении решения сходимость имеет место при равномерной дискретизации профиля. В соответствии с результатами анализа погрешности аппроксимации интегрального оператора можно заключить, что уровень ошибки не возрастает при различии длин соседних панелей не более чем в 2 раза [19].

Для построения дискретного аналога ГИУ используем метод Галеркина, требуя ортогональности невязки уравнения (6) базисным функциям [20, 21].

Схема Т0 с кусочно-постоянным представлением решения. Результатом требования ортогональности функциям {ф0}^=1 невязки ГИУ (6) после подстановки в него решения в форме (8) является линейная система относительно неизвестных {у 0 }^=1:

1 м 11

— 1у0 I ^ | А(г,4)2у0= — I /т(г)Мг, (= 1,...,N. (10)

} = 1 К^ К] 2 К^

Здесь учтено, что а(г) = 1/2 всюду, кроме концов панелей, в которых образуется угловая точка; значения в них не влияют на результат интегрирования. Обе части каждого уравнения (10) поделены на длину соответствующей панели для улучшения обусловленности системы.

При моделировании обтекания одиночного профиля единственное решение ГИУ (6) обеспечивается условием (7), которое с учетом вида решения (8) принимает вид

N

„0г. -

= г. (11)

i = 1

Для решения системы (10), (11), включающей N +1 линейное уравнение относительно N неизвестных, вводится регуляризирующая переменная [10], в результате получается хорошо обусловленная система линейных алгебраических уравнений (СЛАУ) вида

[А00] + [D00] [IN}>

Шт 0

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Г1У0}Л

v R у

{b0}

v Г у

где [ А00] — квадратная матрица размером N х N; Р00] — диагональная матрица; {^} — матрица-столбец, составленная из единиц; {¿}т — матрица-строка из длин панелей; {у 0} = {у °)};[=1 — столбец неизвестных; Я —

регуляризирующая переменная; {Ь0} = {Ь°};[= — столбец правой части. Компоненты [А00], Р00] и |Ь0} вычисляются как

аО0= -1 | й1г | Рт(г, 4) ; С= "1/2; Ч К К]

Ь0 = -1 | / (Г) й\г, I, ] = 1,..., N.

Ч К/

Если рассматривается задача о моделировании обтекания системы из нескольких профилей, задаваемых контурами К «, 5 = 1,..., Np, и представляемых многоугольниками из панелей К(5) с длинами Ч(5), I = = 1,..., N(s), то построение схемы Т0 происходит аналогично. Вместо одного ГИУ (6) теперь необходимо решать систему интегральных уравнений вида

NP

X § Р (г, 4) у(4) ^-а(г) у(г) = / (Г), Г е К(5), 5 = 1,..., Np.

Я = 1 К (Я)

Дополнительные условия вида (7) ставятся для каждого контура по отдельности:

§ у(4) = г(5), 5 = 1,..., Np.

К 5)

Дискретные аналоги этих интегральных уравнений имеют вид, аналогичный (10), (11) (в обозначениях интенсивностей вихревого слоя на панелях верхний индекс (0) опустим):

1 Np 1

-15) X Ху? I ^ I Р(Г, 4)dk-1 у(5) =

Ч Я = 1 ] = 1 К(5) К(?) 2

1 ./

= I /т(г) ^, I = 1,..., N(5),

с ^

N(5)

X у;5)Ч(15) = Г(5), 5 = 1,..., Np.

I = 1

Регуляризированная СЛАУ относительно неизвестных величин у(я)

и дополнительных переменных Д(5) , соответствующих каждому профилю, в матричной форме принимает вид

[ А((1)] + [öf1)] [A«]

[ANip ]

[M(1)]

[ Ai(j)] [A(j)] + [D(j)]

ANp)]

[A,(Np)] Z

p

[M( j)]

[A(NP)] + [n(NP)] \7 ] ... [An/ ] + [DNp ] [ZnP ]

[M(NP)]

[O]

1у(1)| > ' т Л

1y( j)} {Bi}

lr(NP)} {BNp }

№ , 1B} J

Блоки [А-у^ матрицы составлены из компонент, учитывающих влияние завихренности с Р-й панели у-го профиля на а-ю панель г-го профиля, г, ] = 1,..., Ыр, а = 1,..., Ы(г), р = 1,..., Ы('

Г Л j)

[ Ai

(j h -

П,1

а

г«. j)

N(i ),1

а

(i, j) 1, N( j)

(i, j)

N(i), N(j)

а.

(i, j)_. 1

k® k(j)

a,ß

— J dlr J P (r, 4) d\,

блоки [у^ — нулевые матрицы при г ^ у и диагональные со значениями (-1/2) при г = у. Блоки столбца искомых величин {у(у)}, у = 1,..., Ыр, составлены из интенсивностей вихревого слоя на панелях у-го профиля:

{у(у)} = |у((У), ...,у(у(у) | . Блоки столбца правой части имеют вид

В } = (#,...,/2о )т, .Й0^ I Л«^, г = 1,..., Ы(г).

ка

Все компоненты матриц (г) х , г = 1,...,Ыр, равны нулю, кроме компонент г-го столбца, состоящего из единиц. В матрицах [М(хм(у),

у = 1,..., Ыр, в у-х строках находятся длины панелей соответствующих профилей, а остальные компоненты нулевые. Матрица [0]ыр х Ыр — нулевая; последний блок правой части {Вг} состоит из циркуляций Г( у).

Следовательно, процедура моделирования обтекания системы профилей принципиально не отличается от таковой для одиночного профиля, однако запись численной схемы оказывается весьма громоздкой. В связи с этим дальнейшее описание расчетных схем будем вести для случая одиночного профиля, имея в виду их возможное обобщение.

Cхема Т с кусочно-линейным представлением решения. При

кусочно-линейном представлении решения (9) дискретный аналог ГИУ (6) представляет собой СЛАУ вида

Г [А00] + [Р00] [А01] + [Р01] {^} У^П Г {Ь0}^ [А10] + [Р10] [А11] + Р11] {ON}

v {¿°}т № 0 у

{у 0} {у1}

V Ä У

{b1}

vry

где [АрЯ ] — квадратные блоки размером N х N; [РрЯ ] — диагональные блоки, р, ? = 0,1; {^}, {О^^} — N-мерные столбцы из единиц и нулей;

{Ч0}т, {¿1}т — N-мерные строки; {у0}, {у1} — ^мерные столбцы искомых коэффициентов в представлении решения; {Ь0}, {Ь1} — N-мерные векторы, образующие правую часть системы,

аРрЯ =1 ^р^Щ I Рт(г, 4)фЯ(4) ¿РЯ = I фр(г)фЯ(г)й/г,

Ki

K;

' Кг

Чр = I фр (г) й/г, Ьр = Ч I /т (г)фр (г) й/г, р, Я = 0,1, I, ] = 1,..., N.

К Ч К

Учитывая вид ядра Рт(г, 4) и постоянство щ, тг- вдоль панелей, для компонент [ АрЯ ] имеем

pq 1

j _ -г л

L'i

J фР(r) dlr J Q(r -4M(4) d\

Ki

Ki

■—л-^ pq

Li i •

Компоненты диагональных матриц [РрЯ ] и строк {Чр }т равны й00 = -1/2, 41 = -1/24, йО1 = й10=0, Ч0= Ч, Ч1=0, I = 1,..., N. Начала и концы панелей соответствуют направлению обхода контура К: при движении по нему область течения остается справа.

Компоненты блоков правой части представим в виде сумм Ьр =

= (Ьу )р + (Ьо )р + (Ьаи )р :

(Ьу )р = -1I (а(г)и к (г) - У® ) фр (г) й/г,

Ч К

(Ьо)р = --Ч-тг I фр(г)й/гКкх д(г-4))о(4) ,

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Ч К 5

1 ( > (Ъаи )Р = - - Ъ I фр (г) й\т $ (к х 0(г-4))уа« (4) ^ +$ 0(г-4) ^ (4) ^

ь' Кг V К к

Интегралы в слагаемых (Ъу )р могут быть вычислены точно: | (а(г )и к (г) - V» ) Ф0(г) й\г = Г 1 и к (с) - V»] Ь,

| (а(г )и к (г) - V» ) ф- (г) = -2 (и к (рг) - и к (с)), Кг 12

где с/, рг — радиус-векторы центра и начала г-й панели; ик (рг) - ик(сг) = = (-Ьгтг + Ьшгпг)/2; Ьг — скорость изменения длины г-й панели; Шг — ее угловая скорость.

При вычислении интеграла по области течения в слагаемых (Ъо )р будем полагать, что завихренность представлена множеством вихревых эле-

тЧО)

ментов — точечных вихревых нитей, имеющих циркуляции Г щ и расположенных в точках с радиус-векторами р^, щ = 1, ...,Ыо,

N О

О(р) = Хг1?)5(р-р^),

щ = 1

где 8(р) — дельта-функция Дирака. Тогда для интегралов, входящих в (Ъо )р, имеем

|фР (г) й\г |(к х д(г-4) )О(4) ^ =

к/ Б

N О (N О „ }

= Z J (kXQ(r-pw))r«V(r)dlr = -k;

w = 1K

Z IiP (Pw )r(wn)

J

p/

V w = ! J

Здесь учтено, что 0(р-4) = -0(4 - р), и введено обозначение 1 (р) =

= | 0(р-4)фР(4)й14.

кг

При вычислении интегралов, входящих в (Ъа^ )г, интенсивности присоединенных слоев вихрей и источников для согласованности примем кусочно-линейными:

N /

у^(г)= Е ((ум)0Ф0(г) + (ум)1 ф-(г));

г = 1

N , ч

Ча"(г)= )0ф0(г) + )-ф-(г)), гек.

г=1

Коэффициенты (уш )р и (яаи )р зависят от скорости движения центра г-й панели иК (ci), угловой скорости ее вращения ш/ и скорости изменения ее длины Ч:

(уШ )0 = и к (ci)-Т, (у)1 = Ч; (Яа")0 = и к (ci) - П, (яаíí )1 = -LiШi.

Тогда оба интеграла в (Ьац )р можно выразить через ранее введенные величины 1рЯ:

У

N / \

I фр (г) й/г §(к х д(г-4) ) (4) й/4 = к хх(1р,0(уй* )° +1р ру^),

К К ] = 1

N / \

I фр (г) й/г § д(г-4) (4) й/4 = х )0 + !р;1(яаи )}).

к к ]=г

Необходимые соотношения для коэффициентов схемы Т0 получаются из приведенных выше, если рассматривать лишь значения р = 0, я = 0.

Вычисление коэффициентов Т-схем. Представим формулы для вычисления интегралов

1р (р)= I 0(р-4) фр(4) й/4 и 1ря = I фр (г) й/г I д(г -4)ф? (4) й\.

к к к]

Введем вектор А], совпадающий с панелью К] (его орт т]), и векторы 8 и р (рис. 1, а), тогда

10(г) = аи0 х к + Аи°, 11 (г) = аи1 х к + Аи - 1/(2л) т],

где использованы обозначения ) — главное значение аргумента

комплексного числа г):

и0 = ш(т], Т], Т] ) = т], и1 = 1 ш (р + 8, Т], Т] ),

2| А ] |

а = — Z(p, 8), А = —¡п18-1, 2л 2л | р |

ш (а, Ь, с) = (а - Ь)с + (а х Ь) х с, Z(p, 8) = аrg(p - 8 +1 (р 8 к)).

Для двойных интегралов 1рЯ с учетом обозначений (рис. 1, б) справедливы выражения:

10/ = -100= I й!г I д(г -4) й/4 =

к к]

= (щу00 - а2V00 + азV03°) х к + (Xу00 - А2 V°° + АзV°°),

10/ = -1)0= j dlr j Q(r-4)ф1 (4)d\ = ((a! +аэ)у?31 -(«2 +as)v23)*k +

Kj Kj

+ ( ( + À3 )v033 - (^2 + )v03 ) -1/ (4л) | di | Ijj1 = -133 = j dlr j Q(r-4M (r)ф3(4)dl% =

T j,

Ki Kj

= ((ai +a3)v!3 -(a2 +a3)v23 -a3v33 )xk +

+ (i + À3 )v)3 - (À2 + À3 )v23 - À3v33 ) + +1/ (24 л) (| d j | Tj +1 di | t j - 2ш (sb Tj, t j )).

. j-я панель f-я панель^/ ( j-я панель

Точка р _ d,-

с радиусом г

а

Рис. 1. Панели и вспомогательные векторы В приведенных формулах скалярные коэффициенты имеют вид

«1 = ¿(s2, si), «2 ¿(s2, pi), аэ = -1- ¿(pi, Р2), 2л 2л 2л

2л | S2 | 2л | S21 2л | Pi |

и введены векторные множители v00 = ш (si, Tj, т j ), v 2°= © (d ,•, Tj, T j ),

700 _

v00 = © (p2, Ti, T j ),

v13

01 =—rг~| (((pi + si) -T j ) © (si, Tj, T j )-1 Si |2 Tj ) 2| d j|

v03 © (si + P2, T j, T j X 2|dj |

vi3 -

i2|dj ||d

j

( 2(si - © (si-3p2 ,Tj ,T j )) © (si,Tj ,T j )-1 si |2 (si-3p2) ) - ) © (si,Tj ,T j ),

11 i | d/1 ,, , 11 i | d j |

v23 =^тт7,Tj,TjЬ v33 ©(dj,T,T)-i2 | dj | i2 | d/1

Величины ai и Ai или аэ и Аэ не определены, если панели K¡ и Kj являются соседними (имеют общую точку). При этом в Ipq внешние интегралы становятся несобственными; при j = i +1 (когда si = 0, p2 Ф 0) следует принять ai = Ai = 0, а при i = j +1 (когда p2 = 0, si Ф 0) необходимо принять аэ = Аэ = 0. Если i = j, то интегралы становятся сингулярными, их главные значения по Коши формально получаются, если принять все коэффициенты ak = Ak = 0: I00 = Ii = 0; I01 = —= -i / (4 л) di.

Схема TpEM с кусочно-линейным непрерывным представлением решения. Схема Т порождает СЛАУ в 2 раза большей размерности по сравнению со схемой Т0. Учет непрерывности решения на гладких профилях позволяет построить схему, обеспечивающую точность схемы Т при сохранении размерности СЛАУ схемы Т0. В этом случае приближенное решение можно искать в виде непрерывного распределения:

N

y(r)= УУ i ф i (r), (i2)

i = i

где финитные кусочно-линейные базисные функции ф¿ (r), отличные от нуля только на панелях Ki — i и K¡, аналогичны функциям формы первого порядка в методе конечных элементов и могут быть выражены через введенные ранее функции ф0(г) и ф1 (r):

д/ч Г 0,5 ф0—i(r) + ф1—i(r), r е Ki—i, ф i (r) = i

[0,5 ф0(г) — ф}(r), r е Ki, i = i,..., N, K0 = Kn.

Неизвестные {y¿}N=i в (i2) — интенсивности вихревого слоя в началах панелей. Применяя процедуру Галеркина и выполняя регуляризацию, получаем систему вида

' [ A] {IN } у {У} Yf {b} ^ {L}T 0 =

Компоненты блоков [ A] можно выразить через компоненты блоков [Apq ] из схемы Т (индексы суммирования p, q, m и n принимают значения 0 и i; dpq = 0 при p Ф q):

п- = 2 у (—i)p(i—m)+q(i—n)2p+q—2j (( + dpq \

ui) ~ J , J Zj v L> А ч — myu-i — m,j — n^ui — m,j — nf-

У

V R У

vry

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Li -1 + L

i p, q, m, n

Компоненты строки {Ь}т имеют вид Ь/ = (Ь1-1 + Ь/ )/2, а компоненты столбца правой части {Ъ} также выражаются через компоненты столбцов {Ъ0}, {Ъ1} схемы Т1

Ъ = -1+- Е (-1)р(т --)2р -1 и -тЪР-т, р, т = 0,1.

Ьг -1 + Ьг р,т

Для описанной схемы будем использовать обозначение Т рЕМ, имея в виду, что в смысле представления решения она принципиально не отличается от схемы Т-.

Другие варианты численных схем. Выбор в качестве базисных и проекционных функций 5-функций Дирака позволяет построить схему кол-локаций, в которой завихренность из вихревого слоя представляется сосредоточенной в виде точечных вихрей. Такая схема не имеет существенных преимуществ перед широко известной схемой МДВ [10] и даже несколько уступает ей по точности на решениях модельных задач.

При построении дискретного аналога ГИУ можно обеспечивать равенство нулю проекций невязки на семейство функций, отличающихся от базисных. Такой подход соответствует методу Петрова — Галеркина. Так, в [22, 23] описана схема, похожая на схему Т ЕЕм с непрерывным кусочно-линейным решением, однако удовлетворение ГИУ обеспечивается лишь в отдельных точках (метод коллокаций); такая схема основана на ^модели (5). Кроме того, вместо решения СЛАУ с квадратной матрицей, получающейся после регуляризации введением дополнительной переменной, в указанных работах проводится поиск псевдорешения переопределенной системы, что резко повышает сложность алгоритма.

Удовлетворение ГИУ в [24, 25] выполняется интегрально по панелям, но снова в рамках ^модели. Отметим, что в [24] указано на возможность применения Т-модели, но для схемы коллокаций и также с решением переопределенной системы.

Смежные вопросы. При включении Т-схем в алгоритмы расчета нестационарного обтекания профилей возникает проблема преобразования завихренности в форме вихревого слоя на границе профиля в вихревые элементы (частицы), моделирующие распределение завихренности в области течения. Обычно она решается «стягиванием» завихренности с панелей в точечные вихревые частицы. При этом результаты решения модельных задач [26] показывают, что такие вихревые частицы должны иметь сравнительно малые циркуляции, не превосходящие по модулю некоторое

значение Гтах, что означает необходимость образовывать из содержащейся на панели завихренности более чем один вихрь.

Пусть у0, у1 — среднее значение и вариация (изменение) интенсивности вихревого слоя на ;-й панели длиной Ц (для схемы Т0 необходимо принять у} = 0). Если интенсивность вихревого слоя на панели не меняет знак, т. е. | с | < 2, где с = у - / у0, то число образуемых на ней вихревых элементов щ = || у0 | Ц / Гтах , где операция [•] означает округление вверх. Вновь генерируемые вихри получают циркуляции | у0 | Ц / щ и помещаются в точки с радиус-векторами ггк = рг + 0,5(щ _ 1 + щ )Цт; + 5Цпг, к = 1,..., щ,

где pi — радиус-вектор начала г-й панели; т;-, пг — направляющий орт и орт нормали; 8> 0 — малое число, вводимое, чтобы вихри оказались внутри, а не на границе области течения (в практических расчетах можно принимать 5 ~10-4); величины

2к 1 п Щк =—7-, ч , к = 0,..., щ,

щ ((1 - 0,5 сг) + 4 (1 - 0,5 с )2 + 2 к с / щ )

означают доли длины панели, между которыми заключена равная завихренность.

При | сг | > 2, когда интенсивность вихревого слоя на г-й панели меняет знак, описанная процедура выполняется двукратно для частей исходной панели, разделяемых точкой щ* = (с -2)/ (2с) (если у0 = 0, то формально сг =и щ*=1/2). Для первой из полученных «подпанелей» р1 = рг, Ц = щ*Ц, у0,1 = 0,5 (у0 -у1), у1,1 = -2у0,1, с1 = -2; для второй — р/1 = рг +щ*1;т;, I? = (1-щ*)Ц, у?,* = 0,5(у0 +у}), у},11 = 2у0,11, с/1 =2.

Далее вихревые частицы перемещаются в области течения, однако при вычислении поля скоростей вблизи границы профиля вместо влияния образованных на текущем шаге расчета вихрей следует учитывать в (1) влияние породившего их вихревого слоя:

У(г) = Ую + X (кх(г-р„)) +

Щ = 1

N , , N

+ X кхХI?(Г)( + (уй")?)+ X XI?(г)(?й*)?. ; =1 ? } =1 ?

Здесь Г^Р), рщ — циркуляции и радиус-векторы вихревых частиц в области течения, исключая вновь образованные; для модели круглого вихря

Рэнкина 08(г) = г / (2л тах{| г |2, в2}), где в — малый радиус вихря, вводимый для ограничения скоростей, который индуцируется вихрем; ц = 0 для схемы Т0 и ц = 0,1 для схем Т1 и Т ЕЕм.

Конвективные скорости вновь образованных вихревых частиц, расположенных вблизи точек контура гк, с учетом свойств контурных интегралов можно вычислять как

V(r к ) = и к (г к) +1 у(г к) т (гк).

2

Вычислительный эксперимент. Относительную погрешность численного решения у(4) для ГИУ (6), полученного с использованием схем Т0, Т1 и ТрЕм оценим в норме Ь-:

8у = Г* § | у(4) - у*(4) I ^, Г*= 11 у"(4) | .

г к к

Введенная мера погрешности позволяет учесть не только точность определения суммарной завихренности на панелях, но и точность восстановления ее распределения вдоль панелей; поэтому схема Т0 не обеспечит порядок точности выше первого, схемы Т1 и Т ЕЕм — выше второго. Следует учесть, что точное решение у * задано на границе исходного профиля к, численное — на многоугольнике, поэтому для оценки погрешности необходимо спроецировать точное решение на многоугольник, аппроксимирующий профиль (или наоборот). Для рассмотренных ниже тестовых задач точные аналитические решения для интенсивности вихревого слоя на профиле известны [27].

Результаты моделирования обтекания кругового профиля и эллиптического профиля с соотношением полуосей 4 : 1, установленных под углом атаки л /6 к набегающему потоку, разбитых на N панелей близкой длины к, приведены на рис. 2. Вихревой след отсутствует, фактически решается задача моделирования потенциального обтекания при нулевой суммарной циркуляции Г.

При наличии вихревого следа, моделируемого точечными вихрями*, особенно если эти вихри находятся вблизи профиля, точность решения интегрального уравнения несколько снижается. Результаты моделирова-

* Чтобы избежать возникновения неограниченных скоростей, применяют различные способы сглаживания особенности, в результате обычно используется модель круглого вихря Рэнкина или Ламба малого радиуса.

10

10

10"

10

10"

10

1-3

,-5

,-7

100 200

400 800 1600 а

N

100 200

400 800 1600 б

N

Рис. 2. Убывание погрешности при моделировании потенциального обтекания кругового (а) и эллиптического с отношением осей 4 : 1 (б) профилей для схем Т (.), Т1(д) и Т^м (♦)

ния обтекания тех же профилей, но в присутствии одного вихря, находящегося на расстоянии 2 % от радиуса окружности или большой полуоси эллипса при отсутствии набегающего потока, приведены на рис. 3. Согласно рисунку, схемы Т и ТЕЕм с кусочно-линейным представлением интенсивности вихревого слоя обеспечивают второй порядок точности, тогда как схема Т0 с кусочно-постоянным представлением решения — только первый порядок.

10"

10

,-2

10

10"

,-3

д.

I.

д

• 0(h)

д

o(h2)

А

д

„ д

100 200

400 800 1600 а

10"1 д »•

10~2 -

10"3 - /

♦ д 1СГ4 - Oih2)

N 100 200 400

O(h)

д

д

д

800 1600 N

Рис. 3. Убывание погрешности при моделировании течения, индуцированного вихрем вблизи кругового (а) и эллиптического с отношением осей 4 : 1 (б) профилей для схем Т0 (•), Т1(д) и Т^Ем (♦)

Разбиение границы профиля на большее число панелей ведет к увеличению числа обусловленности матриц соответствующих систем (рис. 4), однако они остаются хорошо обусловленными. Отметим, что число обусловленности матрицы линейной системы, возникающей в схеме ТЕЕм, близко к такому для схемы Т0, особенно при сравнительно больших значениях N.

сопсЫ

cond А

N

100 200 400 800 1600 б

N

100 200 400 800 1600 а

Рис. 4. Зависимость числа обусловленности матрицы линейной системы от числа панелей для кругового (а) и эллиптического с отношением осей 4 : 1 (б) профилей для схем Т0 (•), Т1(д) и Т рЕМ (♦)

Несмотря на удвоенную размерность матрицы СЛАУ, использование схемы Т1 и тем более схемы ТрЕм оказывается более «экономичным», поскольку обеспечение той же точности требует многократно меньшего числа панелей, чем в схеме Т0. Число панелей (округленное в большую сторону), на которые требуется разбить границы эллиптических профилей с различным отношением полуосей при использовании численных схем Т0, Т1 и ТрЕм для достижения относительной погрешности 5у = 10-3 и 5у = 10-4 в задачах, аналогичных рассмотренным выше*, приведено в таблице.

При моделировании потенциального обтекания во всех случаях наивысшую точность обеспечивает схема ТрЕм, однако при наличии вихря вблизи профиля она уступает по точности схеме Т1. Наблюдаемое в последнем случае уменьшение требуемого числа панелей при рассмотрении эллипса с отношением полуосей 4 : 1 по сравнению с эллипсами с отношениями 2 : 1 и 1 : 1 можно объяснить изменением структуры решения, которое имеет место при переходе к моделированию обтекания сравнительно тонких профилей.

Число панелей, необходимое для достижения точности 5у при обтекании эллипсов с различным отношением полуосей

Полуоси эллипса Положение вихря 8у = 10 3 для схем 8у = 10 4 для схем

т0 т1 т 1 1 11 FEM Т0 ТЧ Т1 1 11 FEM

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Потенциальное обтекание 1 : 1 | - | 1 600 | 50 | 44 | 15 800 | 160 | 140

* При моделировании потенциального обтекания и обтекания в присутствии вихря, находящегося на расстоянии около 2 % большой полуоси эллипса.

Окончание таблицы

Полуоси эллипса Положение вихря 8у = 10 3 для схем 8у = 10 4 для схем

Т0 Т1 Т1 1 FEM Т0 Т1 Т1 Т FEM

2 : 1 Потени иально( 1 610 обтек 50 ание 44 16 100 320 260

4 : 1 2 400 250 200 23 400 780 630

10 : 1 5 200 920 750 51 800 3 100 2 500

1 : 1 Обтекание (0,50, 0,89) в прису 20 500 тстви 1 200 и вихря 1 900 205 000 3 800 6 100

2 : 1 (0,60, 0,42) 24 300 1 100 1 700 243 000 3 300 5 200

4 : 1 (0,70, 0,20) 20 500 850 1 400 205 000 2 700 4 200

10 : 1 (0,35, 0,109) 27 200 1 300 1 800 272 000 4 200 5 700

Заключение. Для моделирования обтекания гладких профилей в вихревых методах вычислительной гидродинамики разработаны Т-схемы первого и второго порядка точности для численного решения ГИУ относительно интенсивности вихревого слоя, выражающего граничное условие прилипания. Получены необходимые расчетные формулы, которые позволяют включать созданные схемы в новые и существующие программные коды, реализующие метод вязких вихревых доменов, а также другие модификации вихревых методов.

ЛИТЕРАТУРА

[1] Cottet G.-H., Koumoutsakos P.D. Vortex methods: theory and practice. Cambridge Univ. Press, 2000.

[2] Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Cambridge Univ. Press, 1991.

[3] Branlard E. Wind turbine aerodynamics and vorticity-based methods. Cham, Springer, 2017. DOI: https://doi.org/10.1007/978-3-319-55164-7

[4] Головкин М.А., Головкин В.А., Калявкин В.М. Вопросы вихревой гидромеханики. М., ФИЗМАТЛИТ, 2009.

[5] Mimeau C., Mortazavi I. A review of vortex methods and their applications: from creation to recent advances. Fluids, 2021, vol. 6, no. 2, art. 68.

DOI: https://doi.org/10.3390/fluids6020068

[6] Leonard A. Vortex methods for flow simulation. J. Comput. Phys., 1980, vol. 37, iss. 3, pp. 289-335. DOI: https://doi.org/10.1016/0021-9991(80)90040-6

[7] Sarpkaya T. Computational methods with vortices — the 1988 Freeman scholar lecture. J. Fluids Eng., 1989, vol. 111, iss. 1, pp. 5-52.

DOI: https://doi.org/10.1115/L3243601

[8] Дынникова Г.Я. Аналог интегралов Бернулли и Коши — Лагранжа для нестационарного вихревого течения идеальной несжимаемой жидкости. Изв. РАН. МЖГ, 2000, № 1, с. 31-41.

[9] Дынникова Г.Я. О присоединенной массе в модели вязкой несжимаемой жидкости. ДАН, 2019, т. 488, № 5, с. 493-497.

DOI: https://doi.org/10.31857/S0869-56524885493-497

[10] Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М., Янус, 1995.

[11] Дынникова Г.Я. Лагранжев подход к решению нестационарных уравнений Навье — Стокса. ДАН, 2004, т. 399, № 1, с. 42-46.

[12] Kempka S.N., Glass M.W., Peery J.S., et al. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA report SAND-96-0583, 1996.

[13] Быховский Э.Б., Смирнов Н.В. Об ортогональном разложении пространства вектор-функций, квадратично суммируемых по заданной области, и операторах векторного анализа. Труды МИАН СССР, 1960, т. 59, с. 5-36.

[14] Morino L. Helmholtz decomposition revisited: vorticity generation and trailing edge condition. Comput. Mech., 1986, vol. 1, no. 1, pp. 65-90.

DOI: https://doi.org/10.1007/BF00298638

[15] Wu J.C., Thompson J.F. Numerical solutions of time-dependent incompressible Navier — Stokes equations using an integro-differential formulation. Comput. Fluids, 1973, vol. 1, iss. 2, pp. 197-215. DOI: https://doi.org/10.1016/0045-7930(73)90018-2

[16] Марков В.В., Сизых Г.Б. Эволюция завихренности в жидкости и газе. Изв. РАН. МЖГ, 2015, № 2, с. 8-15.

[17] Ogami Y., Akamatsu T. Viscous flow simulation using the discrete vortex model — the diffusion velocity method. Comput. Fluids, 1991, vol. 19, no. 3-4, pp. 433-441. DOI: https://doi.org/10.1016/0045-7930(91)90068-S

[18] Бэтчелор Дж. Введение в динамику жидкости. М., Мир, 1973.

[19] Марчевский И.К. Разработка и реализация Т-схем численного решения граничных интегральных уравнений в математических моделях вихревых методов вычислительной гидродинамики. Дис. ... д-ра физ.-мат. наук. М., МГТУ им. Н.Э. Баумана, 2021.

[20] Кузьмина К.С., Марчевский И.К., Морева В.С. Определение интенсивности вихревого слоя при моделировании вихревыми методами обтекания профиля потоком несжимаемой среды. Математическое моделирование, 2017, т. 29, № 10, с. 20-34.

[21] Кузьмина К.С., Марчевский И.К. О влиянии вихревого слоя и точечных вихрей при приближенном решении граничного интегрального уравнения в двумерных вихревых методах вычислительной гидродинамики. Прикладная математика и механика, 2019, т. 83, № 3, с. 495-508.

DOI: https://doi.org/10.1134/S003282351903010X

[22] Morgenthal G. Aerodynamic analysis of structures using high-resolution vortex particle methods. Ph.D. thesis. Univ. of Cambridge, 2002.

[23] Milani G. Concepts of adaptivity for vortex particle methods and applications to bluff body aerodynamics. Ph.D. thesis. Bauhaus-Universität Weimar, 2018.

[24] Lin H., Vezza M. Implementation of a vortex method for the prediction of separated incompressible flows. Aero Report no. 9425. Univ. of Glasgow, 1994.

[25] Taylor I., Vezza M. Prediction of unsteady flow around square and rectangular section cylinders using a discrete vortex method. J. Wind. Eng. Ind. Aerodyn., 1999, vol. 82, iss. 1-3, pp. 247-269. DOI: https://doi.org/10.1016/S0167-6105(99)00038-0

[26] Kuzmina K., Marchevsky I., Soldatova I., et al. On the scope of Lagrangian vortex methods for two-dimensional flow simulations and the POD technique application for data storing and analyzing. Entropy, 2021, vol. 23, no. 1, art. 118.

DOI: https://doi.org/10.3390/e23010118

[27] Kuzmina K., Marchevsky I., Ryatina E. Exact solutions of boundary integral equation arising in vortex methods for incompressible flow simulation around elliptical and Zhukovsky airfoils. J. Phys.: Conf. Ser, 2019, vol. 1348, art. 012099.

DOI: https://doi.org/10.1088/1742-6596/1348/1/012099

Марчевский Илья Константинович — д-р физ.-мат. наук, доцент, профессор кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1); старший научный сотрудник ИСП РАН (Российская Федерация, 109004, Москва, ул. А. Солженицына, д. 25).

Сокол Ксения Сергеевна — канд. физ.-мат. наук, доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1); старший научный сотрудник ИСП РАН (Российская Федерация, 109004, Москва, ул. А. Солженицына, д. 25).

Измайлова Юлия Андреевна — студент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1); старший лаборант ИСП РАН (Российская Федерация, 109004, Москва, ул. А. Солженицына, д. 25).

Просьба ссылаться на эту статью следующим образом:

Марчевский И.К., Сокол К.С., Измайлова Ю.А. Т-схемы для математического моделирования генерации завихренности на гладких профилях в вихревых методах. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2022, № 6 (105), с. 33-59. DOI: https://doi.org/10.18698/1812-3368-2022-6-33-59

T-SCHEMES FOR MATHEMATICAL MODELLING OF VORTICITY GENERATION ON SMOOTHS AIRFOILS IN VORTEX PARTICLE METHODS

I.K. Marchevskii1' 2 iliamarchevsky@bmstu.ru

K.S. Sokol1' 2 kuz-ksen-serg@yandex.ru

Yu.A. Izmailova1' 2 yulia.izmailova@mail.ru

1 Bauman Moscow State Technical University, Moscow, Russian Federation

2 Ivannikov Institute for System Programming, Russian Academy of Sciences, Moscow, Russian Federation

Abstract

New numerical schemes are proposed for solving the boundary integral equation that arises in CFD vortex particle methods of when simulating a plane flow around smooth airfoils. They are based on considering the 2-nd kind integral equation with respect to vortex sheet intensity with a bounded or absolutely integrable kernel instead of traditionally solved singular integral equations of the 1-st kind with Hilberttype singularity. To solve it, the Galerkin approach is used. It is shown that when approximating the airfoil boundary with a polygon, it is possible to develop schemes of the 1-st and 2-nd order of accuracy, considering a piecewise-constant or piecewise-linear (discontinuous or continuous) distribution of the solution along the panels. The necessary formulae are presented for calculating the components of the matrix and the right-hand side of the system of linear algebraic equations, that is a discrete analogue of the integral equation. They are suitable for modelling of the vorticity generation when simulating the flow around either single airfoil or system of airfoils, including moving and/or deformable ones. The developed schemes can be used in the framework of the viscous vortex domains method as well as other modifications of vortex particle methods, since they only concern the convective velocities of the flow near the airfoil and are not related to methods for modeling viscous diffusion of vorticity

Keywords

Vortex particle methods, vorticity generation, boundary integral equation, Galerkin method, the second order of accuracy

Received 19.07.2022 Accepted 05.08.2022 © Author(s), 2022

This research was supported by the Ministry of Science and Higher Education of the Russian Federation (agreement no. 075-15-2020-808)

REFERENCES

[1] Cottet G.-H., Koumoutsakos P.D. Vortex methods: theory and practice. Cambridge Univ. Press, 2000.

[2] Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Cambridge Univ. Press, 1991.

[3] Branlard E. Wind turbine aerodynamics and vorticity-based methods. Cham, Springer, 2017. DOI: https://doi.org/10.1007/978-3-319-55164-7

[4] Golovkin M.A., Golovkin V.A., Kalyavkin V.M. Voprosy vikhrevoy gidromekhaniki [Issues of vortex hydromechanics]. Moscow, FIZMATLIT Publ., 2009.

[5] Mimeau C., Mortazavi I. A review of vortex methods and their applications: from creation to recent advances. Fluids, 2021, vol. 6, no. 2, art. 68.

DOI: https://doi.org/10.3390/fluids6020068

[6] Leonard A. Vortex methods for flow simulation. J. Comput. Phys., 1980, vol. 37, iss. 3, pp. 289-335. DOI: https://doi.org/10.1016/0021-9991(80)90040-6

[7] Sarpkaya T. Computational methods with vortices — the 1988 Freeman scholar lecture. J. Fluids Eng., 1989, vol. 111, iss. 1, pp. 5-52.

DOI: https://doi.org/10.1115/L3243601

[8] Dynnikova G.Ya. An analog of the Bernoulli and Cauchy — Lagrange integrals for a time-dependent vortex flow of an ideal incompressible fluid. Fluid Dyn., 2000, vol. 35, no. 1, pp. 24-32. DOI: https://doi.org/10.1007/BF02698782

[9] Dynnikova G.Ya. On the added mass in a viscous incompressible fluid. Doklady Akademii nauk, 2019, vol. 488, no. 5, pp. 493-497 (in Russ.).

DOI: https://doi.org/10.31857/S0869-56524885493-497

[10] Lifanov I.K. Metod singulyarnykh integralnykh uravneniy i chislennyy eksperiment [Methods of singular integral equations and numerical experiment]. Moscow, Yanus Publ., 1995.

[11] Dynnikova G.Ya. Lagrange method for Navier — Stokes equations solving. Doklady Akademii nauk, 2004, vol. 399, no. 1, pp. 42-46 (in Russ.).

[12] Kempka S.N., Glass M.W., Peery J.S., et al. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA report SAND-96-0583, 1996.

[13] Bykhovskiy E.B., Smirnov N.V. Orthogonal decomposition of the space of vector functions square-summable on a given domain, and the operators of vector analysis. Trudy MIAN SSSR, 1960, vol. 59, pp. 5-36 (in Russ.).

[14] Morino L. Helmholtz decomposition revisited: vorticity generation and trailing edge condition. Comput. Mech., 1986, vol. 1, no. 1, pp. 65-90.

DOI: https://doi.org/10.1007/BF00298638

[15] Wu J.C., Thompson J.F. Numerical solutions of time-dependent incompressible Navier — Stokes equations using an integro-differential formulation. Comput. Fluids, 1973, vol. 1, iss. 2, pp. 197-215. DOI: https://doi.org/10.1016/0045-7930(73)90018-2

[16] Markov V.V., Sizykh G.B. Vorticity evolution in liquids and gases. Fluid Dyn., 2015, vol. 50, no. 2, pp. 186-192. DOI: https://doi.org/10.1134/S0015462815020027

[17] Ogami Y., Akamatsu T. Viscous flow simulation using the discrete vortex model — the diffusion velocity method. Comput. Fluids, 1991, vol. 19, no. 3-4, pp. 433-441. DOI: https://doi.org/10.1016/0045-7930(91)90068-S

[18] Batchelor G.K. An introduction to fluid dynamics. Cambridge Univ. Press, 1973.

[19] Marchevskii I.K. Razrabotka i realizatsiya T-skhem chislennogo resheniya granich-nykh integralnykh uravneniy v matematicheskikh modelyakh vikhrevykh metodov vychislitelnoy gidrodinamiki. Dis. d-ra fiz.-mat. nauk [Development and implementation of T-schemes for numerical solution of boundary integral equations in mathematical models of vortex methods of computational fluid dynamics. Dr. Phys.-Math. Sc. Diss.]. Moscow, BMSTU, 2021 (in Russ.).

[20] Kuzmina K.S., Marchevskii I.K., Moreva V.S. Vortex sheet intensity computation in incompressible flow simulation around an airfoil by using vortex methods. Math. Models Comput. Simul., 2018, vol. 10, no. 3, pp. 276-287.

DOI: https://doi.org/10.1134/S2070048218030092

[21] Kuzmina K.S., Marchevskii I.K. On the calculation of the vortex sheet and point vortices effects at approximate solution of the boundary integral equation in 2D vortex methods of computational hydrodynamics. Fluid Dyn., 2019, vol. 54, no. 7, pp. 9911001. DOI: https://doi.org/10.1134/S0015462819070103

[22] Morgenthal G. Aerodynamic analysis of structures using high-resolution vortex particle methods. Ph.D. thesis. Univ. of Cambridge, 2002.

[23] Milani G. Concepts of adaptivity for vortex particle methods and applications to bluff body aerodynamics. Ph.D. thesis. Bauhaus-Universität Weimar, 2018.

[24] Lin H., Vezza M. Implementation of a vortex method for the prediction of separated incompressible flows. Aero Report no. 9425. Univ. of Glasgow, 1994.

[25] Taylor I., Vezza M. Prediction of unsteady flow around square and rectangular section cylinders using a discrete vortex method. J. Wind. Eng. Ind. Aerodyn., 1999, vol. 82, iss. 1-3, pp. 247-269. DOI: https://doi.org/10.1016/S0167-6105(99)00038-0

[26] Kuzmina K., Marchevsky I., Soldatova I., et al. On the scope of Lagrangian vortex methods for two-dimensional flow simulations and the POD technique application for data storing and analyzing. Entropy, 2021, vol. 23, no. 1, art. 118.

DOI: https://doi.org/10.3390/e23010118

[27] Kuzmina K., Marchevsky I., Ryatina E. Exact solutions of boundary integral equation arising in vortex methods for incompressible flow simulation around elliptical and Zhukovsky airfoils. J. Phys.: Conf. Ser., 2019, vol. 1348, art. 012099.

DOI: https://doi.org/10.1088/1742-6596/1348/1Z012099

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Marchevskii I.K. — Dr. Sc. (Phys.-Math.), Assoc. Professor, Professor of Department of Applied Mathematics, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation); Senior Researcher, Ivanni-kov Institute for System Programming, Russian Academy of Sciences (A. Solzhenitsi-na ul. 25, Moscow, 109004 Russian Federation).

Sokol K.S. — Cand. Sc. (Phys.-Math.), Assoc. Professor, Department of Applied Mathematics, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation); Senior Researcher, Ivannikov Institute for System Programming, Russian Academy of Sciences (A. Solzhenitsina ul. 25, Moscow, 109004 Russian Federation).

Izmailova Yu.A. — Student, Department of Applied Mathematics, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation); Senior Laboratory, Ivannikov Institute for System Programming, Russian Academy of Sciences (A. Solzhenitsina ul. 25, Moscow, 109004 Russian Federation).

Please cite this article in English as:

Marchevskii I.K., Sokol K.S., Izmailova Yu.A. T-schemes for mathematical modelling of vorticity generation on smooths airfoils in vortex particle methods. Herald of the Bauman Moscow State Technical University, Series Natural Sciences, 2022, no. 6 (105), pp. 33-59 (in Russ.). DOI: https://doi.org/10.18698/1812-3368-2022-6-33-59

i Надоели баннеры? Вы всегда можете отключить рекламу.