Особенности вихревой схемы при моделировании дозвукового отрывного обтекания с полубесконечным эквивалентным телом Текст научной статьи по специальности «Физика»

УДК 531.6.011.32:532.582.4:517.958

DOI 10.18698/2309-3684-2017-4-7391

Особенности вихревой схемы при моделировании дозвукового отрывного обтекания с полубесконечным эквивалентным телом

© В.Н. Тимофеев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: математическое моделирование, дозвуковое отрывное обтекание, эквивалентное тело, вихревая схема, метод дискретных вихрей, донное давление

Введение. Решение задач дозвукового отрывного обтекания тел является актуальным для широкого круга практических приложений, включающих в себя вопросы проектирования и эксплуатации транспортных средств, архитектурных сооружений, промышленных установок, а также объектов самого разного назначения и уровня сложности. Для получения аэродинамических характеристик при наличии отрыва потока необходимы большие вычислительные мощности даже при использовании моделей не самого высокого уровня сложности, в том числе моделей, основанных на теории идеальной среды [1-6].

Для снижения затрат вычислительных ресурсов при математическом моделировании пространственного обтекания тел дозвуковым потоком газа с отрывом, локализованным в окрестности линии донного среза, была разработана методика, основные положения которой изложены в работах [7, 8]. В соответствии с данной методикой на одном из этапов полномасштабной реализации концепции вязко-невязкого взаимодействия [9] или при приближенном решении задачи о дозвуковом отрывном обтекании [10, 11] будем считать, что поверхность переднего участка Е(1) эквивалентного тела совпадает

с поверхностью Е (0) обтекаемого тела, а поверхность Е(2) хвостового участка, моделирующего влияние спутного следа, образующего-

ся вниз по потоку за линией Ьс стыковки уходит в бесконечность,

непрерывно деформируясь.

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

участка начинают уменьшаться. Считаем, что указанное уменьшение происходит вплоть до граничной линии Ь,. Вниз по потоку от этой

линии форма и размеры поперечного сечения не меняются, а образующие хвостового участка эквивалентного тела становятся параллельными вектору скорости V набегающего потока. Часть Е (21) поверхности хвостового участка, расположенная между линиями Ьс

и Ь,, имеет конечную протяженность. Другой участок Е(22) поверхности эквивалентного тела, расположенный вниз по потоку от граничной линии Ь,, также является составной частью поверхности

Е(2) хвостового участка (назовем его полубесконечным). Полная поверхность Е эквивалентного тела представляет собой объединение поверхностей переднего и хвостового участков Е(1)и Е(2) или

Е(1) и Е(21) и Е(22).

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

обтекаемого тела. Границей каждого из таких четырехугольников является П-образная ломаная линия.

На следующем этапе численного моделирования происходит формирование вихревой схемы течения (рис. 1), т. е. размещение на границах многоугольников Ък набора вихревых нитей Ьк (к = 1,...,Ы,).

На части Е (21) поверхности хвостового участка размещаются вихревые многоугольники конечных размеров. Границы удлиненных четы-

рехугольников на полубесконечном участке Е(22) заменяем на П-образ-ные вихревые нити. В данной статье рассмотрим некоторые особенности формирования вихревой схемы течения для эквивалентного тела с полубесконечным участком.

Рис. 1. Вихревая схема течения: Су — у-я контрольная точка, расположенная на поверхности

Математическая модель, вихревая схема течения и алгоритм

решения. Изучим установившееся обтекание тела при умеренных дозвуковых скоростях набегающего потока газа. Среду будем считать несжимаемой, однородной и невесомой жидкостью, динамика которой описывается уравнением неразрывности и уравнениями движения в форме Эйлера [12]. Требуется выполнить граничное условие непротекания на поверхности эквивалентного тела и граничное условие затухания возмущений на бесконечности. Для построения используется изображенная на рис. 1 прямоугольная декартова система координат Oyz, связанная с обтекаемым телом. Предположим, что вне поверхности эквивалентного тела течение является потенциальным, а на самой поверхности Е существуют предельные значения потенциала, его градиента и нормальной производной [13]. В любой точке М0(х0, уо, z0), расположенной вне указанной поверхности, разность V(М0) - , называемую вектором возмущенной скорости, считаем равной градиенту потенциала ф(Мо) возмущенных скоростей, здесь V (Мо) — вектор скорости потока в точке Мо,

поэтому найдем его как сумму вектора скорости набегающего потока и градиента потенциала возмущенных скоростей

