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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Том XVIII 1987 №5

УДК 533.6.013.2.011.32

РАСЧЕТ БЕСЦИРКУЛЯЦИОННОГО ОБТЕКАНИЯ ПРОИЗВОЛЬНЫХ ТЕЛ

С. М. Белоцерковский, И. К. Лифанов, А. А. Михайлов

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

Метод дискретных вихрей развит для нахождения характеристик бесциркуляционного обтекания идеальной несжимаемой жидкостью произвольных тел как гладких, так и имеющих углы. К изучению такого рода схем приводит упрощенный учет некоторых нестационарных воздействий жидкости на тело, когда вихревой след не оказывает заметного влияния на его обтекание. Тогда получается схема, в которой нет вихревого следа, а нестационарное воздействие описывается с помощью присоединенных масс [1, 2]. В результате сильно упрощаются не только задачи аэродинамики, но и динамики движения и аэроупругости.

Эффективный численный метод решения бесциркуляционных задач для гладких тел произвольной конфигурации развит в работе [3] на основе потенциала простого слоя. Здесь скорости и давления конечны во всем пространстве, и вопрос о выделении особенностей не возникает. Практический и принципиальный интерес представляет расчет бесциркуляционного обтекания тел произвольной формы с угловыми точками, в которых скорости и давление обращаются в бесконечность. Данная; работа обобщает метод дискретных вихрей на широкий класс бесциркуляционных задач, применимый ко всем плоским и пространственным конфигурациям (тонкие крылья, тела конечной толщины как гладкие,, так и с изломами).

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

тодом. Кроме того, упрощаются формулы для потенциала и для скачка потенциала на поверхности тел. ----

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

1. Общие положения. Течение везде вне тела считается потенциальным (вихревой след за телом не образуется). Задача сводится к нахождению потенциала возмущенных скоростей Ф, удовлетворяющего уравнению Лапласа вне тела и условию непротекания на поверхности тела, т. е. к решению задачи Неймана:

ДФ = 0 , М £(£>-),

дФ

а«ж„

-/(Л10), М0£а.

(1)

Здесь ф_)—дополнение к области ф), занятой телом, пмо—направление внешней нормали к поверхности а тела (£)) в точке М0, а функция

/ определяется видом движения.

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

= л. (2)

Здесь — возмущенная, V?* — переносная скорости; п — единичная

нормаль к поверхности тела; г — радиус-вектор точки, в которой выполняется условие непротекания.

Потенциал скоростей можно представить в виде линейной комбинации параметров движения

6

Ф = ^о X ~~Ь ^0 у ?2 + Цр г Уз + Яс ?4 + &2 Уб = у1. ^1 •

/=1

Тогда в выражения для сил и моментов, действующих на тело, войдут коэффициенты Хгз (так называемые присоединенные массы), имеющие вид

\и = — р

ч v ” Ti дп

где р — плотность среды, а интегрирование ведется по поверхности тела а [1].

Вычисление коэффициентов %ij сводится к нахождению потенциалов dlf:

ф,-, a -q— определяется граничными условиями (2) при соответствующем виде движения:

У = 1, .... 6;

дп ‘

fi = cos nx ; /2 =: cos пу \ /8 = cos nz ; j\ = у cos nz — z cos ny ; fb — z cos nx. — x cos nz ;

/6 =» x cos ny — у cos nx .

2. Сингулярные интегральные уравнения. Решение задачи (1) будем искать в виде потенциала двойного слоя (традиционно его искали в виде потенциала простого слоя):

Ф =

А

и =

дп

м

1п

4тс

1

-1 -2,

— —2

где 7 — плоскопараллельное течение; 2—трехмерное; а — кусочно-ляпу-новский контур в плоском случае (поверхность — в пространственном).

Функция Ф(-Мо) гармонична в (I)-) и в (В+) и терпит разрыв на границе [5], равный g(Mo)

Ф+ (Л*„) = -у ё (М0) + Ф (М0) ;

