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

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Том 150, кн. 1

Физико-математические пауки

2008

УДК УДК 532.5.031:533.6.011.32:533.694.2

ОБ ОДНОЙ ЗАДАЧЕ ОПТИМИЗАЦИИ ДВУХЭЛЕМЕНТНОГО КРЫЛОВОГО ПРОФИЛЯ

П.А. Волков, Н.Б. Ильинский

Аннотация

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

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

Введение

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

Предложенный метод построения формы двухэлементного крылового профиля (биплана, профиля с закрылком или предкрылком) с максимальной подъемной силой в классе безотрывно обтекаемых профилей основан на решении обратной краевой задачи аэрогидродинамики (см.. например. [1]). Суть метода состоит в отыскании в рамках выбранной математической модели «оптимального» распределения скорости, заданного вдоль поверхности неизвестного контура профиля и соответствующего максимально возможному значению коэффициента подъемной силы на искомом профиле.

Традиционные способы (см.. например. [2]) оптимизации заключаются в многократном решении прямой задачи и подборе крылового профиля с желаемыми аэродинамическими характеристиками. Метод, предложенный в нашей работе, позволяет непосредственно определить оптимальную форму профиля из фиксированного класса контуров по параметрическому распределению скорости вдоль его поверхности. Углы атаки и положение профилей друг относительно друга при этом не задаются, а находятся в процессе решения. Модель течения построена таким образом, что появляется возможность пренебречь эффектами сжимаемости внутри пограничного слоя.

1. Постановка задачи

В физической плоскости г (рис. 1, о) искомые непроницаемые крыловые профили АкВк (к =1, 2), составляющие двухэлементный профиль биплана, обтекаются безотрывно плоским установившимся потенциальным дозвуковым потоком с заданными числами Рейнольдса Дето и Маха Мто та бесконечности. Контуры Ьк

двухэлементного профиля: индексы 1 и 2 соответствуют Li и L2

профилей считаются гладкими за исключением задних кромок Bk , где внутренний к области течения угол равен 2п.

Направление оси абсцисс в физической плоскости z выбрано по направлению скорости vTO набегающего потока. Периметры профилей известпы и равны lk (k = = 1, 2).

При условии безотрывного обтекания и малости толщины пограничного слоя дуговые абсциссы sk контуров профилей и полу тел вытеснения на участках Bk 'AkBk " будем считать совпадающими и отсчитывать их от пуля в точках Bk' до lk в Bk'' так, чтобы щи возрастании sk область течения оставалась слева.

Распределения приведенной скорости Ak вдоль участков Bk 'Ak Bk " полутел вытеснения в параметрическом внде есть

A k = vk (sk, dj )/a* = A k (sk,dj), sk G [0, lk], k = 1, 2, j = 1,..., 5,

где dj - свободные параметры, a* - критическая скорость (рис. 1,6). Функции Ak -кусочно-гладкие, обращающиеся в нули в точках Ak разветвления потока sak, и в этих точках оии непрерывно дифференцируемы. Распределения приведенной скорости выбираются из класса гидродинамически целесообразных (ГЦРС) [3], что обеспечивает безотрывное обтекание в рамках принятой математической модели. Эти распределения берутся с так называемой «полкой» (рис. 2, а), то есть с участком постоянной скорости, так как именно «полочные» распределения являются экстремальными в задаче максимизации площади эпюры распределений скорости.

Задания распределений скорости вдоль искомого контура не достаточно для получения физически реализуемого решения профиля заданной толщины, ограниченного замкнутым контуром без точек самопересечения. Дополнительные геометрические ограничения, налагаемые на Lk (рис. 2, 6) и позволяющие зафиксировать класс искомых решений, имеют вид:

min ck (t) > c', (1)

te[xki, xk2]

c''' > max ck(t) > c'', (2)

te[xfci, xfc2]

ck ( t) k

от хорды этого профиля в точке x* = t; участок хорды [xki, xk2] и величины c', c'', c''' задаются. Первое из ограничений выражает условие отсутствия точек самопересечения контура, второе эквивалентно определению границ изменения максимальной толщины профиля.

