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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м IX 197 8

№ 2

УДК 629.7.015.4:533.6.013.422:629.7.023.2

ПАНЕЛЬНЫЙ МЕТОД РАСЧЕТА НАГРУЗОК, ДЕЙСТВУЮЩИХ НА КРЫЛО, СОВЕРШАЮЩЕЕ ГАРМОНИЧЕСКИЕ КОЛЕБАНИЯ В ДОЗВУКОВОМ ПОТОКЕ

П. М. Гостев, А. С. Кутин, В. В. Мозжилкин,

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

Наиболее распространенным методом определения нестационарных характеристик летательных аппаратов при дозвуковых скоростях является метод С. М. Белоцерковского [1]. Однако, как отмечено в работах [2, 3], он требует внимательного выбора расчетной схемы и сложен в численной реализации при конечных числах Струхаля. Поэтому в последнее время возрос интерес к численным методам, основанным на непосредственном решении интегральных уравнений линейной теории крыла, в частности» к численным методам решения интегрального уравнения для потенциала скорости [4—10]. Видимо, это связано с тем, что потенциал скорости является кусочно-непрерывно-дифференцируемой функцией, тогда как скачок коэффициента давления Дср имеет разрывы I рода на передней кромке и шарнире руля [11]. Поэтому численные методы для потенциала скорости отличаются большей точностью и лучшей сходимостью по сравнению с методами для Дср.

1. Математическая постановка линейной задачи о вычислении скачка потенциала скорости на плоской несущей поверхности, совершающей гармонические колебания в дозвуковом потоке со скоростью О на бесконечности, формулируется следующим образом [9].

Введем систему координат (х, у, г), связанную с крылом таким образом, что х параллельно потоку, V лежит в плоскости крыла,

% перпендикулярно плоскости хОу. Центр системы координат находится в вершине передней кромки крыла. В безразмерных переменных

где <р — скачок потенциала скорости, р — частота гармонических колебаний крыла, 5 — полуразмах крыла, Ф удовлетворяет интегральному уравнению

Здесь 5'— поверхность крыла, ^ — след, а ЧУ—безразмерный ■скос потока, связанный с размерной составляющей скорости по оси г на крыле т(х, у, 0) соотношением

В силу симметрии крыла Ф (X, У, 0) = ®^, —У, 0). Тогда интегральное уравнение (1) можно записать в виде

Здесь 5 и полукрыло (У>- 0) и его след.

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

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

Построен алгоритм, позволяющий производить расчет с произвольными трапециевидными панелями и следами.

яг

