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

Электронный журнал «Труды МАИ». Выпуск № 51

www. ша1. ги^аепсеЛгМу/

УДК 533.6

РАСЧЁТ СТАЦИОНАРНЫХ АЭРОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК ТОНКОЙ НЕСУЩЕЙ ПОВЕРХНОСТИ С УЧЁТОМ ЕЁ УПРУГИХ ДЕФОРМАЦИЙ

Д. А. Болосов

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

аэродинамика; аэроупругость

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

Данная задача относится к классу сопряженных и включает две задачи [1]:

- расчёт поля течения около несущей поверхности и распределённых нагрузок (давления);

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

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

1. Расчёт поля течения и распределённых нагрузок

Для расчёта течения около несущей поверхности используется система стационарных уравнений Эйлера:

dE dF dG л ^

— +-+ — = 0, E =

dx dy dz

pu Pv pw

pu2 + p puv ; F = puv pv2 + p ; G = puw pvw

puw pvw pw2 + p

где t - время;

p - плотность;

(и, V, м?) - проекции вектора скорости V на оси декартовой системы координат (х, у, z);

p - давление.

Рис. 1. Обобщённый алгоритм решения сопряжённой задачи

Уравнение энергии в стационарном случае представляет собой интеграл Бернулли и

для совершенного газа с показателем адиабаты у имеет вид:

2 2 2 Y p u + v + w ---1--= const.

Y-1 P 2

Допущение об отсутствии зон возвратного течения (отрывов) позволяет использовать «маршевый» метод [3] интегрирования с применением схемы Мак-Кормака [2] вдоль бортовой хорды (координаты x):

Шаг предиктор: Е= Е^ 1 -¥)>)-ДХ1 -^),

Ду

Шаг корректор: Е= 0.5

Е + е -Дх (Р+1 - )-Дх (а% - 1)

' 1 к

'М Ду\ М

Шаг по маршевой координате (координате х) определяется из условия устойчивости КФЛ (условия Куранта-Фридрихса-Леви):

Дх = тахДхв,Дхс , Дхв <Ки-

Дп

ТВ

тах А

Дхс < Ки-

ДС

тах А

где Дп, ДС - шаг сетки вдоль осей параметрических координат;

Ав, Ас - собственные значения матриц Якоби системы стационарных уравнений Эйлера, записанной в неконсервативной форме:

АдУ _д¥ _д¥ п

А-+ В-+ С-= 0,

дх дп дС,

Р и

V

w р

где V =

вектор примитивных переменных.

Матрицы Якоби имеют следующий вид:

А =

и 1 1 Р ! 0 ! о 1 ! о 1

Л 1 ±. 1 1 1 1

0 и 0 ! о

---- 1 — 1

0 0 \ 1 и ! 0 1 0 1

0 "ОТ 1 1 0 1 1 и 1 1 ! о 1 1

0 1 2 1 ра 0 ! о 1 1 1 1 и 1 1

в

ПхМ + Пу* 1 ПхР Пу Р 0 0

1_

0 1 ! пхи + Пу* 1 ' 0 0 1 пх-р _ 1

0 ! 0 1 пхи + Пу* 0

0 ! 0 1 0 пхи + Пу* 0

г г г г

0 1 2 \ ПхРа 1 Пх ра2 0 пхи + Пу*

с =

Си+Сг™ J Сх Р 1 1 1 1 1_ 0 J с р \ 0 1

0 Схи + Сг™ 1 1 1 0 0 ! с - 1 Ьх

_________ 1 1— ______ !____ _р___

0 0 1 + Сг™ 0 ! 0 1

0 0 1 1 1 0 Схи+Сг™ 1 ^ х

1 ! р

0 Сх Ра2 г 1 1 1 1 п 0 С Ра' "Г Схи + сг™ 1 _

Приведём систему стационарных равнений уравнений Эйлера в неконсервативной форме к виду: ЗУ ЗУ ЗУ

Зу + в * — + с — = 0.

дх Зу Зх

Матрицы В *, С* определяются по зависимостям: В* = А-1В, С* = А-1С, или в развёрнутом виде:

Пи - *

В*

