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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XIX 1 9 8~8

М 5

УДК 629.735.33.015.3.025.73

МЕТОД ЛОКАЛЬНЫХ ВАРИАЦИЙ ПОВЕРХНОСТИ В ЗАДАЧАХ ОПТИМИЗАЦИИ ФОРМЫ ПРОФИЛЯ

Н. П. Бузоверя, Д. П. Кроткое

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

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

Начало развитию данного подхода положили работы Хикса и его-сотрудников [1—3].

В ранних работах (1, 2] для описания контура профиля использовались полиномы, при этом коэффициенты полинома являлись проектными переменными. В последующих работах, например [3], полиномы были заменены более общими аналитическими функциями для описания форм вариаций контура профиля. Аналитические функции для описания форм вариаций контура, но несколько другого вида, использовались также Ляпуновым С. В. в задаче построения профилей минимального волнового сопротивления ¡[4]. В этих случаях геометрия контура профиля определялась следующим образом:

У(х) = Уо (*) + 2 а)

)

Здесь х, у — соответственно продольная и поперечная координаты профиля, отнесенные к его хорде; функция Уо{х) определяет форму базового профиля; а, — весовые коэффициенты, используемые в качестве проектных переменных в программе оптимизации; ¡¡(х)—функции, описывающие формы вариаций контура базового профиля (в дальнейшем именуемые функциями формы).

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

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

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

Аналитическое описание контура профиля осуществляется в классе параметрических сплайн-функций со сглаживанием в некотором переменном коридору А (я), где в — параметр, по которому строится сплайн (например, х (я) = з!п2 5/2, —тс-^О) [5]. Алгоритм сглаживания основан на определении методом локальных вариаций коэффициентов сплайна, минимизирующих функционал

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

где А] — амплитуда у-й возмущающей функции формы; £ (в) — текущее значение длины, отсчитываемой вдоль контура профиля от задней кромки верхней поверхности; — координата управляющей точки (точки максимальной амплитуды); ¿+(в|), —

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

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

J=<§ (у" (*))2 ¿5,

/Л* 00) =

при в <

при 5;- < Я < хт ,

сглаживания, рис. 1. Получаемая при этом вариация формы профиля определяется следующим образом:

§Уу (* (*)) = У) (х («) - Уо (х («)) — ajfj (х («)).

Здесь у у {х (з)— полученные после сглаживания ординаты возмущенного контура профиля; у0 (х ($)) — ординаты базового профиля; а1 — весовые коэффициенты.

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

Их вид зависит от координат точек управления Xj(Sj), _У;($у), выбранных на контуре профиля, границ распространения возмущения справа х/ (5/), у^(з^) и слева х]~ ($У), уу («7) от управляющей точки, от амплитуд А}, от коридора сглаживания и, наконец,

от самой геометрии базового профиля.

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

Задача проектирования оптимальной формы профиля сводится к задаче нахождения условного экстремума функции нескольких перемен-

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

срг (в)= / = 1,2,..../, 1

?i(a)>ci, i = / Н- 1, р, )

и доставляющего экстремум целевой функции

min (max)/(a).

{aj)

В дальнейшем для определенности будем считать, что рассматривается задача на минимум.

Введем малое положительное число е и будем считать, что точка

а является допустимой, если она удовлетворяет ограничениям задачи

(1) с заданной точностью

|срДа) — с,|<е, /= 1 , 2,..., I, j

?,(а)>с, — в, i = l+ ]

Поиск оптимума начинается с произвольной точки а, лежащей в допустимой области. Если такая точка не принадлежит допустимой области, то осуществляется вход в допустимую область (смотри ниже). Обозначим через Ма множество активных ограничений (I ?i (о) — ci I s, IKiKp), в которое войдут все ограничения типа равенств и некоторые ограничения типа неравенств.