V (Мо ) = Va¡>+V ф(Мо), (1)

^ ^ д ^ д ^ д где V = г ---ъ ]--ъ к--оператор Гамильтона [14].

дхо

дУо дzо

В принятых допущениях с учетом свойств упомянутых выше предельных значений потенциал возмущенных скоростей должен быть решением внешней задачи Неймана для уравнения Лапласа:

Аф(М0) = 0, М0 £Е; дф(М 0 )

дп

= • п

(М 0 )

М0 еЕ;

lim ф(М0 ) = 0;

lim Уф(М0 ) = 0,

(2)

(3)

(4)

(5)

где А =

д2

д2

д2

дх0 дУ0 дг0

— оператор Лапласа; п(М0)|Е — орт

вектора нормали к поверхности 2 в точке М0; Я — модуль радиус-вектора точки М0.

Решение задачи (2)-(5) найдем в виде потенциала двойного слоя

ф(М0 ) = ± я (М )йа

4л

(6)

Здесь М^Х0, У0, ^0) — точка, в которой вычисляются потенциал возмущенных скоростей и вектор скорости потока; п(М) и g (М) — орт вектора нормали к поверхности 2 и поверхностная плотность потенциала двойного слоя, определяемые в точке М(х, у, г), принадлежащей элементу поверхности й 2 с площадью йъ; г — вектор, соединя-

I 2 2 2

т\ = ^(х - х0) + (у -у0) + (г - г0) —

модуль вектора г.

Определим поверхностную плотность потенциала двойного слоя методами численного моделирования. Поверхность 2 эквивалентного тела заменим конечным числом многоугольников (панелей) Ек (к = 1,...,N3), где N3 — суммарное количество аппроксимирующих многоугольников. Число сторон к-го многоугольника Ек обозначим через Тк. На каждой из панелей Ек поверхностную плотность потенциала двойного слоя g(M) считаем равной постоянной величине gk (к = 1,...,Ns).

Учтем [15], что градиент потенциала двойного слоя, размещенного на панели Ек с постоянной поверхностной плотностью gk, равен

Е

вектору скорости, индуцированной замкнутой вихревой нитью Ьк, которая расположена на границе дЪк указанной панели и имеет постоянную циркуляцию

Гк , к = 1,..., М3. (7)

Скорость, создаваемую замкнутой вихревой нитью, определим с помощью формулы Био — Савара.

На основании свойства линейности градиента, аддитивности и линейности поверхностного интеграла (6) и упомянутого выше равенства преобразуем формулу (1) к виду

V(Мо) = Гкщ (Мо).

(8)

к=1

Здесь тк (М0) — векторы функции скорости, индуцированной к-й вихревой нитью Ьк в точке М0(х0, у0, г0). Каждый из этих векторов определяется в соответствии с формулой Био — Савара:

(Мо ) = £ \ ^ (9)

Ьк

т

где дя — элемент вихревой линии Ьк; г — вектор, соединяющий точку М(х, у, z), принадлежащую элементу ds, с точкой М0;

I —I /22 2 —

г = г =у] (х - х0) + (у - уо) + (г - г0) — модуль вектора г.

Рассмотрим вихревой многоугольник Ьк (к = 1,..., Ы8), расположенный на границе панели Ък, как объединение вихревых отрезков

Тк

Ьк1 (у = 1,..., Тк), т. е. представим его в виде Ьк = и Ьку. В результате,

У=1

например, при Тк = 4, получается вихревой четырехугольник

Ьк = и Ьк:, изображенный на рис. 2.

У=1

В силу аддитивности криволинейного интеграла в правой части соотношения (9) вектор тк (М0)

функции скорости вихревого многоугольника найдем по формуле

- (иж \ 1 Г ^хг т (М0 /— •

у=1 Ьи

Рис. 2. Вихревой четырехугольник Ьк

4

Введем векторы функции скорости, индуцированной в точке М0 вихревым отрезком Ь], определяемые равенствами

& х г

4л

и

к]

(М0 ) =

] = l, Тк.

(10)

' Ьк]

После этого вектор функции скорости вихревого многоугольника считается равным сумме векторов функций скорости составляющих

тк

его вихревых отрезков: ик (М0 ) = ^ ик] (М0 ).

]=1

Поиск векторов ик] (М0) осуществляется после формирования вихревой схемы течения, в ходе которой устанавливается положение в пространстве вихревых многоугольников Ь] (к = 1,...,) путем

