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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XVII 1986

№ 5

УДК 532.527

ПРОЕКЦИОННЫЙ МЕТОД РАСЧЕТА ХАРАКТЕРИСТИК ОТРЫВНОГО ОБТЕКАНИЯ ТЕЛ ИДЕАЛЬНОЙ ЖИДКОСТЬЮ

А. В. Воеводин, Г. Г. Судаков

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

В настоящее время для расчета отрывного обтекания крыльев малого удлинения широко используется теория удлиненных тел. В этом приближении исходная трехмерная задача сводится к плоской нестационарной задаче об отрывном обтекании деформирующегося контура, представляющего собой поперечное сечение исходного тела. При этом роль времени играет продольная координата х. Для расчета возникающих двумерных задач имеется два основных метода: итерационный метод (например, [1 —3]) и метод установления (например, [4—6]). Существенным моментом обоих подходов является использование аналитического выражения для конформного отображения поперечного сечения тела на стандартную область (например, круг). Это, с одной стороны, позволяет точно удовлетворить условию непротекания на поверхности крыла и фюзеляжа, а также аналитически сформулировать условие Жуковского на острой кромке крыла. С другой стороны, необходимость иметь аналитическое выражение для конформного отображения приводит к существенному сужению класса исследуемых тел. Решенные с помощью указанных^подходов задачи в основой касаются либо крыльев нулевой толщины без поперечной кривизны, либо комбинаций таких крыльев с фюзеляжем простой формы. Имеется только один пример расчета отрывного обтекания толстого крыла, а именно: в работе [7] для расчета обтекания крыла ромбовидного поперечного сечения использовалась простейшая модель вихревой

пелены (вихрь — разрез). Кроме того, в работе [8] приведены результаты расчета отрывного обтекания тела полукруглого поперечного сечения на основе более точной модели вихревой пелены.

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

1. Постановка задачи. Пусть их — скорость набегающего потока, а ^ 1 •—угол атаки, t = x/uoa. Теория удлиненных тел позволяет в первом приближении свести исследование отрывного обтекания тел малого удлинения к отысканию потенциала возмущенного течения f (z, у, t) в плоскости поперечных переменных Z, у. Функция ср должна удовлетворять двумерному уравнению Лапласа, граничному условию д<?/ду = а при (z2 +у2)112 -* оо, условию не-протекания на поверхности тела, условию отсутствия скачка давления и нормальной составляющей скорости при переходе через пелену, а также условию Жуковского при наличии острых кромок. Введем комплексный потенциал скорости W (о t), где a = z+iy, и Real(W) = <p. Поле скоростей можно представить в виде суммы скорости набегающего потока и скоростей, индуцируемых особенностями, распределенными на контуре тела Ly (источники или вихри) и линии тангенциального разрыва скорости Z,2 (вихри):

где С2 — действительный параметр (интенсивность особенностей). Из условия непротекания следует

где Г (г, у, £) — 0 — уравнение контура тела п— пг-\- 1пу, пг и пу — компоненты вектора единичной внутренней нормали к контуру тела в поперечной плоскости гу.

Эволюция вихревой пелены 12 в плоскости комплексной переменной а определяется уравнением

2. Описание численного метода. Для решения задачи (1)—(3) сформулируем процедуру, основанную на проекционном методе [9].

dw

д о