Для перехода на т-й итерации от точки а (m> в допустимой области к точке а(т+1) в а(т) находится направление «S(m),

C(m)_ Vf(tt(m))

° — I I ’

где F(a(m'>) = I(afi (ß(m)) — функция Лагранжа.

Здесь MW— множество активных ограничений на т-й итерации. Если точка а(ш) не лежит строго на границе множества (1), то осуществляется поворот вектора вокруг этой точки к границе области (1). В этом случае множители Хг, определяющие направление вектора S{m\ находятся из условия

(yF («(«»)) уш (а<т))) = а8<р? (a(m))t q (. (3)

где

8?в(а(в,)) = <Р,(«(т)) —с,. д<1,

' С .

, / + 1 <д<.р,

г. („с.,=|’>'<а"))-с*при

1 0 при <р?(а(т))cq

а а — эмпирический весовой коэффициент.

В развернутой форме соотношение (3) имеет вид

X ХI (V?/ («(т)) • У?д (0(т))) = - № (а(т)) ■ уср9 (а(т))) + а8?? (а<т>),

?^Г>- (4)

В направлении вектора 5(/п) делается очередной шаг от точки а{т) к точке а(т+1) = а(т) + А<т)5(т). Величина шага &(т) определяется в процессе оптимизации и равна тт{№™\ ^т)), где — значение

шага, при котором луч а{т) + 5(ш) пересекает область допустимых значений (2), а —значение шага, при котором точка а</п+1)

принадлежит области допустимых значений (2) и целевая функция достигает минимума на луче а(т+1) — а(т). Если процесс оптимизации осуществляется из точки а, лежащей вне допустимой области

(2), то в выражении (4) Цд = <?д(а)— сд, Ма, а член (у/-у<р9) полагается равным нулю, т. е. изменение целевой функции в этом случае не учитывается. При движении из недопустимой области в направлении вектора 5 на каждом шаге по £ контролируется суммарная невязка по ограничениям

^Ма

Если рассматриваемые ограничения линейны относительно проектных переменных, то попадание в допустимую область осуществляется за одну итерацию. В противном случае одной итерации может оказаться недостаточно и как только на каком-либо шаге по Ы суммарная невязка ДСу, начинает увеличиваться, осуществляется переход на следующую итерацию. На каждой последующей итерации определяется новое направление вектора 5<ш), и в зависимости от того находится ли текущая точка в допустимой области или нет, оптимизация идет либо по одной, либо по другой ветви. Итерационный процесс заканчивается при достижении необходимой точности оптимизации. Для вычисления значения градиента каждой из независимых переменных дается малое приращение и вычисляется приращение функции и ограничений с помощью программы прямого расчета аэродинамических характеристик профиля.

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

давления на профиле от заданного /=-"]* {ср— и так далее. Гео-

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

В приведенных примерах данной работы для расчета аэродинамических характеристик трансзвукового обтекания профиля использовалась программа Блынской — Лифшица, основу которой составил метод решения полных уравнений для потенциала течения невязкой жидкости [7].

Минимизация волнового сопротивления. Для иллюстрации работоспособности и эффективности предложенного метода была рассмотрена

задача оптимизации формы профиля, обеспечивающего минимальное волновое сопротивление. В качестве базового был взят один из профилей, с относительной толщиной 12%, рис. 2, а. Минимизация коэффициента волнового сопротивления свх велась на режиме обтекания М<х> = 0,8 при ограничении на относительную толщину профиля с = 0,12. Такой режим обтекания был выбран специально для большей иллюстративности, так как при данном числе М«, базовый профиль при всех значениях коэффициента су имеет достаточно большое волновое сопротивление и процесс оптимизации становится ярко выраженным. С целью исследования влияния коэффициента подъемной силы су на форму и характеристики получаемых оптимальных профилей оптимизация проводилась при различных значениях су = 0,2, 0,4, 0,6. В каждом случае на профиле выбиралось до восьми зон управления.

Результаты оптимизации приведены на рис. 2, а. Анализ представленных зависимостей съх (су) показывает, что полученные оптимальные профили являются остро настроенными. Каждому заданному режиму обтекания (Моо, Су) соответствует свой оптимальный профиль. На рис. 2,6 для одного из оптимальных профилей, а именно профиля, спроектированного на М<х, = 0,8, су=0,4, приведен график зависимости с® (Моо)- Четко прослеживается его острая настройка на режим Моо = 0,8. На том же рисунке для данного профиля представлена иллюстрация сходимости итерационного процесса оптимизации. Она отражает поведение коэффициента с® в зависимости от величины шага улучшающей вариации k на каждой итерации. Видно, что уже на третьей итерации процесс оптимизации близок к завершению. Такая же картина наблюдается и в других рассмотренных примерах.