а)

Л

А

<¿1

/

в >

0.1 0.0

ш

0.4 0.6

х*

Рис. 2. Схема (а) задания свободных параметров = 1,..., 5) в распределениях

приведенной скорости; геометрические ограничения (, б), налагаемые на Ь (к = 1, 2)

В соответствии с принятой математической моделью течения существует комплексный потенциал потока т(г) = <р(х, у) + гф(х, у), где г = х + гу. Считаем, что известен расход ч* между контурами Ьк полутел вытеснения и разность потенциалов между точками Л2 и Л\. Положив т(г) равным нулю в точке Л\, будем иметь:

^(«а1)=0, 'ф(8а1)=0, ) = <£>*, Ф($а2) = Ч* ■

Требуется найти форму двухэлементного крылового профиля, обладающего максимальной подъемной силой.

В качестве характерного линейного размера выберем величину I = 0.5(11 + 12) полусуммы периметров искомых профилей.

Учитывая безотрывный характер обтекания, воспользуемся предположениями, что для больших (порядка 105 — 107) значений числа Дето вязкость будет сказываться лишь в достаточно тонком пограничном слое около профиля, а для дозвуковых скоростей внешнего потока сжимаемостью внутри слоя можно пренебречь. Решение обратной задачи сведется при этом к нахождению контуров полутел вытеснения. обтекаемых потоком газа Чаплыгина, зная которые находим форму искомых профилей, отступив па участках Вк 'ЛкВк " внутрь полутел на толщину вытеснения 6*(з).

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

2. Решение обратной задачи

в

У |А*(т)|ь-1 ат + /к > /с,

где /к(в) - формпараметр, /о, а, Ь - эмпирические постоянные, характеризует вклад в /к ламинарных участков. Если пограничный слой на профиле полностью турбулентный, то /¿к = 0, а = 1.17, Ь = 4.75, /0 = —2.

В качестве канонической области во вспомогательной плоскости и = £ + щ выбран прямоугольник со сторонами и ^ = гп/2. Область изменения комплексного потенциала в плоскости т для газа Чаплыгина сохранит тот же вид, что и для идеальной несжимаемой жидкости [4]:

ги(и) = <£>(£, г)) + гф(£, 1]) = А((и — га) + А((и + га)+

Г + Гт а (и — га)

+ ^тг^ 1п + +«71 + /?,

2пг а (и + га)

-—(■и) = П(и) = —Ар(и — га) — Ар(и + га) -\--——т—[С('ы ~~ ~~ С('ы + + 71;

аи 2пг

<Р1(£1) = ^(£1,пг/2), Ф1(£1) = Ф(£1,пг/2), ) = 0), фт(Ь) = Ф(Ь,0), ¿1,6 е [0,^1],

где ¿1,^2 е [0,^1^ Л = Л1 + гЛ2', в = в1 + гр2, а - неизвестные постоянные; р(и), ((и), а (и) - функции Вейерштрасса с полупериодами пг/2ш ш1/2] п1 = ((и+ + Ш1) — ((и) - константа, зависящая от ^¡циркуляции скорости Гк по каждому

1

а ?i\ —

контуру Lfc равны fk(lk) - fk(0): 71 = — -Г1 + (Г1 + Г2)--(А + А)щ

ш>1 L п

константа, зависящая от ш1, а.

Модели газа Чаплыгина соответствует вспомогательная функция

х(и) = in-2|Л('Ы)| - гв, с2 = 0.296,

1+ у/1 + 4с2 А2 («)

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

х(и) = Х*(и) +

. п(и - Ual) . п(и - Ua2) sin- sin-

Ш1 Ш1

аналитическую функцию х*(и) = $ * — гв*, действительная часть которой на границах прямоугольника известна:

1+ v/l + 4c2A2(Sfc(6))

- In

тг(бг - Cal - Нк - 1)тгг) . - £a2 + ^(2 - k)ni) sin--- sin--

Ш1 Ш1

восстановим по формуле Вилля (см., например, [6]):

W1

П1 - 2

X» = - \ [Si (£К (« -i-Ьг)- S2(C)C Ы - О] de + + т,

пг J 2 2ш1

о

где P2 - произвольная вещественная постоя иная, а для P1 имеем

UJl UJi

P1 = j S№ dt = j s2(t) dt (3)

оо

Условие (3) является условием однозначности функции х*(и)-

Координаты контуров полутел вытеснения находятся интегрированием выражения

(¿г = ехр[—х(г()]№ '(и) с1и — с2 ехр [%(«)] го '(и) с1и, (4)

где черта сверху означает комплексное сопряжение. Форму искомых профилей находим, отступив по нормали на участках В к 'ЛкВ к '' внутрь полутел на величину толщины вытеснения 3*(в). Функция 3*(в) определяется при помощи любого из известных методов, например, однопараметрического метода Кочина Лойцянского [5].

Контуры Ьк искомого двухэлементного профиля будут замкнуты, если соответствующие контуры полутел вытеснения, определяемые из (4), разомкнуты на величину

Агк = —г^ок ехР (г^ок), (5)

где 0ок - аргумент вектора скорости в точках Вк'и Вк ''; 3*к = 3*(0) + 3* (1к); 3*(0), 3* (Iк) — значения толщин вытеснения в точках В к ' и Вк '' соответственно. Кроме того, должно выполняться условие (3) однозначности функции х* и условие совпадения величины приведенной скорости набегающего потока, найденной в ходе решения задачи, с заданным значением Ато = уж/а*:

Иех(га) =1п{2Ато [1 + {1 + 4с2А^}1/2р (6)

3. Нахождение аэродинамических характеристик

Коэффициенты сопротивления Схк и подъемной силы Сук каждого из профилей двухэлементного профиля определяются формулами

Сх к = Ф (cf к ¿х + ср к ¿у), Су к = Ф (cf к ¿у — ср к вх), ■у Ьк ■■'Ьк

где с^к - местный коэффициент трения, срк - коэффициент давления для адиабатического течения

ср к

^0

М2

1-Л^/[(^0 + 1)^0-1)]

щ = 1.41 для воздуха.

Сх Су

намическое качество К двухэлементного профиля находятся по формулам:

Сх = Сх1 + Сх2, Су = Су1 + Су2, К = Су/ Сх-

4. Решение задачи оптимизации

Это решение состоит в оптимальном выборе исходных распределений приведенной скорости А к = А к (в к,вз), при которых функционал Г (¿^) = Су (А к (в )) принимает максимальное значение с учетом условий разрешимости обратной задачи (3), (5), (6), геометрических ограничений (1), (2) и дополнительных линейных ограничений на свободные параметры:

Л1в = дь (7)

Л2в < д2,

(8)

0 < а < (9)

где матрицы А1? А2 и векторы 01, д2, 0, 0т € Дп заданы.

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

) = е1Дг + е2Дс* + е3Дс** - Су, (10)

где Дг - сумма квадратов невязок условий разрешимости обратной задачи, а величины Де*, Де** отыскиваются из выражений:

Де* = Д1е* + Д2е*,£Де** = Д1е** + Д2е**,

Ак о*

(ок - о') > ок < о 0,

к —

(ок* - о 'А2

fc -L ) > Lfc

< ок <

^о** = <0, о'' < ок* < о''',

,(ок* - о '")2: „

ок = min о к (t), ок * = max о к (t).

te[xkl,xk2] te[xkl,xk2]

Для заданных величин штрафа £1 > 0, £2 > 0 и £3 > 0 решается задача минимизации функционала (10) с линейными ограничениями (7) (9). Численная реализация решения задачи осуществлена на языке Фортран при помощи процедуры. основанной на алгоритме «TOLMIN» [7].

Решение задачи минимизации (10) выполняется многократно для различных значений штрафа £1, £2, £3- Эти величины подбираются так, чтобы уменьшение значений функционала достигалось (в соответствии с поставленной задачей) за счет увеличения коэффициента подъемной силы Cy. Затем построенные контуры с максимальным Cy, которые, как правило, получаются разомкнутыми, замыкаются путем постепенного увеличения ^.Одновременно увеличивая £2, £3, добиваемся выполнения геометрических ограничений (1), (2). При этом следует из-

Cy

описанного процесса подтверждена численными экспериментами.

5. Результаты расчетов

Для численной реализации предложенного метода решения оптимизационной задачи разработан комплекс программ. Рассматривалось обтекание двухэлементного профиля полностью турбулентным сжимаемым потоком при Мто = 0.3 и Дето = 1.45 • 107.

В качестве исходных задавались параметрические распределения А & = = А & (в ), к = 1, 2, приведенной скорости «полочного» типа из класса ГЦРС. Из соображений безотрывности обтекания осуществлялся выбор длины полки (участок А = Атах) и достраивалось распределение А & на участке торможения. Оптимальный выбор А & осуществлен подбором свободных параметров (7 = 1,..., 5) (рис. 2, а).

Пример численного решения оптимизационной задачи, когда геометрические ограничения (1, 2) не учитывались, приведен на рис. 3 (штриховая кривая на рис. 3, а соответствует распределению приведенной скорости А 2 = А2/Ато основного профиля, сплошная кривая - распределению А1 = А1/Ато закрылка). Ясно,

Рис. 3. Распределения приведенной скорости (в) и контур (б) оптимального двухэлементного профиля без учета геометрических ограничений

Табл. 1

Л' профиля Ь1 С*1 сГ ъ2 <=2 сГ Су К

1 0.19 0.2 7.4 0.78 0.7 5.9 0.0086 1.54 179

2 0.19 1.0 15.8 0.76 6.8 20.2 0.0130 1.47 113

а)

1.0} 0.0

\

\ ч ч >

-1.0 Г----^

I___.__]_

0.4

" Г "

_ ± _ 0.6

в

- 4 -

У ^ б) „т—-о 2

0.2^— у —1-----

I V ! г--.-

0.1 I------- —___!___----1

!■ 1 1 1

X

-0.6

_ ± _ -0.4

Рис. 4. Распределения приведенной скорости (в) и контур (б) оптимального двухэлементного профиля с учетом геометрических ограничений

что физическая реализация таких тонких профилей на практике затруднительна. Действительно, в результате оптимизации максимальная толщина каждого из элементов профиля не превысила 8% хорды элемента, а с приближением к задней кромке их толщина снизилась и составила в точке 0.85 хорды 0.2% и 0.7% соответственно (см. табл. 1. профиль 1).

Рис. 6. Расчетная сетка: области, прилегающие к контурам

На рис. 4 приводен результат решения оптимизационной задачи в полной постановке для о ' = 1%, о " = 10% о = 25%, xk1 = 0.1, xk2 = 0.85, к = 1, 2. Здесь основной профиль и соответствующее ему распределение А 2 = А2/Ато изображены штриховой кривой, а закрылок и распределение А[ = А1/А1Х> - сплошной.

Полученный двухэлементный профиль удовлетворяет всем требованиям класса искомых решений. Увеличение толщины профиля до приемлемых значений привело к снижению коэффициента Су подъемной силы на 4.5% по сравнению с профилем 1 (рис. 3) и значительному (более 50%) увеличению коэффициента Сх сопротивления, что негативно сказалось на величине К аэродинамического качества (см. табл. 1, профиль 2).

6. Тестовый расчет

Для обоснования достоверности полученных результатов был проведен прямой расчет обтекания оптимального профиля 2 (рис. 4) установившимся дозвуковым потоком вязкого газа. Дискретизация области течения осуществлена структурированной мультиблочной сеткой с прямоугольными ячейками (рис. 5, 6: размеры указаны в хордах: Ь = Ь1 + Ь2 ).

Рис. 8. Сравнение результатов расчета cpk с данными программы FLUENT

Характеристики набегающего потока следующие: MTO = 0.3, Дето = 1.45 • 107, плотность = 1.29 кг/м3, деление = 101325 Па, температура = = 273 К, скорость = 100 м/с, динамический коэффициент вязкости = = 1.7894 • 10-5 кг/мс-1, характерный линейный размер l = 2 м.

Расчет проводился для однопараметрической модели турбулентности Спалар-та Аллмараса с использованием программы FLUENT. Результаты расчета представлены на рис. 7, 8.

На рис. 7 изображены линии тока.

Расчет поля чисел Маха показал, что в областях разрежения давления, примыкающих к верхней поверхности основного профиля и закрылка и соответствующих «полочным» участкам распределений приведенной скорости Хк, течение — дозвуковое, а величина M варьируется в пределах 0.55 ^ 0.58.

На рис. 8 приведено сравнение графиков коэффициентов давления cpk, k = 1, 2, полученных из решения обратной задачи, с данными, вычисленными при помощи программы FLUENT (значения cpk, найденные при помощи программы FLUENT, отмечены точками). В целом графики удовлетворительно согласуются. Некоторое рассогласование коэффициентов давления cp1 на «полочном» участке закрылка объясняется влиянием «следа» за основным профилем на характеристики закрылка, которое в обратной задаче не учитывалось. Последнее обстоятельство, хотя и является одной из причин несовпадения величины коэффициента сопротивления Cx = 0.0253 с величиной Cx = 0.0130, полученной ранее в ходе решения задачи, не оказало значительного влияния на расчет коэффициента подъемной силы Cy = 1.44 (разница в 2% по сравнению с Cy = 1.47, см. табл. 1).

Заключение

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

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

Summary

P. A. Volkov, N.B. Ilinskiy. About t.lie Problem of Two-Element Airfoil Optimization. A numerical solution is presented for the problem of finding an optimal two-element, wing airfoil shape with maximum lift, coefficient, flowed by a subsonic viscous non-separation stream. The solution method is based 011 the theory of inverse boundary-value problem of aerohydrodyuamics for analytical functions. The effects of viscosity and compressibility within the framework of boundary layer theory and the Cliaplygin's gas model are taken into account..

Key words: optimization, boundary layer, Cliaplygin's gas, inverse boundary-value problem of aerohydrodyuamics, two-element, airfoil.

Литература

1. Елшаров A.M., Ильинский И.Б., Поташсч А.В. Обратные краевые задачи аэрогидродинамики. М.: Наука, 1994. 440 с.

2. Елшаров A.M., Ильинский И.Б., Поташеа А.В. Обратные краевые задачи аэрогидродинамики // Итоги пауки и техники. Механика жидкости и газа. М.: ВИНИТИ, 1989. Т. 23. С. 3 115.

3. Степа,нов Г.Ю. Гидродинамика решеток турбомашип. М.: Физматгиз, 1962. 512 с.

4. Седов Л.И. Плоские задачи гидродинамики и аэродинамики. М.;Л.: Гостехиздат, 1950. 444 с.

5. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1987. 840 с.

6. Ахиевер И.И. Элементы теории эллиптических функций. М.: Наука, 1970. 304 с.

7. Бар'теиьев О.В. Фортран для профессионалов. Математическая библиотека IMSL. Ч. 2. М.: Диалог-МИФИ, 2001. 320 с.

Поступила в редакцию 11.12.07

Волков Павел кандидат физико-математических паук, младший научный сотрудник отдела краевых задач НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.

Ильинский Николай Борисович доктор физико-математических паук, профессор, заведующий отделом краевых задач НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.

Е-шаП: Nikolmj.IlinskiyQksu.ru