ф-№=-уг(л0)+ф(іИ0).

(3)

Здесь Ф+(М0), Ф_(Л1о)—значения потенциала на внешней и внутренней стороне границы а.

Плотность g(M0) будем находить из условия

<5Ф

дп

м.

/Шо).

Нормальная производная потенциала непрерывна на границе а [7], поэтому в равенстве

дФ

дпм„ А

Г ё(М) — йо

А дпъ) дп~м

можно перейти к пределу при М0, стремящейся по нормали Пм„ к границе сг. В результате получим уравнение

(4)

Если плотность двойного слоя удовлетворяет уравнению (4), то функция Ф (М) будет решением исходной задачи Неймана с точностью до константы. Когда сг является замкнутой поверхностью, то согласно [6]

5/(к)*°м=о-

(5)

В равенстве (4) знак нормальной производной можно внести

М0

под интеграл в силу непрерывности нормальной производной потенциала двойного слоя.

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

-тГ в>)Л=ЯЛ,о), М0£а. (6)

А .3 дпма V дпм /

9

Интеграл в (6) надо понимать как интеграл от обобщенных функций [6].

Для уравнения (6) константа является собственной функцией, так

как

Г и ) с1ам = 0 .

3 дпм0\дпм )

<3

Поэтому решение этого уравнения определено с точностью до константы и, следовательно, единственное решение можно выделить заданием значения g(M) в некоторой точке или значения интеграла от функции g(M) по сг.

Задачу (1) и уравнение (6) можно рассматривать и для незамкнутых поверхностей а. При этом для функции /(Л40) в общем случае уже не будет выполняться равенство (5).

Для плоскопараллельного течения уравнение (6) имеет вид

±5 г(т '(Ъ'О-ТО Г'Ъ) а^лм,), *,€.. (7,

а Г Если а имеет углы, то М0 не должна совпадать с ними.

В пространственном случае уравнение (6) можно преобразовать к виду

> |С НМ ) м (8)

4 тс Г5

а

В частном случае, когда граница а —квадрат [—1,1]Х{—1,1] на плоскости Охг, будем иметь

г = (х — х0)7 + (г — г0)~] , пм — пм0=], й? = йх<1г и уравнение (8) приобретает вид

— Г Г £■(•*> 2)---------<1хЛг_----- = / ^ ^

4* ] ] * [(х-х0У + {.г-га??'2

—1—1

3. Численное решение задачи методом дискретных вихрей. Потенциал скоростей, индуцированный замкнутой вихревой нитью постоянной интенсивности, равен потенциалу двойного слоя поверхности, натянутой на нить, с плотностью, равной интенсивности вихревой нити [2]. Условие непротекания от вихревой нити интенсивности Г будет иметь вид (6), в котором я(Л1) надо заменить на Г.

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

ность разобьем на р штук дискретных указанных вихревых образований интенсивности Гй, £=1,...,р. Условие непротекания будем выполнять в расчетных точках М0ь, взятых в центре 6-го вихревого образования (рис. 1). Тогда задача нахождения неизвестных интенсивностей Гл

Рис. 1

сведется к решению следующей системы линейных алгебраических уравнений:

р

'£Тк™Ьт==/т> /Л-1,... ,Р, (9)

й = 1

где ™кт = {Укт, п0т); Укт - скорость, индуцируемая, в точке М0т &-й вихревой рамкой единичной интенсивности;/т = /(/И0т). Систему (9) можно получить из уравнения (6) следующим образом. Поверхность о разбивается на части ой, /г= 1, .... , р, ограниченные &-м вихревым образованием. На функция g (М0) заменяется константой £(М0*), и уравнение (6) переходит в систему

1]'<Мт1г4-(^“)Л*Н<л,'-ь (10)

к = 1 а к °т

Если в этой системе обозначить

£ (ЛГ0 к) = Г* ; -у С -4г~ и ) = п)кт ,

А дпмот \дпм /

”и

то (10) примет вид (9).