На рис. 2, в в качестве примера приведено распределение коэффициента давления ср(х) для оптимального профиля, спроектированного на = 0,6. Наблюдаемое «полочное» распределение давления на верхней поверхности характерно и для других полученных оптимальных профилей. Геометрические обводы оптимальных профилей для всех трех значений су = 0,2; 0,4; 0,6 приведены на рис. 2, г. На этом же рисунке представлены таблицы аэродинамических характеристик профилей до и после оптимизации. С ростом коэффициента подъемной силы су максимальная толщина нижней поверхности профиля увеличивается, а верхней, наоборот — уменьшается, причем верхняя поверхность профиля становится более уплощенной. «Подрезка» нижней поверхности профиля в районе хвостика с ростом су увеличивается. Следует отметить, что оптимизация формы профиля в трансзвуковом диапазоне скоростей без учета вязкости на практике не может гарантировать снижения сопротивления. Практическое проектирование необходимо вести, привлекая более совершенную модель обтекания, учитывающую эффекты вязкости.

Во втором примере рассматривается задача определения оптимальной формы профиля минимального волнового сопротивления на режиме Моо = 0,76; су = 0,4 при ограничениях на относительную толщину с = 0,12 и коэффициент момента тангажа —0,175. В качестве базового был

взят один из профилей с относительной толщиной 12%. Коэффициент волнового сопротивления базового профиля, вычисленный на указанном режиме по программе Блынской—Лифшица, был равен 0,0037, коэффициент момента тангажа mz = —0,174. На профиле было выбрано семь зон управления (четыре сверху и три снизу). Проделана всего лишь одна итерация, в результате которой был получен профиль, на котором реализуется практически бесскачковый режим обтекания

с, =0,1 с*гО,Ч сч=0,6

БазоНый профиль сё 0,0197 0,0162 0,0176

mz ■0,3173 -0.3833 -0А512

Оптапамьный профиль ci 0,0025 0,003 0,0081

nii 0,1723 -о,2т -0.37SJ

(с® =0,0005; тг = —0,171). Результаты решения задачи оптимизации представлены на рис. 3, а.

Произошло перераспределение толщины профиля в сторону ее увеличения на верхней поверхности и соответствующего уменьшения на нижней поверхности. Наблюдаемый у разового профиля скачок давления на верхней поверхности в районе л: = 0,2 ликвидировался. Полученный оптимальный профиль имеет достаточно протяженную местную сверхзвуковую зону с почти изэнтропическим сжатием потока. Местные числа Маха невелики и не превышают числа М =1,1.

Сравнение эпюр давления оптимального профиля в вязком и невязком потоках при су = 0,4 представлено на рис. 3,6. Там же представлен и сам профиль. В данной работе результаты расчетов вязкого обтекания профилей выполнены с помощью программы [8], в основу кото-

Рис. 3

рой были положены метод расчета потенциального обтекания [9], с некоторыми его модификациями, и метод [10] для расчета турбулентного пограничного слоя.

Расчет вязкого обтекания оптимального профиля получен для числа Re = 2-107 и координаты точки перехода л:п = 0,06. Влияние вытесняющего действия пограничного слоя более заметно сказывается на распределении давления по верхней поверхности. Однако, как видим, существенной перестройки эпюры давления по сравнению со случаем невязкого обтекания не произошло. Наблюдаемое незначительное перераспределение нагрузки в сторону носовой части профиля привело к уменьшению абсолютного значения коэффициента момента тангажа mz примерно на 9%. Значение коэффициента волнового сопротивления оптимального профиля при обтекании его потоком вязкого газа на указанном режиме равно с\ =0,0008.

ВЯЗ

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

sk

/= í Ms) ~cp3(s))2ds.

sh

Здесь s — параметр, по которому строится сплайн при описании геометрии профиля, sH, Sfc — начальное и конечное значение параметра, а cP3(s), c.P(s)—соответственно заданное и получаемое распределение коэффициента давления на профиле.

Минимизация целевой функции велась на режиме Мос = 0,8; су = 0,4. Заданное распределение срз(х) соответствовало одному из профилей с относительной толщиной 12% при расчете его обтекания на данном режиме. В качестве базового также был взят профиль с относительной толщиной 12%. Его обводы и соответствующее распределение коэффициента давления ср(х) приведены на рис. 4,а. Там же представлено и заданное распределение срз(х).

На профиле выбиралось до десяти зон управления. Было проведено пять итераций. Причем на каждой итерации использовался свой набор параметров, определяющий число и вид каждой из функций формы. Значение целевой функции удалось уменьшить с / = 0,4065 до значения 1 = 0,0077. В результате полностью произошла перестройка эпюры давления в нужную сторону. Получен профиль, обводы которого представлены на рис. 4,6. Здесь же представлено сравнение полученного для данного профиля распределения коэффициента сР(х) и заданного распределения срз(х). Как видим, согласование эпюр вполне приемлемо для того, чтобы закончить процесс оптимизации.