— —I (х ■

(1)

L,

L. + L,

(2)

(3)

Введем разбиение контуров I, и Ь2 на N отрезков ДА — [а (С^г-и и определим на этих отрезках конечно-элементный базис 1:

[ 1» 3 (:

!**■«?)= { -0, о£Д*, '*=1, ЛГ- 1,

[ 8 [о — а ((2дг)], о^Длг,

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

(Я, 0 = .МОМС).

тогда

где

N-1

а (С?, Ь) = [ак (<3 — (?*_!) + [Д.* ((?) , (4)

й= 1

*=1...... ^-1.

У* — V*—1

Из формулы (4) следует, что часть пелены вне ядра, а также крыло и фюзеляж аппроксимируются прямолинейными отрезками (панелями). Подставляя (4) в выражения (1) для комплексно-сопряженной скорости, получим:

Л?1

дИУ ' —

д а

<5>

+ V Г | 1 в*-влг-1

2тсг Л а* (о — о') 2яг а—а«г

й= ЛГ!+1 Ай я ' Л'

где Л^1 — количество отрезков Ьи на которых введено распределение источников.

Далее введем второе разбиение контуров Ьх и Ьг на ./V* отрезков Дт = [о (^т_х), о(<Зт)] и базис 2, вообще говоря, не совпадающий с первым:

‘ 1 гТ

I г __г > 3 С ^т >

| ^т—1

*т («) =

I.

О, а£Ьт, т—\, ... , /V* — 1,

8 [з - о (С^)], о £ АЛ',

где Сл = в(фт). Пусть Дт = Д* при /и = & = 1, ... ,7^ и т-Л^ —А2 = = /г ——Л/2= 1, ... , Л^з, где Д^2 — количество отрезков/-!, на кото-

рых введено распределение вихрей, Л'! + N2 = Р — количество отрезков второго разбиения Ьи N3 — количество отрезков на Проектируя уравнение (2) на базис 2, получим

* ъ(С) йс

д^\2 11/2

+ *т(С)<К = 0,

что эквивалентно системе уравнений

С) ^с=о,

О а

(6)

т=\, ..., Р.

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

Производя аналогичную операцию с уравнением (3) для вихревой пелены, получим систему уравнений

Формула (7) означает, что точка середины панели Дт на вихревой пелене движется со скоростью Фш, равной значению среднеинтегральной скорости по панели Ат.

Оставляя в левой части уравнений (6) лишь слагаемые, зависящие от интенсивности особенностей, распределенных на контуре тела, получим систему линейных уравнений для определения Д = — С}ь-1 на при заданной вихревой пелене. Уравнения (7)

позволяют рассчитать перемещение точек вихревой пелены при известном распределении особенностей на контуре тела. Если вихревая пелена сходит с острой кромки крыла, необходимо потребовать выполнения условия Жуковского. Учитывая (4), при произвольном задании плотностей распределения особенностей на прилежащих к кромке панелях и возникает особенность в скорости. Поэтому плотность завихренности на прилежащей к кромке панели Ь2 должна определяться из условия отсутствия этой особенности (см. п. 4).

Система уравнений (6), (7) совместно с условием в точке отрыва позволяет рассчитать течение при при заданной в момент ^

вихревой пелене. Вычислительная процедура состоит из двух этапов. На первом этапе из условия непротекания (6) при заданной вихревой пелене находится распределение особенностей по поверх-

(7)

ности тела. Одновременно из условия в точке отрыва определяется циркуляция участка пелены, сошедшего с тела на предыдущем шаге по времени. Затем решаются уравнения (7) и определяется положение точек в следующий момент времени. Для интегрирования (7) используется метод Эйлера второго порядка. Исключение составляет участок вихревой пелены, сошедший с тела за время Д£, так как определение величины Ср+1(£ + Д^) не сводится к решению уравнения (7). Эта величина определяется по формуле

где Сотр — координаты точки отрыва, М—номер отрезка контура тела для которого См = С0Тр. вещественные коэффициенты сх и с, находятся из условия

Две последние формулы получены из условия, что конфигурация вихревой пелены в окрестности точки отрыва асимптотически описывается уравнением у = кх312 (см., например, [5]).

Восстановление носителей базиса 2 (отрезков Дт) на 12 проводится по формулам

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

3. Исследование устойчивости численного метода. Рассмотрим уравнение (3) в случае уединенной вихревой пелены. Тогда (3) можно переписать в виде

Ср+2 = С0Тр + ''-И (С-, + ^С‘>) •

где (2 — текущая циркуляция вихревой пелены Ь2. Рассмотрим далее малое возмущение вихревой пелены

а = а0 + 8а ,

где |3о[<|о0|. Тогда из (8) следует, что в линейном приближении

Уравнение (9) является основным при исследовании устойчивости метода. Будем искать решение уравнения (9) в виде гармонической волны длиной Х = 2тс/<в, малой по сравнению с характерным размером вихревой пелены. Тогда приходим к задаче об эволюции гармонического возмущения на плоском бесконечном тангенциальном разрыве скорости (задача 1), а именно:

о0 = г , 8о = е (^) еЫг . (10)

Подставляя (10) в (9) и вычисляя интеграл, получим уравнение

(11)

которое является математическим описанием неустойчивости Гельм-гольца. Здесь и далее ^0 = —— плотность циркуляции невозмущенного тангенциального разрыва. Уравнение (11) приводит к экспоненциальному (по времени) росту возмущений:

I /л I “То ,

£ ( ) 1 = е2 , (12)

I • (0) I ' '

причем амплитуда волн длины Х-^0 растет бесконечно быстро. Это означает, что задача Коши (8) некорректна по Адамару. Для решения некорректных задач (8) разработана специальная техника (см. [10]), которая существенно базируется на методе дискретных вихрей и требует введения регуляризатора. Роль регуляризатора сводится к подавлению роста возмущений с короткими длинами волн.

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

Разобьем плоский бесконечный тангенциальный разрыв на систему одинаковых отрезков ДА с циркуляцией Д <3Й (6 = 0, +1, ...). Тогда процедура п. 2, примененная к (8), дает

<13,

<Ц 2пг А (5„ ^__со ап дй ъ С

* 1

где Д6 —прямолинейный отрезок [С*, С*-ц]> ^к = — (С* + С*+1> .

Положим

С * = АА + е (*)«*» <*-!)*, |е|«1,

тогда, вычисляя интегралы в правой части (13) и оставляя в (13) лишь линейные по е члены, получим

Уравнение (14) является аналогом уравнения (13), при этом в (14) используется лишь процедура аппроксимации скорости, индуцированной вихревой пеленой (аппроксимация по пространству, но не по времени). Характерно, что метод дискретных вихрей приводит также к уравнению (14), но с другой функцией /{ч>Н) =

= ^ . графики зависимости приведены на рис. 1,

где цифра 1 соответствует /= 1 (точное решение), цифра 2 — /—

= 1—— (метод дискретных вихрей), 3 — проекционный метод

[формула (14)]. Видно, что проекционный метод дает значительно лучшую аппроксимацию точного решения во всем диапазоне длин волн.

Для численного решения уравнения (14) процедура п. 2 предлагает следующую схему. Пусть используется численный метод интегрирования второго порядка. Тогда

е1 (^ + Д і) = е (і)— іАгЦ)Аі, )

®2 + А і) — є (і) — іА — [є (і) + є1 {і Д і) ] М. I

Наконец, процедура построения нового базиса по возмущенным значениям узлов приводит к окончательному результату

4 — точное решение [формула (14)], 5 — метод дискретных вихрей [5], 6 — данный метод [формула (17)]. Как легко видеть, данный метод подавляет возмущения с длинами волн и значительно

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

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

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

Таким образом, базис 2 содержит на один элемент больше, чем базис /. Это обстоятельство позволяет определять плотность циркуляции прилегающего к кромке участка вихревой пелены из решения системы уравнений (6) и автоматически удовлетворить условию Жуковского на острых кромках. На рис. 2 и 3 приведены

£ и + Д0 = [ £ {*) + -у [е (і) + е1 (£ + Д/)] д/1 соэ2^ =

| е(<) + іАІ(0Д<+ у М*) Д**} . (16)

Отделяя в (16) действительную и мнимую части, получим

Іт є (і + М) 1га £ (і)

(17)

На рис. 1 приведены значения р+=р+

о

для —— = 2 : к

Сл?,+1 = 3^ + 1 ,

£р + 2 = аР+ 2 •

результаты расчета формы вихревой пелены, образующейся при обтекании тонкого треугольного крыла при а/5 = 1; 1,5; 2, и комбинации тонкого треугольного крыла с фюзеляжем в виде кругового конуса при а/5 = 1, 2 и 5*/5 = 0,625, где а — угол атаки, 8 — полуугол при вершине крыла, ,8* — полуугол при вершине конуса. Результаты расчетов хорошо согласуются с результатами расчетов, полученными итерационными методами [1, 3] (сплошная линия).

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

( ^м+1 ~ | ^еа' ('*N ~*м) I ^ |

а^==ам+1-----------------:—...— ------- ,

I *ЛГ-1| Кеа1 (^Ж + 1 \и) I '1Я +1 I

Сл1 — '^отр ■

Форма вихревой пелены для крыльев ромбовидного поперечного сечения с относительной толщиной х ==0,1; 0,2 и 0,4 при а/5 = 1 показана на рис. 4. На рис. 5 показана рассчитанная форма вихревой пелены для комбинации толстое крыло — корпус, изображенной на том же рисунке, при а/8 = 1; 1,5; 2.

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

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

Im J п д£ v/W+1 (С) rfC = 0 , С.и = £„тр ,

Ам+\

которое решается совместно с системой (6).

На рис. 6 приведено сравнение результатов расчетов формы вихревой пелены для угла отрыва 60 == 132° и а/5* = 2,52 (сплошная линия) с данными работы [6] (штрихпунктирная линия) и работы [11] (пунктирная линия), в которой 60г=131 и а/8* = 2,5.

ЛИТЕРАТУРА

1. Smith J. Н. В. Improved calculations of leading-edge separation from slender delta wings. Proc. Roy. Soc., 1968, A 306.

2. P u 1 1 i n D. J. A method for calculating inviscid separation from about conical slender bodies.—ARC RM, 1973, N 140.

3. Воеводин А. В. Исследование неединственности решения задачи об отрывном обтекании системы крыло — фюзеляж малого удлинения.— Ученые записки ЦАГИ, 1979, т. X, № 1.

4. Sacks А. М., Lundberg R. Е„ Hanson С. W. A theoretical investigation of aerodynamics of slender wing-body combinations axhibiting leading-edge separation.— NASA CR, 1967, 719.

5. Судаков Г. Г. Расчет отрывного течения около тонкого треугольного крыла малого удлинения.— Ученые записки ЦАГИ,

1974, т. V, № 2.

6. Захаров С. Б. Расчет невязкого отрывного обтекания тонкого кругового конуса на больших углах атаки.— Ученые записки ЦАГИ, 1976, т. VII, № 6.

7. С 1 а г k R. W., S m j t h J. H. В., T h о m p s о n С. W. Some series-expansion solutions for slender wings wilh leading — edge separation.—

ARC RM, 1975, N 3785.

8. Воеводин А. В. Отрывное обтекание конической комбинации крыла малого удлинения с несимметричным фюзеляжем.—Ученые записки ЦАГИ, 1983, т. XIV, № 4.

9. М а р ч у к Г. И., А г о ш к о в В. И. Введение в проекционносеточные методы.— М.: Наука, 1981.

10. Молчанов В. Ф. О реализации метода плоских сечений в нелинейной теории крыла.— Ученые записки ЦАГИ, т. V, Л6 2,

1974.

11. Smith J. Н. В. Achievements and problems in modeling highly— swept flow separations. Numerical methods in aeronautical fluid dynamics.— Academic Press, 1982.

Рукопись поступила 30/IV 1985 г.