При моделировании обтекания тонкой несущей поверхности, когда <т в уравнении (6) незамкнутая кривая (например, отрезок) или незамкнутая поверхность (например, часть плоскости), то доказано, что система (9) невырождена. В случае отрезка также доказана и сходимость численного решения к точному [7].

Если а в уравнении (6)—замкнутая кусочно-гладкая линия (телесный профиль) или замкнутая кусочно-гладкая поверхность тела конечной толщины, то система (10) линейно зависима. Действительно, сумма элементов в т-й строке (т = '1,...,р), равна нулю (она равна зна-

чению правой части уравнения (6) при £(.М) = 1) и, таким образом, столбцы матрицы этой системы линейно зависимы. Так как в расчетах для Хц значения Гь, к=\,...,р нужны с точностью до константы, то можно считать, что Гд равна нулю, и (9) приобретает вид:

Г*®*т=/т > т— 1, .... р . (11)

к = 1 к + я

Система (11) переопределена (число неизвестных меньше числа уравнений), поэтому для получения невырожденной определенной системы вводится регуляризирующая переменная уор, которая при л->оо стремится к нулю

р

Тор+ тк™кт=/т, т=\,...,р. (12)

к = 1

к+ч

В случае, когда а — замкнутая кусочно-гладкая кривая, ниже доказывается сходимость решения системы (12) к точному.

Приближенное значение потенциала на внешней стороне поверхности а в точке М0т находится по формуле

ф(^от)=~У г*^—- Аа*~Т , (13>

л гкт 1

а присоединенные массы

Х„=- р V Ф‘ (М0т) Дат . (14)

ОПЛ1о т

Формула (13) получается из (3) заменой интеграла квадратурной суммой, причем, через <]/кт обозначен угол между внутренней нормалью в точке Мп и радиус-вектором # = К (М0т, Мк), а через Аза — площадь элемента поверхности оА. Индекс г или у над Ф соответствует тому, что в системе (9), (12) справа берется функция, соответствующая г'-му (или у-му) движению.

В случае расчета присоединенных масс тонких поверхностей методика расчета присоединенных масс упрощается, так как вместо (14) в этом случае будет

ч=-р 2АФ‘(уИ0-)^-д5-

т = 1 М0 щ

где АБт — часть тонкой поверхности 5, ограниченная т-й вихревой рамкой. Скачок потенциала ДФ*(М0т) легко определяется из (13)

ДФ ‘ (М0 т) = Т‘т . (15>

Используя равенство (15), выражение (14) можно записать в виде

р дФ>

= -р2г'^45.. (16)

, дпм

т=1 7И" ">

4. Обоснование численного метода для профиля. Пусть параметрическое уравнение контура о имеет вид

* = *(е). о<е<2тг,

У-У(в),

причем, х' (б) и у' (0) принадлежат классу Гельдера Н{а) на [0,2тг] и х 2 -\-у'2 ф 0 на [0,2тс], л: (0) и у {В) периодические с периодом 2я. Тогда уравнение (7) можно записать в виде

2тс

g (0) + с кг (0О, 6) g (0) £*0 =/(0о), е0 £ [0>] , (17)

в1п2

где функция (0о, 0) принадлежит классу Я (а) на [0,2 я], если х"(%) и г/ (0) принадлежат этому классу. Первый интеграл слева в уравнении (17) надо понимать как определенный интеграл от обобщенной функции [6].

Для численного решения уравнения (17) надо рассматривать следующие системы линейных алгебраических уравнений.

Возьмем на отрезке [0,2 я] две системы точек

{в* = (Л— 1) А} и { 0Ой = 0л +А|, &= 1, ... , п, А = -^ .

Если известно g(q) =0 в точке qeff, то

Тол + А. ёп (ео*) ^ +

п 2 *= 1 £„(бо*) [*в

Ьфя

п X £».(0о*) Лкт=/(В0т),

й=1 **9

где (0О *) —приближенное значение функции g(д) в точке 0ОЙ при разбиении а на п частей. Параметризация контура выбрана так, что точка является образом одной из точек 0От при т^=-кч> т. е. Хд = X (00 кд), Уд = У{%кд)> причем, Акт имеет один из трех видов

1) К (0От, 0О*)- — ;

2)| М60 т, 0)^0;