Модификация носовой части. В этом примере разработанным методом заново проектируется передняя часть верхней поверхности профиля NACA 64А109 с целью уменьшения его волнового сопротивления. В практике такие задачи могут возникать довольно часто, когда в силу конструктивных или каких-либо других соображений (например, чтобы не испортить диффузорную часть профиля) разрешается модифицировать только некоторую часть поверхности профиля.

Минимизация коэффициента волнового сопротивления с® велась на режиме Моо = 0,82, а = —0,05°. Модифицировалось только 35% верхней поверхности профиля, при этом относительная толщина профиля оставалась постоянной, с = 0,09. В процессе оптимизации на профиле ис-

Моа-0,8 > С^ОЛ

о

0,1--------заданное распределение срз

------ полученное » Ср

Рис. 4

пользовалось до пяти зон управления. Результаты, полученные после четырех проведенных итераций, представлены на рис. 5, а. Носовая часть у модифицированного профиля более затупленная, что привело к смещению звуковой точки ближе к критической и к более быстрому расширению потока до местной сверхзвуковой скорости. Из сравнения эпюр давлений модифицированного и исходного профилей видно, что интенсивность скачка заметно снизилась. Коэффициент волнового сопротивления профиля стал равным с\ =0,0029 по сравнению с с®= = 0,0068 для исходного профиля. За счет перераспределения нагрузки вдоль хорды уменьшилось абсолютное значение коэффициента момента тг с —0,107 до —0,093. При этом коэффициент су увеличился с 0,194 до 0,22.

Эпюры давления для исходного и модифицированного профилей, а также зависимость коэффициента сопротивления сх(М«,), рассчитанные

а)

1,0 х

-профиль МАСАБЧА109

0,5-г----подисрицироданный профиль

данные из работы [¿]

-0,5-------------------1,0 х

Мо* = 0,82, а. = -0,05°, Ле=210 7; х^ООб

-0.5

Ю

1,0х

■ профиль НАСАВШОЭ модифицированный профиль

Рис. 5

с учетом вязкости при 1?е = 2-107, л:п = 0,06, представлены на рис. 5,6. Проведенная модификация профиля увеличила число М» резкого нарастания сопротивления.

Отметим, что последняя из рассмотренных здесь задач решалась и методом Хикса в работе [2]. Проведенное сопоставление результатов показало, что уровни понижения волнового сопротивления в обоих случаях примерно одинаковы (Дс®~ 0,004). Качественно схожи также формы профилей и эпюры давлений, рис. 5, а.

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

1. Hicks R. М., Vanderplaats G. N. Design of low-speed airfoils by numerical optimization. — SAE Paper 750524, 1975.

■2.'Hicks R. М., Vanderplaats G. N., Mur man E. М., К i n g R. R. Airfoil section drag reduction at transonic speeds by numerical optimization. — SAE Paper 760477, 1976.

3. H i с к s R. М., Vanderplaats G. N. Aplication of numerical optimization to the design of supercritical airfoils without drag-creep. — SAE Paper 770440, 1977.

4. Ляпунов С. В. Построение профилей минимального волнового сопротивления с различными ограничениями. — Ученые записки ЦАГИ, 1986, т. 17, № 4.

5. Б у з о в е р я Н. П. Об аппроксимации таблично заданной функции кубическим сплайном со сглаживанием. — Труды ЦАГИ, 1983, вып. 2181.

6. Rosen J. В. The gradient projection method for non-linear programming, Part II, Non-linear constraints. — J. SIAM 9, 1961.

7. Л ифшиц Ю. Б. К теории трансзвуковых течений около профиля. — Ученые записки ЦАГИ, 1973, т. 4, № 5.

8. Б у з о в е р я Н. П., Кроткое Д. П., Маркин В. С. Некоторые особенности обтекания модели профиля в трансзвуковой аэродинами-

ческой трубе с проницаемыми стенками. — Труды ЦАГИ, 1985, вып. 2298.

9. С а г 1 s о n L. A. Transonic airfoil flow field analysis using cartesian coordinates. — NASA CR-2577, August 1975.

10. Nash G. F., Macdonald A. G. J. The calculation of momentum thickness in. a turbulent boundary layer at Mach number up to unity.—

Aero. Res. Concil C. P. 1967, iN 963.

Рукопись поступила 301IV 1987