ср( х, у, г, і) = изФ(Х, У, Z) ехр [г (ХА' + шТ)],

0)

и> ехр [— І (IX + и?)]

: Щ

На следе крыла справедливо соотношение

Ф (X, У, О) = Ф [Хт (Г), У, 0) е~ь {Х~Х? <у>! V = ш/Э», (2)

здесь X — Хт(У) — уравнение задней кромки.

Из интеграла Коши — Лагранжа [1| следует, что

(3)

5+^

г± = г(Х, ±У, г), (1,71 £5).

Применены эффективные [квадратурные формулы для вычисления интегралов по крылу и следу.

Точки коллокации выбираются в центре средней хорды панели.

Приведены расчеты плоских колеблющихся крыльев.

2. Алгоритм расчета нагрузок заключается в следующем.

Разобьем полукрыло N—1 прямыми У = Ym = const (т = 1,

N — I) (фиг. 1). Передние и задние кромки полосы У £ \Ут, Ym+1] (т = 0, . . ., N—1) аппроксимируем прямыми с углами стреловидности

xL(Ym+’) - XL(Ym)

у л*+1_ ут

где X = XL{У) — уравнение передней кромки.

, (Yn+1) — XAYm)

arctg- ут+\_^т

Каждую прямую У = Ут разобьем на J частей. Соединяя узлы {Ху, Ут\ с одинаковыми индексами / отрезками прямых, получаем искомое покрытие крыла панелями. Введем обозначения: 5?1 панель с верхним левым углом {Х?, Ут); \Рт — елея полосы [Ут, Ут+1]. Такой способ разбиения позволяет легко учитывать элероны и изломы кромок крыла [1].

На каждой панели функцию Ф (X, У, 0) считаем постоянной

Ф(Х, У, 0) = Ф» (X, Г £5;).

Тогда интегральное уравнение (4) с учетом соотношения (2) преобразуется к виду:

4.UГ (5, ч) = - £ {2 ф? V? + 17т~1 + *2 (foj + frr1)] +

m=0 I j=0

+ ф? (7Го + /Гот“1 + In + !ит~Х )

(5)

с/« 6/

*У-5 -ш?

чГП

Л/+1

1__та1-* -<+г_^г1-* е-^ут+1

у"»+1_1) 1 г">+1

'/+1

_т+1 1

+

хт+1

xj+\

+ 1к ! е-^т+1 йХ [ -[(*?+1 - 5) + (71 - Ут) с!ё 7*+1]-1 х

X'

т+1

^7+1*

X 1 ,т+! [(^т+1 - ч) + та1 - У С1КТ/*,]

гт+1

0+1

- «г?

т

Я-1

0+1

[(У- -г,) + (**+1 - 6) с!е 7«+1] +

ут+1

+

/А

Т}+1

| е~ш^

ут

-шу+1

йУ

+[(*7 - .6) + (ч “ *"”) с4В Т/1]-1 X

х \—^тг- к^т+1 - ч) + (*я+1 ~«)чЯ ■

кт+1

—1кг

-г-рг-[(1”-ч)+ (-*?-«)с%1;] + Дг / «-"''Л' •

1 I

^ “ХМ, К!

£>/72

У

= || 1±**1 <Г‘ 1*+’<*-*,>]

^т+1

/* = вьл?+1 ^ е-^ах [ 1+^.0 д-аг

(5)

Здесь ч?, 7^_, — углы стреловидности передней и задней кромок

панели 5”; Ф/ — значение скачка потенциала скорости на задней кромке; Дт — треугольник, образованный задней кромкой ^(Г) при К® < У <; Кт+1 прямыми У = Ут и Аг = А'у+1.

Выберем точки коллокации (£, -г)) в центре средней хорды панели. Тогда система (2.1) представляет собой систему NJ уравнений относительно ДГ(/ + 1) неизвестных. Замыкающее соотношение

для Фу определяется аналогично в работе [9] аппроксимацией Ф^, пт, 0) в окрестности задней кромки полиномом Лагранжа второй степени с узлами в точках коллокации боксов 5^_2,

и в точке £“ = ХТ(У1т), Т1т = {Ут -\-Ут+1), где Ф(£”, у\т, 0) = Фу.

На задней кромке перепад давления равен нулю и, следовательно, при Х{У1т)-*-Ъ™ дФ/дХ -)- ЬФ -* 0. Подставляя в это соотношение полином Лагранжа, находим искомое замыкающее соотношение:

ф т —

Ф£_2-*»©£,■

1 - Ь (57_2 - X?) - й" [ I - /V (Є»., - *тт)] ’

( е? О-X? У

(6)

\ =/-1 “ Лт }

Рассмотрим вопрос о вычислении интегралов (5). Вычисление ь

интегралов вида J е~ш йх, входящих в I”1, выполняется по двух-

а

точечной формуле Гаусса [12]:

I/(х) йх = -Ц-^- [/{х,) + / (я,)].

X,

Н'-Н)

> X, = а + 2~ (і + р=) .

(7)

Вычисление интегралов и если точка коллокации не принадлежит панели 5™, осуществляется последовательным применением формулы (7) для вычисления внутреннего и внешнего интегралов. Если точка коллокации принадлежит панели 5™, то интеграл содержит особенность О < которая устраняется переходом к полярной системе координат (г, ф) с центром в точке коллокации:

2 *-(?»),

У + 1

б® :

/ ехр [— /Аг(ф)]гіф,

0+1

где /-(ф) — уравнение контура панели 5^*:

г(ф) =

0+1

5іп( т7+і-^+і)

*п(^+1-ф) Гт+1 -ч

віп ф ’

зт(

®т(Т7-ф) Кт — г]

БІЛ ф ’

Ф7+і<Ф<ФуяЙ1. Ф^,1<Ф<Ф7+1, Ф7+1 < ф < ф?, Ф"1<Ф<Ф™.1 4- 2 тс.

Интеграл 0“- вычисляется последовательным применением квадратурной формулы (7).

В применим к внутреннему интегралу формулу интегрирования по частям:

ут+1

Г 1+1Ы 1кгау__ 1 ( +1

J гЗ е | гт+1 в —

у/Я

у/71 + 1 \

I

Применяя к оставшемуся интегралу формулу (7), получим:

V /м "Ь 1

/т __

✓ оо со

/• [Аг"1 + !+ ^д:] г -г[Агот+чх]

Чтп----------------ах~ гт,у ы (IX +

л гт+1 (Л'—ё)з гт(х—ф

■ут+1 ут+1

-V Л/

00

/ Г е-ф7Г+1,Х] ‘ Г ^-фг^+.х]

+ ИгЬтI J (Х _ ?)2 J ■ _;)2 )[ ‘

\Лт+1 хт + 1 /)

Здесь

■ й*= (^ш+1 - Ут), гТ = г (X, 7Г) (г = 1, 2);

п* = У *+&«*■ (1 У?=Ут+Ьт(1+у=у

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

00

е -Цкг+чС)

оо

и,=

Р л ' I (ВТ + Уч)

и

а

]-£—Р------------^ (С = *-5).

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

С со

и, = и; + и;== |/,(*)</*+|/,(.*) (*=1,2).

ас

При С^с г»С и, следовательно,

00 00

С р *-*<* + ») С

и;'=]-^—л, и*=] ——^—л.

с с

Эти интегралы легко вычисляются по формулам:

И1 =4 [ -с-{-т- - **) - а2£(ас).

и;-

— Іа.Е (ас) (а = k + v),

Е (ас) =

= — Ei (— lac) = — Сі (olc) + І Si (ас).

Интегралы И' и И2 заменой переменных т = С + Мг преобразуются к стандартному виду определенных интегралов от быстро осциллирующих функций, которые вычисляются с помощью формулы Филона [12]. Константа с выбирается из критерия [4]:

с = max XT + 4 max (XT - X0m).

m m

Интеграл E(ac) расходится при a 0. Поэтому при малых частотах интегралы Hj и И2 вычисляются по асимптотическим формулам:

Иг

\г (а)

(У-уI)2 I а

НМ

(г-ч)

In

(Y-ri)+r(a)

И,=

3. Данный метод запрограммирован на языке ФОРТРАН-4 для ЭВМ „Минск-32“. На фиг. 2 — 4 приведены рассчитанные распределения Ф и Аср на крыльях, совершающих тангажные колебания около оси, проходящей через середину корневой хорды. Расчеты проведены с равномерным разбиением полукрыла на панели.

0,12

о,в

ЯеФ

/ ^>2

/

1~¥*¥ 2-6*8 3-8*8

1т*

V

0,10 0,08 0,06 0,01 0,02 О -0,02

0,2 0,t 0,6 0,8 х/Ь

Фиг. 2

3—Ученые записки № 2

0,6

4*

0,2

О

-0,2

0 1-¥*V 2-6*6 3-8*8 о-Г71

А ,3 ТтАСр

0

/ь. /а 0 Tie Аср

^оо

0,2

0,£ 0,8 Фиг. 3

0,8 х/Ь

АСр

0,4-

0,2

0 -Д-, 1*2 :т

О 0,2 0,4 0,6 0,8 Х/Ь

Фиг. 4

На фиг. 2 и 3 показаны распределения Ф и Дср соответственно на хорде у!з = 0,15 прямоугольного крыла с относительным удлинением, равным 2, М = 0, ш = 2 при Л^Х •/= 4 X 4, 6X6, 8X8. Легко видеть, что метод обладает быстрой сходимостью по Ф, однако вследствие численного дифференцирования, основанного на кусочно-параболической аппроксимации, сходимость по Дср к результатам расчетов из работ [7, 13] хуже, причем ухудшение сходимости отмечается лишь в окрестности передней и задней кромок.

На фиг. 3 символами „Д“ и „0“ обозначены результаты расчета Дср при ЛГХ-/=4Х4 и 6X6 по формулам численного дифференцирования с учетом особенности распределения Ф: Ф =/(Х) Р(Х, У), где /(X) — весовая функция, и последующей кусочно-параболической интерполяции функции Р(Х, У). Сходимость по Дср значительно улучшается. Результаты, приведенные на фиг. 4 и 5, получены с помощью этой методики.

На фиг. 4 приведено распределение Дср на хорде .у/е = 0 стреловидного крыла с относительным удлинением, равным 3; 71=38°30'; 7т = 52°30'; М = 0,8 при ш = 0 и сравнение с результатами работы [14]. Следует отметить, что при Л^Х-/ = 2Х4- совпадение удовлетворительное. На фиг. 5 приведено распределение Дср по хорде _у/$ = 0,9 прямоугольного крыла Л/? = 3 при М = 0,24, ш = 1,41. Форма колебаний задается соотношением [14]

2 = 0,180431 (уі$) | + 1,702551 (у^у \ 1,136881 (_у/«)31 + 0,253871 (уі^) |.

Совпадение с результатами работ [7, 14] хорошее.

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

ЛИТЕРАТУРА

1. Белоцерковский С. М., Скрипач Б. К., Табачников В. Г. Крыло в нестационарном потоке газа. М., „Наука", 1971.

2. Федорова С. И. К вопросу о точности метода дискретных вихрей при расчете флаттера. „Ученые записки ЦАГИ“, т. 3,

№ 4, 1972.

3. НабиуллинЭ. Н. О двух методах расчета нестационарных аэродинамических нагрузок в несжимаемом потоке. „Ученые записки ЦАГИ', т. 6, № 5, 1974.

4. S t а г k V. J. Е. A subsonic oscillating-surface theory for wings with partial-span controls. „А1АА Pap.“, N 72 — 61, 1972.

5. H a vi 1 a n d J. K., Yoo Y. S. Downwash-velocity potential method for oscillating surfaces. „А1АА J.“, vol. 11, N 5, 1973.

6. M о r i n о L. Unsteady compressible potential flow around lifting bodies: General theory. „А1АА Pap.‘, N 73—196, 1973.

7. Mori no L., Kuo С. C. Unsteady subsouic compressible flow around finite thickness wings. „А1АА Pap.“, N 73 — 313, 1973.

8. Mori no L., Kuo С. C. Subsonic potential aerodynamics for complex configurations: A general theory. „А1АА J.*, vol. 12, N 2, 1974.

9. Jones W. P., Moore J. A. Simplified aerodynamic theory of oscillating thin surfaces in subsonic flow. „А1АА J.‘, vol. 11, N 9, 1973. -

10. Захаров А. Г. Численный метод расчета аэродинамических характеристик плоских и неплоских крыльев в сверхзвуковом потоке. „Ученые записки ЦАГИ“, т. 4, № 3, 1973.

11. La ndahl М. Pressure-loading functions for oscillating wings with control surfaces. „А1АА J.“, vol. 6, N 2, 1968.

12. Хемминг P. В. Численные методы. М., „Наука", 1972.

13. Laschka В. Die instatinaren luftkrafte an harmonisch schwin-genden tragfliigeln endlicher spannweite bei Unter-und Oberschallgesch-windigkeit. „Jahrbuch 1961 der WGL*, 1961.

14. A1 b a n о E., R о d d e n W. P. A doublet-lattice method for calculating lift distributions on oscillating surfacea in subsonic flows. „А1АА J*., vol. 7, N 2, 1969.

Рукопись поступила 4jII 1977 г.