задания координат (с1к], у 1к)-, г11д-) точек начала Л1к]- и координат

(х2к], У2к],72к] ) точек конца А2] для всех вихревых отрезков Ь] (к = 1,...,N; ] = 1,...,Тк). Затем определяются координаты векторов г10к] и г20к] с началами в точках А1к], А2к] и с концами в точке М0, а также координаты векторов г12к], имеющих начало в точке А1к] и концы в точке А2] (к = 1,...,; ] = 1,...,Тк). В дальнейшем во избежание громоздкости индексы к] или схожие индексы у всех или у части объектов, с которыми осуществляется оперирование, записаны не будут. Например, на рис. 3 точки А1к] и А2к] изображены как

точки А1 и А 2, а векторы г}

12 к]

г10 к] и г20 к]

как векторы г12, г

10

и г

20

Рис. 3. Вихревой отрезок Ь:

^12 , ?10, ?20 — оРты векторов ^ ^

- -(2)

г 20; **10 — предельное положение орта

?10 для вихревого отрезка второго вида;

-(1)

20

— предельное положение орта для вихревого отрезка первого вида

20

После параметризации отрезка Ьк] интегрированием по формуле

(10) получим вектор функции скорости вихревого отрезка в следующем виде [16]:

w

<(М)=^-Г-ГТ11-^ , (Ц)

4л г 2 г 2 -г12г10

"(г12 • >10 )

гу> • гт ^ • г-><

г10 г20

где г12, г10 и г20 — модули соответствующих векторов. Здесь индекс к] вихревого отрезка Ь] указан лишь у вектора Wkj (М0) и не

записан у всех векторов и их модулей в правой части соотношения.

Как было отмечено во введении, полубесконечный участок эквивалентного тела моделируется П-образными вихревыми нитями. Каждый из таких П-образных вихрей образован тремя вихревыми отрезками, два из них оказываются полубесконечными.

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

а конец устремлен в бесконечность А2^ю. У вихревого отрезка второго вида, наоборот, положение А 2 конца остается неизменным, а начало — А1 ^го. Направление от начала к концу для полубесконечных отрезков обоих видов характеризуется ортом ^.

Для выполнения предельного перехода рассмотрим вихревой отрезок конечной длины. Введем орты соответствующих векторов:

- _ 1 г - _ 1 г - _ 1 г -12 _ г12, -10 _ г10, -20 _ г20. г12 г10 г20

Соотношение (11) преобразуем к одной из двух эквивалентных форм

w

к

(М0 )_ ^ /12*-1° ^ ((2 • -10 --12 • -20 ); (12)

4л

г10 (1 - ^12• ?10)2)

w

к]

(М0)_ ^ , ^Х^ ^ (12 • -0 - -2 • -20 ). (13)

4л

г20 ( 12 • -20) )

В соотношениях (12), (13) были выявлены параметры, которые при упомянутом выше предельном переходе сохраняют конечные значения, и определены предельные значения остальных параметров. Для полубесконечного вихревого отрезка первого вида Ь^] с фиксированным

началом А1 остаются конечными вектор г 10 и его модуль г 10. Условие А2 ^да устремления конца такого отрезка в бесконечность означает, что г20 ^ да. Следовательно, для орта 52о предельное положение ^

совпадает с ортом -512.

Вектор функции скорости полубесконечного вихревого отрезка первого вида получим из формулы (12) подстановкой вместо орта ^о его предельного положения :

(М0) = -^ >Х^ ^ (1 + 512-?10)

4- гш (1 - (512• 5Ш)2)

и последующим преобразованием к виду

< (мо )=4- (Ч )). (14)

к 4- г1о (((512 • 51о))

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

А1 ^да позволяет получить для орта предельное положение

51(о) = 512, а для вектора функции скорости из формулы (13) выражение

у Мо )=4- (А)у (15)

4- г2о I1 + (512 • 52о))

Вихревая схема течения формируется так, чтобы в каждом П-образном вихре ЬПк = ЬП к ЬП к 2 ^ Ьик 3 первым следовал полубесконечный вихревой отрезок ЬП к 1 второго вида, затем вихревой отрезок ЬПк 2 конечной длины, а третьим — полубесконечный вихревой отрезок ЬПк 3 первого вида (рис. 4). Геометрические характеристики П-образного вихря полностью определяются координатами точек начала А П1к (х П1к, У П1к, г П1к ) и конца А П2к (х П 2к, У П 2к, г П 2к ) втоГо-