I

VI цм - V

— 1 о—-

\ У 2 2

и \ и - а

--1—

I

2

0 ! -

Пха - иV

2 2 и - а

0 | ра

2 Пи - *

2 2 и - а

ПуР

и

22 и - а

Пуа

! 0

0

2 2 1 I- и - а _[

V и

0 0 1 1 1 0 1 ~ 1 1 V 1 1 — 1 0

1 1 ! и !

ПуР

и

22 и - а

! 0

и(и2 - а2) рри2 - а2_)

Пу_ ри

Пха - и*

2 2 и - а

0

0

C* =

w i Zxu - WW

u

P

2 2 u - a

1 > 2" -""1 i Zxa - uw

0 i -

22 u - a

0 !

0 !

0

0 ' pa:

Cxu - W

2 2 u - a

0 0

w u

0 0

zp

u

22 u - a

j Zxu - W

I ../ 2 IT

I

u (u2 - a2)

Zza2

22 u - a

'i i

i p(u2

Cxu - W

(u 2-a 2J_

WW

u

Pu

С p2

22 u - a

1 ^ 2 -i Zxa - uw

I -jUx-

I 2 2

i u - a

где V = цхи + цуу, й = £хи + С^ .

Найдём собственные значения матриц В *, С*, используя зависимости

B *-AB E

= 0,

C*-AC E

= 0, в результате получим:

в , v в , uv ± asl

¿1,2,3 =Пх +ПУ~, Аз =Пх + Пу

u

uv ± a^v2 + u2 - a2

22 u -a

¿2,3 =Zx + С,~, л:,з =Zx +Zz u

■ ±

uw ± aylw2 + u2 - a2

22 u - a

В данной задаче несущая поверхность полагается бестелесной (имеющей нулевую толщину профиля). В плоскости x = const строится прямоугольная сетка, границы которой удаляются на необходимое расстояние, исключающее влияние последних. Для возможности рассмотрения несущих поверхностей различной формы в плане геометрия несущей поверхности не выделятся, а требуемая точность вычисления суммарных (аэродинамических характеристик - АДХ) и давлений на наветренной и подветренной стороне достигается заданием требуемого количества узлов вдоль осей z и У . В соответствии с этим граничное

условие «непротекания» [3] (V ■ nn = 0, где V - вектор скорости, nn - единичная нормаль к

поверхности крыла) определяется только в тех узлах, которые локализованы контуром несущей поверхности.

В результате «маршевой» процедуры решения исходной системы уравнений получаем область, показанную на рис. 2а.

0

0

0

u

Рис. 2. Сетка в результате решения

Для верификации изложенного метода рассмотрим прямоугольную несущую поверхность (передняя и задняя кромки образованы линиями х = 0 и х = 10 соответственно, а концевая хорда - линией, г = 10 линия симметрии в данном случае г = 0) с удлинением X = 2 (рис. 3а), расположенную под углом атаки а = 5° в потоке с числом Маха М = 2 .

Согласно теории идеальной жидкости на поверхности такого крыла образуется две зоны (I, II). Зона I, расположенная вне характеристического конуса (конуса Маха), характеризуется постоянной разностью давлений, определяемой зависимостью: Ар = 2и2 • р-а ■ tgв, где и - скорость набегающего потока; р - плотность набегающего потока; а - угол атаки;

в - угол при вершине характеристического конуса. Зона II ограничена конусом Маха (рис. 2а) с вершиной в начале концевой хорды.

-о-о-о-о маршевый метод на основе уравнений Эйлера [3] <ккн> теория идеальной жидкости [4] а) зона влияния концевой хорды б) сопоставления расчёта распределения разности

давлений

Рис. 3. Разность давлений наветренной и подветренной стороны консоли прямоугольной несущей поверхности (X = 2, а = 5°, М = 2)

Рассматривая течение около треугольного крыла, выделяют два основных случая: обтекание крыла с дозвуковой передней кромкой (рис. 4) и сверхзвуковой передней кромкой (рис. 5).

10 о

■маршевый метод на основе уравнений Эйлера [3] ■теория идеальной жидкости [4] а) зона влияния концевой хорды б) сопоставления расчёта распределения разности