2л

в зависимости от следующих причин:

1) если ядро к\ (0о, 0) удовлетворяет условию Гельдера и не содержит в себе дроби, дающей неопределенность 0/0 на диагонали, т. е. при 0о=О;

2) если для к\ (0от, 0) легко находится первообразная;

3) когда (0о, 0) содержит в себе указанную в (1) дробь.

Если известно значение с интеграла от ^ОЭ) на отрезке [0,2 л], то

00 щ — 0ь

-с£ 0т2-- +

П

+ Х ёп^ок) ^я = /(еои) - т= 1, ...

п ,

к= 1

(18)

Во всех этих случаях справедлива следующая теорема.

Если функции /(0) и к (во, 0) удовлетворяют условию Гельдера степени а, то между точным решением уравнения (17) в точках 0От и решением системы (18) выполняется соотношение:

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

Если контур о имеет углы то точки 0ь, £= 1,..., п выби-

раем так, чтобы углы были в числе этих точек (этого можно добиться выбором способа параметризации). В остальном вся численная схема остается такой же. В этом случае ядро &1(0о, 0) в уравнении (17) будет иметь разрывы первого рода, и поэтому сходимость приближенного решения к точному будет той же вне окрестностей углов и интегральной по всему контуру. Порядок сходимости в интегральной метрике будет

5. Примеры расчетов. Вихревая схема для расчета присоединенных масс цилиндрических тел бесконечного удлинения прямоугольного сечения представлены на рис. .1,0. Вместо замкнутых четырехугольных вихревых рамок в плоском случае используются их аналоги — пары вихрей с равнопротивоположными циркуляциями Гь и —Гй.

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

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

Схема расположения вихревых рамок и расчетных точек для вычисления присоединенных масс параллелепипеда квадратного сечения показана на рис. 1,6. Общее количество вихревых рамок на поверхности параллелепипеда было 200—250.

(19)

Рис. 2

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

рЛЧ

10^4 р А» /з

^22 -----

55

Л22 р Ь?1

Ю-А»

^33 —

р А*/з

рАЧ

где /— длина параллелепипеда, К — его высота. Результаты расчетов представлены на рис. 3.

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

Введем безразмерный коэффициент присоединенной массы

____ 4^-22

[1 =

яр Vі I

Присоединенная масса Я22 рассчитывалась по формуле (16). Результаты расчетов для прямоугольных крыльев различных удлинений 1/Ь представлены на рис. 4 в виде зависимости коэффициента (х от 1/Ь.

Рис. 4

Здесь же нанесены экспериментальные данные, заимствованные из работы [8].

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

ЛИТЕРАТУРА

1. Ламб Г. Гидродинамика. — М.; ОГИЗ, 1947.

2. Седов Л. Г. Механика сплошной среды. — М.: Наука, т. II 1976.

3. Маслов Л. А. Произвольное движение продолговатого тела в идеальной жидкости. — Изв. АН СССР, МЖГ, № 966, № 6.

4. Белоцерковский С. М., Лифанов И. К., Михайлов А. А. Моделирование на ЭВМ отрывного обтекания профилей с угловыми точками.—ДАН СССР, 1985, т. 285, № 6.

5. Тихонов А. Н., Самарский А. А. Уравнения математической физики. — М.: Наука, 1966.

6. Курант Р. Дифференциальные уравнения с частными производными.— М.: Наука, 1965.

7. Белоцерковский С. М., Лифанов И. К. Численные методы в сингулярных уравнениях.—М.: Наука, 1985.

8. Р и м а н И. С., К р е п с Р. К. Присоединенные массы тел различной формы. — Труды ЦАГИ, 1947, № 635.

Рукопись поступила 6/\Ш 1986 г.