го вихревого отрезка ЬПк 2, а также ортом , задающим направление, по которому конец третьего отрезка Ьж3 устремляется в бесконечность.

Для того чтобы найти в точке Мо( Хо, Уо, ¿о) векторы функции скорости, индуцированной к-м П-образным вихрем Ьш, вычислим координаты векторов, гП1ок и гП2ок с началами в точках АП1к и Ашк и кон-

цами в точке М0, а также координаты векторов , имеющих начала в точках Атк и концы в точках Ашк. После этого найдем модули

ГП10к , ГП20к, ГП12к и оРтЫ ?П10к > %20к> ?П12к укаЗЭЖЫХ выше векторов. Первый вихревой отрезок £щ;1 П-образного вихря Ьт, являющийся вихревым отрезком второго вида, в первоначальных обозначениях рассмотрим как вихревой отрезок Ь^Щ) с началом и

охарактеризуем параметрами ?12к1 =-?Пк; гП20к1 = гП10к; ?20к1 = %10к.

ШЗ

Рис. 4. П-образный вихрь £

Пк

В соответствии с формулой (15) для вектора функции скорости, индуцированной первым вихревым отрезком ¿ПИ' получилось следующее выражение:

1 х 5П10

Пк1 (М 0 ) = -

4ж ГП10(1 -и' -П10)

Здесь в правой части не указан индекс к для параметров П-образного

вихря £Пк.

Аналогичным образом находились векторы функций скоростей, индуцированных вторым и третьим вихревыми отрезками П-образного вихря Ц^. Для второго вихревого £]Пк 2 отрезка, имеющего конечную длину, справедливы соотношения:

-12к2 = %12к, г10к2 = гП10к, -10к2 = %10к, -20к2 = -П20к

и, как следует из формулы (12), равенство

Пк 2

(М0 ) =

1

-П12 х 5П10

4л