давлений

Рис. 4. Разность давлений наветренной и подветренной стороны консоли треугольной несущей поверхности с дозвуковой передней кромкой (%ПК = 60 °, а = 5°, М = 1.5)

10 о

-о-о-о-о маршевый метод на основе уравнений Эйлера [3] оооо теория идеальной жидкости [4] а) зона влияния концевой хорды б) сопоставления расчёта распределения разности

давлений

Рис. 5. Разность давлений наветренной и подветренной стороны консоли треугольной несущей поверхности со сверхзвуковой передней кромкой (%ПК = 600, а = 5°, М = 3)

В случае дозвуковой передней кромки на крыле имеется всего одна зона, в которой давление возрастает от корневой хорды к передней кромке (рис. 4б). В случае несущей поверхности со сверхзвуковой передней кромкой (рис. 5), как и на прямоугольном крыле, образуется две зоны: зона II, располагающаяся между передней кромкой и границей конуса Маха, характеризуется постоянным давлением, определяемым по зависимости, указанной выше. Распределение давления в зоне I аналогично случаю несущей поверхности с дозвуковой передней кромкой.

Были проведены систематические расчёты прямоугольных несущих поверхностей. Результаты расчётов и экспериментальные данные [4] представлены на рис. 6. Расчеты по углу атаки проводились с шагом 0,5о. При использовании сеток с количеством узлов в сечении по продольной координате 51x51 и больше результаты по интегральным аэродинамическим характеристикам практически не изменяются.

Су

0.20

0.15

0.10

0.05

о.ооу-

- □

- о

- и/

- □ / У / Э

Эксперимент О к= 0,5 О Ь=1,0 □ ^=3,0 1,1, -

| | | |

7 5 6 а, град

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.0С|

- 1 1 1 Эксперимент О Хпк= 45° Л Хпк=60° □ Хпк=70° <> Хпк=75° _

- \

- Л

- э

- □ о о

- о __2

-

_1_ _1_ _1_ _1_ _1_ _1_ _1_ _1_ _1_ _1_ _1_

.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 М

Рис. 6. Сопоставление с известными экспериментальными данными [5] результатов расчёта прямоугольных несущих поверхностей для числа Маха М = 1,7 (маркеры соответствуют экспериментальным данным, сплошные линии - результатам расчётов)

Рассмотренный метод позволяет проводить расчёты АДХ несущих поверхностей различной формы в плане для диапазона чисел Маха М = 1.5 - 5 и углов атаки а = 0 -10°. По результатам расчета получаются не только суммарные АДХ, но и распределённые нагрузки.

2. Расчёт напряжённо-деформированного состояния несущей поверхности под действием распределённых нагрузок

Для упрощения расчёта напряжённо-деформированного состояния и сведения его к двумерному представим несущую поверхность как тонкую пластину и введём допущения (гипотезы Кирхгоффа) [6]:

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

- малость нормальных напряжений в площадках, параллельных срединной поверхности;

- малость толщины пластины по сравнению с её размерами;

- материал пластины однородный, изотропный и подчиняется закону Гука.

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

__ р д4о „ д4о д4о р _ Ек3

УУс =— или —- + 2 2 2 +—- = , В = —г-п,

В дх4 дх 2д22 д24 В 12(1 )

где а - прогиб;

р - давление несущей поверхности;

В - изгибная жесткость пластины;

к - толщина пластины;

/ - коэффициент Пуассона;

Е - модуль Юнга.

Граничные условия для жёстоко закреплённой по бортовой хорде несущей поверхности имеют вид:

л дс

о = 0 и — = 0 - для жёстко закреплённого края;

дп

да д а д

—Т + и—Т = 0 и — дп дs дп

д2 о /„ \д2 а —- + (2-/)—-дп дs

= Р - для свободного края.

Указанное выше уравнение решается методом конечных элементов [7], для чего несущая поверхность разбивается на заранее заданное количество треугольников, число которых определяется необходимой точностью получаемых результатов. Перемещения пластины можно представить в виде:

°=Л ^р^+в ц+А ^ ^+21, ^ )+£ ^ +2 )+

+ в (ЬЬ + 2 Ь ¿2 Ь ) + в (Ь Ь + 2 Ь Ь 4 ^ + в (ЬЬ + 2 Ь ¿2 Ь ) + + в (Ц3 Ц2 + | Ь1Ц2Ц3 |

где Ь1,Ь2, Ь3 - I - координаты треугольного элемента.

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

а1 + Ъ1 х + с1 у а 2 + Ъ2 х + с2 у аз + Ъз х + сз у ^ 1

1 = 2А ' 2 = 2А ' 3 = 2А , = 2

а1 = Х2 уз Х3У2 ,

Ъ1 = У 2 - Уз,

с1 = хз Х2 ,

1 х1 у1

1 х2 у 2

1 хз уз

где Д - площадь треугольного элемента,;

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

}, Ох, = -

{ д} \дУ л

и О ,

где Ох{, Оу{ - углы поворота в узле г вокруг осей х и у соответственно.

Окончательно функцию формы для узла (вершины) 1 треугольного элемента можно записать в следующем виде:

к Г =

-¿1 + ^2 ¿3 Ьз

Г 2 1 1 Г 2 1

Ъ\ Ь, Ь + — ЬЬЬ I- Ъ\ ЬЬ, +- " Хл

2

2

Г 2 1 11 2 1

" Хл Хл

2 V \ 2 Функции для узлов 2 и 3 получаются циклической перестановкой индексов. На рис. 7 изображена геометрия несущей поверхности и её разбиение на конечно-элементную сетку.

На рис. 8, 9 представлен пример расчёта напряжённо-деформированного состояния для данной жёстко закреплённой по бортовой ходе несущей поверхности постоянной толщины 0,65 мм, нагруженной равномерно распределённым давлением 0,1 МПа.

140

Рис. 7. Геометрия несущей поверхности и конечно-элементная сетка

0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 О.вО 0.65 О.ЭО

Рис. 8. Распределение прогиба несущей поверхности (см) от равномерно распределённого давления 0,1 МПа

Рис. 9. Деформированная несущая поверхность

3. Результаты расчёта стационарных АДХ несущей поверхности с учётом её деформаций

Расчёт аэродинамических характеристик проводился для числа Маха М=2 в диапазоне углов атаки 0 - 12° для геометрии, изображённой на рис. 7.

На рис. 10 представлены основные аэродинамические характеристики несущей поверхности.

—о—э— без учета деформаций —□—□— с учетом деформаций

Рис. 10. Аэродинамические характеристики несущей поверхности в зависимости от угла атаки (коэффициент подъёмной силы Сп, коэффициент момента тангажа Ы2, центр давления Хй )

Как видно из рис. 10, деформации несущей поверхности оказывают влияние на её аэродинамические характеристики, а следовательно, и на аэродинамические характеристики летательного аппарата в целом.

Библиографический список

1. Raymond E. Gordnier, Miguel R. Visbal. Computation of three-dimensional nonlinear

panel flutter. AIAA A0116437 January 2001.

2. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: В 2-х т. Т. 2: Пер. с англ. - М.: Мир, 1990. - 392 с.

3. Аэродинамика ракет: в 2-х Кн. 2. Пер. с англ./Под ред. М. Хемша, Дж. Нилсена. - М.: Мир, 1989. - 512 с.

4. Красильщикова Е. А. Тонкое крыло в сжимаемом потоке. - 2-е изд., доп. - М.: Наука.

Физматлит, 1986. - 288 с.

5. Аэродинамика, устойчивость и управляемость сверхзвуковых самолётов / Под ред.

Г. С. Бюшгенса. - М.: Наука. Физматлит, 1998. - 816 с.

6. Бояршинов С. В. Основы строительной механики машин: Учебное пособие. - М.:

Машиностроение, 1973. - 456 с.

7. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 543 с.

Болосов Дмитрий Александрович, инженер-исследователь ГУП "Конструкторское бюро приборостроения,"

300001 г. Тула, ул. Щегловская засека, д. 59, т. (4872) 46-94-64, е-mail kbkedr@tula.net