ГП10 ( (^Ш2' -П10) )

(

П12' 5П10 -П12' -П20

Для третьего вихревого отрезка Ьщ3, рассматривающегося как вихревой отрезок первого вида (в первоначальных обозначениях Ь^ с

конц°м A 2к3 ^ю), из равенств ^3 = ; гШ3 = гП20к ; S10k3 = ^20*

в соответствии с формулой (14) получим выражение вектора функции скорости

w

(Мо ) =

1 ?П X £

П А П20

7 4л ГП20(1 ' ,П20)

В результате вектор функции скорости, индуцированной П-образ-ным вихрем Ьш, как сумма WПk (М0 ) = WПk1 (М0 ) + WПk2 (М0 ) +

+ WПk3 (М0) представляется в достаточно простом и удобном для компьютерных вычислений виде Г

w

Пк

(М0 ) =

1 4л

5П12 х 5П10

ГП10 ( (,П12' *%10)2)

5П х 5П10 , 5П х 5П20

("^П12' ^10 ^12' ,П20 )

ГП10 (1 ' ^10 ) ГП20 (1 ^' ,П20 ) ,

После определения векторов функций скоростей, индуцированных всеми вихревыми многоугольниками и П-образным вихрями, используем векторное равенство (8) для нахождения циркуляций Г*; у = 1,..., Ы,. На некотором множестве контрольных точек Су, у = 1,...,Ы, где N — суммарное количество контрольных точек, требуется выполнение граничного условия непротекания поверхности эквивалентного тела. Контрольные точки расположены на поверхностях панелей Еу, находящихся вверх по потоку от линии Ь, и имеющих

конечный диаметр. На панелях, размещающихся вниз по потоку от линии Ь,, контрольных точек не было. Радиус-вектор контрольной

точки Су, принадлежащей у-й панели, находим как среднее арифметическое радиусов-векторов гу1, соответствующих 1-м вершинам у-го вихревого Ту -угольника:

Г Ту \

^ 1

гс =—

су Т

у

Е Гу1

I=1

у = 1,..., N.

Циркуляции П-образных вихрей, которые находятся вниз по потоку от линии Ь,, считаем равными циркуляциям вихревых многоугольни-

ков, примыкающим к ним с противоположной стороны указанной линии. Так, если вниз по потоку от линии Ь. в кольцевом слое размещается Ыг П-образных вихрей Ьщ (к = N. -Ыг +1,...,N.), а в кольцевом слое, примыкающем к линии Ь. вверх по потоку, Nr вихревых многоугольников Ьк (к = N. - 2Nr +1,..., N. - Nr) конечного диаметра, то выполняются равенства Гк+^ = Гк, к = N. - 2Nr +1,..., N. - Nr. При

этом суммарное количество контрольных точек N равно числу N. - Nr, поэтому для определения циркуляций векторное равенство (8) заменим соотношением

N

V(Мо) = ^+ХГкЩ (Мо). (16)

к=1

В нем выполняется следующее преобразование векторов функции скорости вихревого многоугольников:

К (Мо ) = щ (Мо), к = 1,..., N - Nr;

щ (Мо) = щк (Мо) + щПк+мг (Мо), к = N - Nr +1,..., N.

После такого преобразования граничное условие непротекания (3) поверхности эквивалентного тела на множестве контрольных точек запишем в виде:

Х Гк (() • я (С,)) = - V. ■ Я(С,), у = 1,..., N.

к=1

Данное соотношение представляет собой систему линейных ал-

N

гебраических уравнений Х ахк Г к = К, ,= 1,..., N, в которой

к=1

= ) • Я (С,), у = 1,..., N; к = 1,..., N

6v=-Vю• Я (С,), У = 1,..., N.

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

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

Ау (Су) = М^ т1(Су) +^^ ^2 (Су), у = 1,..., N. ох1 дг2

дg (Су) (Су)

Здесь частные производные -и —определяются по

дт1 дт2

направлениям двух неколлинеарных единичных векторов Т1(Су) и т2(Су), принадлежащих касательной плоскости, проходящей через контрольную точку Су, и вычисляются численным дифференцированием по значениям циркуляций у-го и соседних с ним вихревых многоугольников с учетом равенства (7).

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

N 1

V (Су) = Ую +Е ГkWl (Су)+- АУ (Су). к=1 2

Статическое давление в любой точке потока, в том числе и контрольной, определим с помощью интеграла Бернулли:

+1 Рю (( -V2), (17)

Р = Р- 2

где рю, Ую — статическое давление и модуль вектора скорости набегающего потока соответственно.

V

и его без-

После вычисления модуля вектора скорости V = размерного значения V / Vю коэффициент давления ср = 2(р-рю)/ /(Р^2) находим из интеграла Бернулли (17) по формуле ср =

= 1-^ / V» )2.

Результаты математического моделирования. Исследовано обтекание цилиндрических тел с головной частью оживальной формы. Характерным линейным размером считался диаметр миделевого

*

сечения цилиндрической части ат. Безразмерные значения длин

и координат получаются после деления соответствующей размерной

*

величины на диаметр ёт. Геометрическая форма обтекаемого тела полностью характеризуется безразмерными значениями длины I & головной части и длины Iс цилиндрической части. Суммарная длина обтекаемого тела, равная I +1 с, обозначается через ¡г.

При нулевых углах атаки и скольжения эквивалентное тело оказывается осесимметричным. Форма поверхности Х(21) хвостового участка эквивалентного тела определяется длиной I. этого участка, равной расстоянию по оси симметрии между линиями Ьс и Ь., и законом изменения его диаметра в зависимости от продольной координаты й = й(х). Проведение математического моделирования с использованием рекомендаций и приемов, описанных в работах [7, 8], позволяет найти длину ¡., минимальный диаметр й. сечений, расположенных вниз по потоку от линии Ь . , и распределения безразмерной скорости и коэффициента давления по поверхности исследуемого тела в зависимости от числа Рейнольдса и геометрических параметров обтекаемого тела. Число Рейнольдса Re определено по скорости набегающего потока и по суммарной длине обтекаемого тела ¡г. На рис. 5 показана поверхность эквивалентного тела, соответствующая варианту I ё = 1, ¡с = 2,5 со следующими значениями геометрических параметров: ¡. = 2,7, йт = 0,286 при Re = 5-106. Часть

Х (22) поверхности полубесконечного хвостового участка, имеющая постоянный диаметр йт и расположенная вниз по потоку от граничной линии Ь . , на данном рисунке не отображена.

Рис. 5. Поверхность эквивалентного тела с полубесконечным хвостовым участком ¡^ = 1,5, ¡с = 2,5 при = 5 ■ 106

На рис. 6 и 7 представлены распределения модуля безразмерной скорости V / V. и коэффициента давления ср по поверхности эквивалентного тела для варианта ¡ё = 1,5; ¡с = 2,5 и Re = 5 -106.

В качестве примера применения предложенной методики рассмотрен сравнительный анализ влияния изменения исходных данных из двух наборов вариантов на распределение коэффициента давления по поверхности обтекаемого тела. Сравнение проведено при значении Re = 5 • 106, взятом из срединной части исследованного диапазона чисел Рейнольдса 106...2 -107. Для каждого из наборов построены

зависимости ср = /(х, ¡г) коэффициента давления ср от значения

безразмерной абсциссы х точки поверхности и суммарной длины обтекаемого тела ¡1 (рис. 8, 9).

Рис. 6. Распределение модуля безразмерной скорости V / ¥ю по поверхности эквивалентного тела I = 1,5, I = 2,5 при Яе = 5 • 106

6,935е-001

5,744е-001

4,552е-001

3,361е-001

2,169е-001

9,781е-002

-2,133е-002

-1,405е-001

-2,596е-001

Рис. 7. Распределение коэффициента давления ср по поверхности эквивалентного тела I = 1,5, I = 2,5 при Яе = 5 -106

8 10 12 х

Рис. 8. Зависимость коэффициента дав- Рис. 9. Зависимость коэффициента давления с на поверхности эквивалентно- ления с на поверхности эквивалент-

го тела от безразмерных продольной координаты х и длины обтекаемого тела I при I = 1, Яе = 5' 106:

1 — ¡1 = 4; 2 — ¡1 = 6; 3 — ¡< = !

ного тела от безразмерных продольной координаты х и длины обтекаемого тела ¡( при I с = 6, Яе = 5' 106:

1 — I = 1; 2 — I = 1,5; 3 — I = 2

g g g

4 — ¡1 = 10; 5 — ¡1 = 12

В первом наборе данных длина головной оживальной части сохраняется постоянной ¡ ё = 1, а варьируемая длина цилиндрической части

принимает значения 3, 5, 7, 9, 11, что соответствует суммарной длине обтекаемого тела ¡г, равной 4, 6, 8, 1о и 12, и линиям 1, 2, 3, 4 и 5 (см.

рис. 8). Анализ зависимостей ср = /(х, ¡г) показывает, что при постоянной длине головной части увеличение длины цилиндрической части сопровождается некоторым падением абсолютной величины коэффициента донного давления срй, который достигается на донном срезе

при х = ¡г. Во всем исследованном диапазоне изменения длины цилиндрической части падение модуля коэффициента донного давления носит монотонный характер. Другой характерной особенностью поведения графиков ср = /(х, ¡г) при постоянной длине носового участка

является слабая зависимость от длины цилиндрической части, а следовательно, и от суммарной длины тела значений коэффициента давления на поверхности головной части и прилегающего к ней участка цилиндрической части. При возрастании длины цилиндрической части ¡ с с 5 до 11 существенные отличия значений коэффициента давления

ср начинают проявляться лишь при х > 4. Более того, начиная со значения ¡г = 4 увеличение безразмерной суммарной длины тела ¡г = ¡^ на две единицы до ¡г2 = ¡{1 + 2 приводит к тому, что при х < ¡{1 - 2

значения коэффициентов давления на графиках зависимостей ср = /(х,¡п) и ср = /(х,¡г2) мало отличаются дуг от друга. Различия

проявляются только при х > ¡1 - 2.

Изменение распределения коэффициента давления в носовой части обтекаемого тела становится заметным при рассмотрении зависимостей ср = /(х, ¡г) для второго набора исходных данных, в котором при неизменной длине цилиндрического участка ¡с = 6 варьируется длина оживальной головной части ¡ё, принимающая

значения 1; 1,5 и 2. Указанным вариантам соответствуют значения ¡{ = 7; 7,5; 8 и кривые 1, 2 и 3 (см. рис. 9). Для такого набора геометрических параметров увеличение безразмерной суммарной длины тела ¡г слабо сказывается на распределении давления в донной области и сопровождается лишь незначительным падением абсолютной величины коэффициента донного давления Срй. При обтекании головной части в области, расположенной на некотором удалении от

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

ного абсолютной величине

cp ■

rmm

минимального значения коэффи-

циента давления.

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

ЛИТЕРАТУРА

[1] Белоцерковский С.М., Ништ М.И., Котовский В.Н., Федоров Р.М. Трехмерное отрывное обтекание тел произвольной формы. Москва, ЦАГИ, 2000, 265 с.

[2] Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Cambridge University Press, 2005, 592 p.

[3] Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. Москва, Ин-т механики МГУ им. М.В. Ломоносова, 2006, 184 с.

[4] Щеглов Г. А. Модификация метода вихревых элементов для расчета гидродинамических характеристик гладких тел. Вестник МГТУ им. Н.Э. Баумана. Серия Машиностроение, 2009, № 2, с. 26-36.

[5] Дергачев С.А., Щеглов Г.А. Моделирование обтекания тел методом вихревых элементов с использованием замкнутых вихревых петель. Научный вестник МГТУ ГА, 2016, № 223, с. 19-27.

DOI 10.26467/2079-0619-2016--223-19-27

[6] Плюснин А.В. Моделирование внутреннего и внешнего нестационарного взаимодействия корпуса летательного аппарата с жидкостью методом граничных элементов. Математическое моделирование и численные методы, 2014, № 2, с. 77-100.

[7] Тимофеев В.Н. Математическое моделирование отрывного дозвукового обтекания осесимметричных тел с учетом донного давления. Инженерный журнал: наука и инновации, 2014, вып. 10.

URL: http://engjournal.ru/catalog/mathmodel/aero/1246.html (дата обращения 17.10.2017).

[8] Тимофеев В.Н. Построение полубесконечного эквивалентного тела при математическом моделировании дозвукового отрывного осесимметричного

обтекания. Математическое моделирование и численные методы, 2016, № 4 (12), с. 67-83.

[9] Гогиш Л.В., Степанов Г.Ю. Отрывные и кавитационные течения. Основные свойства и расчет модели. Москва, Наука, 1990, 384 с.

[10] Тимофеев В.Н. Математическое моделирование дозвукового пространственного обтекания тел. Актуальные направления развития прикладной математики в энергетике, энергоэффективности и информационно-коммуникационных технологиях. Труды международной научной конференции, посвященной 180-летию МГТУ им. Н.Э. Баумана. Москва, МГТУ им. Н.Э. Баумана, 2010, с. 198-202.

[11] Тимофеев В.Н., Бушуев А.Ю. Математическое моделирование дозвукового обтекания тел с отрывом потока в донной области. Вестник СГТУ, 2012, № 1 (64), вып. 2, с. 11-14.

[12] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 2. Универсальные законы механики и электродинамики сплошных сред. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 560 с.

[13] Соболев С.Л. Уравнения математической физики. Москва, Наука, 1992, 432 с.

[14] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 1. Тензорный анализ. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 464 с.

[15] Лойцянский Л.Г. Механика жидкости и газа. Москва, Дрофа, 2003, 840 с.

[16] Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент (в математической физике, аэродинамике, теории упругости и дифракции волн). Москва, ТОО «Янус», 1995, 520 с.

[17] Аубакиров Т.О., Белоцерковский С.М., Желанников А.И., Ништ М.И. Нелинейная теория крыла и ее приложения. Алматы, Гылым, 1997, 448 с.

Статья поступила в редакцию 27.09.2017

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

Тимофеев В.Н. Особенности вихревой схемы при моделировании дозвукового отрывного обтекания с полубесконечным эквивалентным телом. Математическое моделирование и численные методы, 2017, № 4, с. 73-91.

Тимофеев Валерий Николаевич — канд. техн. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор более 40 научных и методических работ. Область научных интересов: математическое моделирование, численные методы, механика жидкости и газа, аэродинамика. e-mail: [email protected]

Special features of vortex diagram in simulation of subsonic detached flow around the semi-infinite equivalent body

© V.N. Timofeev Bauman Moscow State Technical University, Moscow, 105005, Russia

The paper introduces some special features of mathematical simulation of subsonic detached flow around the bodies, the flow being localized in the vicinity of the ground shear. The formation of vortex diagram for the semi-infinite equivalent body is examined. The formulae for determining the vector functions of the vortex segments velocity are re-

duced to a form allowing one to easily pass to the limit as the points of the origin or ends of these segments are moved off to infinity. Furthermore, the study shows the relationships for finding the velocity function vectors of semi-infinite vortex segments and U-shaped vortex lines, the relationships being adapted for computer calculations. Findings of mathematical simulation of the flow around cylindrical bodies with the head part of the ogival form are given.

Keywords: mathematical simulation, subsonic detached flow, equivalent body, vortex diagram, discrete vortex method, base pressure

REFERENCES

[1] Belotserkovskiy S.M., Nisht M.I., Kotovskiy V.N., Fedorov R.M. Trekhmernoe otryvnoe obtekanie tel proizvolnoy formy [Three-dimensional detached flow of the bodies of arbitrary form]. Moscow, TsAGI (Central Aerohydrodynamic Institute) Publ., 2000, 265 p.

[2] Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Cambridge University Press, 2005, 592 p.

[3] Andronov P.R., Guvernyuk S.V., Dynnikova G.Ya. Vikhrevye metody rascheta nestatsionarnykh gidrodinamicheskikh nagruzok [Vortex methods of calculation of nonstationary hydrodynamic loads]. Moscow, Institute of Mechanics Lomono-sov MSU Publ., 2006, 184 p.

[4] Scheglov G.A. Vestnik MGTU im. N.E. Baumana. Ser. Mashinostroenie — Herald of the Bauman Moscow State Technical University. Series Mechanical Engineering, 2009, no. 2, pp. 26-36.

[5] Dergachev S.A., Scheglov G.A. Nauchnyy vestnik MGTU GA — Bulletin of Civil Aviation High Technologies, 2016, no. 223, pp. 19-27.

DOI 10.26467/2079-0619-2016--223-19-27

[6] Plyusnin A.V. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2014, no. 2, pp. 77-100.

[7] Timofeev V.N. Inzhenernyi zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2014, no. 10. Available at: http://engjournal.ru/ catalog/mathmodel/aero/1246.html (accessed October 17, 2017).

[8] Timofeev V.N. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2016, no. 4 (12), pp. 67-83.

[9] Gogish L.V., Stepanov G.Yu. Otryvnye i kavitatsionnye techeniya. Osnovnye svoiystva i raschet modeli [Detached and cavitational flows. Basic properties and model calculation]. Moscow, Nauka Publ., 1990, 384 p.

[10] Timofeev V.N. Matematicheskoe modelirovanie dozvukovogo prostranstven-nogo obtekaniya tel [Mathematical simulation of subsonic spatial flow around bodies]. Aktualnye napravleniya razvitiya prikladnoy matematiki v energetike, energoeffektivnosti i informatsionno-kommunikatsionnykh tekhnologiyakh. Trudy mezhdunarodnoy nauchnoy konferentsii, posvyashchennoy 180-letiyu MGTU im. N.E. Baumana [Current trends in the development of applied mathematics in energy, energy efficiency and information and communication technologies. Proceedings of the international scientific conference dedicated to the 180th anniversary of BMSTU]. Moscow, BMSTU Publ., 2010, pp. 198-202.

[11] Timofeev V.N., Bushuev A.Yu. Vestnik SGTU — Vestnik Saratov State Technical University, 2012, no. 1 (64), iss. 2, pp. 11-14.

[12] Dimitrienko Yu.I. Mekhanika sploshnoy sredy. V 4 tomakh. Tom. 2. Universal-nye zakony mekhaniki i elektrodinamiki sploshnykh sred [Continuum mechanics. In 4 vols. Vol. 2. Universal laws of continuum mechanics and electrodynamics]. Moscow, BMSTU Publ., 2011, 560 p.

[13] Sobolev S.L. Uravneniya matematicheskoy fiziki [Equations of mathematical physics]. Moscow, Nauka Publ., 1992, 432 p.

[14] Dimitrienko Yu.I. Mekhanika sploshnoy sredy. V 4 tomakh. Tom 1. Tenzornyy analiz [Continuum mechanics. In 4 vols. Vol. 1. Tensor analysis]. Moscow, BMSTU Publ., 2011, 464 p.

[15] Loytsyanskiy L.G. Mekhanika zhidkosti i gaza [Fluid mechanics]. Moscow, Drofa Publ., 2003, 840 p.

[16] Lifanov I.K. Metod singulyarnykh integralnykh uravneniy i chislennyy eksperi-ment (v matematicheskoi fizike, aerodinamike, teorii uprugosti i difraktsii voln) [The method of singular integral equations and numerical experiment (in mathematical physics, aerodynamics, theory of elasticity and diffraction of waves)]. Moscow, LLP Yanus Publ., 1995, 520 p.

[17] Aubakirov T.O., Belotserkovskiy S.M., Zhelannikov A.I., Nisht M.I. Nelinei-naya teoriya kryla i ee prilozheniya [Nonlinear wing theory and its applications]. Almaty, Gylym Publ., 1997, 448 p.

Timofeev V.N., Cand. Sc. (Eng.), Assoc. Professor, Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. Author of over 40 scientific and methodical works. Science research interests: mathematical simulation, numerical methods, fluid mechanics, aerodynamics. e-mail: [email protected]