О геометрическом генотипе форм аэродинамических профилей Текст научной статьи по специальности «Математика»

Том XLV

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

2014

№ 3

УДК 629.735.33.015.3.025.73 533.6.011.35:629.7.025.73

О ГЕОМЕТРИЧЕСКОМ ГЕНОТИПЕ ФОРМ АЭРОДИНАМИЧЕСКИХ ПРОФИЛЕЙ

А. А. НИКОЛЬСКИЙ

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

Ключевые слова: параметризация, профиль, сплайн, оптимизация.

В общем случае форма контура профиля y (x) задается таблицей координат верхнего контура (x, yui) и нижнего контура (x, yu ), i = 1,..., N, где N — количество точек на контуре

профиля. Совокупность из 2N ординат, определяющих форму профиля, можно считать его геометрическим генотипом. Для задач аэродинамического проектирования необходимо аппроксимировать форму контура с помощью минимального количества параметров, охватывая при этом максимально широкий класс контуров. Однако известные методы параметризации [1—9] сводятся либо к выделению определенного ограниченного класса контуров, представимых например аналитическими кривыми, сплайнами Безье, полиномами Бернштейна, конечными рядами Фурье или сводятся к локальным вариациям заданного контура и не охватывают всего множества возможных контуров профилей.

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

1. Отнесем все контуры аэродинамических профилей к трем основным типам по следующим признакам:

Профили «обыкновенного» типа, а именно: ординаты верхнего контура yu (x) монотонно возрастают до точки максимума yu max, после чего монотонно убывают до задней кромки, а ординаты нижнего контура yi (x) монотонно убывают до точки минимума yi min, после

чего монотонно возрастают до задней кромки.

Трансзвуковые профили с подрезкой у задней кромки нижнего контура и профили лопаток турбин геометрически сходные в том, что монотонность функции yi (x) может нарушаться после точки минимума yi min.

Некоторые вертолетные профили с подрезкой у задней кромки верхнего контура характерные тем, что монотонность функции yu (x)

может нарушаться после точки максимума yu max.

В действительности профили первого типа можно считать частным случаем профилей второго типа и объединить их в общий тип A.

НИКОЛЬСКИМ Александр Александрович

кандидат физико-математических наук, начальник сектора ЦАГИ

Рис. 1. Профили основных типов

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

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

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

2. Без ограничения общности рассмотрим далее профили типа А. Функции, описывающие контуры профилей, представим в виде у, (х) = у,тах/, (х) и уи (х) = уитах/ (х). В качестве эталонной функции, множество значений которой упорядочено так же (рис. 2), как множества значений функций / и / , удобно использовать следующую функцию:

/ = (1 -4х). (1)

Пусть значения функций /г и /и заданы в точках х{, г = 1,..., N. Приравнивая /е (х« ) = / (хI), из (1) определим значения ха и тем самым дискретную зависимость хй (хг). Повторяя то же для функции /и, находим зависимость х^ (хг). Таким образом, набор генов, состоящий из двух последовательностей хи и хш, г = 1,..., N, соответствует некоторому семейству профилей типа А. Добавляя два дополнительных гена уитах и у,тах, получим конкретный экземпляр профиля данного семейства. Обратно, зная эталонную функцию (1), с помощью зависимостей хп (хг) и хиг (хг) находим функции / и /и , определяющие семейство профилей

Рис. 2. Функции, описывающие контуры профиля и эталонная функция

Рис. 3. Зависимости xt (x) и xu (x), образующие семейства профилей RAE2822

и VR

типа A. Все выше сказанное справедливо и для профилей типа B (вместо зависимостей f и fu следует взять зависимости ft и f). На рис. 3 показаны зависимости xt (x) и xu (x) для семейств

профилей RAE2822 и VR [10]. Отметим, что эти функции принадлежат довольно ограниченному классу функций, монотонно возрастающих от нуля до 1.

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

Поскольку мы имеем дело с функциями, заданными таблично, для их аппроксимации будем использовать сплайны. Однако обычный кубический интерполяционный сплайн для аппроксимации контура профиля y (x) не подходит из-за бесконечной производной — (0) в носке про-

dx

филя. По этой причине предлагается использовать альтернативный сплайн, который назовем сплайном контура профиля, сокращенно S(a) (airfoil spline).

3. Снова рассмотрим верхний контур профиля. Сделаем преобразование координат b ="Jx.

Пусть rc — радиус кривизны в носке верхнего контура, тогда dy (0) = <J2rc . Теперь во вспомогательной плоскости by контур профиля можно аппроксимировать обычным интерполяционным кубическим сплайном:

sfa)(b)=a^ + ьг (b-bi)+с (b-b )2+di (b-b )3.

В физической плоскости получим альтернативный сплайн:

sfa)(x) = ai + b (4x -<Jx~) + ci ((yfx -Jx~)) + di (4x -Jx~) .

Зная сплайн S(a), легко вычислить необходимые для многих приложений величины: длину дуги, углы наклона и кривизну контура. При необходимости можно сгладить контур профиля, построив в плоскости by вместо интерполяционного сглаживающий кубический сплайн [11].

Теперь определим, какое минимальное количество точек с достаточной точностью аппроксимирует контур профиля. Последовательность из m точек xi, i = 1,...,m далее для краткости будем называть разбиением и обозначать Pm. Точность аппроксимации контура профиля оценим на основе отклонения ординат аппроксимированного с помощью сплайна контура от ординат

исходного контура: 55 = £(а)(х) — уи (х). В качестве условия необходимой точности примем 55 <5№/, где — точность изготовления модели крыла для испытаний в АДТ. По величине

среднеквадратичной разности коэффициентов давления с =

1 ы /

N ^ (

2=1

с(№) — с ( тЯ

С Р2 Р2 I

по исходному

(ср) ( х)) и аппроксимирующему (сРт) ( х)) контурам будем определять достаточность точности.

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

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

Дальнейший анализ проводится на примере расчетных данных, полученных для вертолетного профиля УЮ2. На рис. 5 показано 55Г =55 -104 — приведенное отклонение ординат верхнего и нижнего контуров профилей вдоль хорды профиля для разбиения Р от исходных орди-

0.08

0.06

0.04

0.02

-0.02

-0.04

м + *• * 1 в + * < 4 а * * 5 +Р29 3 аРю 1- +

^ ез (Н

0Л 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 и* I

_ н в * а 1 * + а ^ а * < - 1 В з 1

Рис. 4. Расположение контрольных точек (разбиения) по контурам профиля

Рис. 5. Приведенное отклонение ординат верхнего и нижнего контуров профиля УР (разбиение Р29) от ординат исходных контуров (разбиение Р41)

Рис. 6. Кривизны верхнего и нижнего исходного и интерполированного контуров

нат контуров профиля (разбиение Р41). Поскольку для трубных моделей

видно, что на большей части контура 5Л <С 5"''. Некоторый всплеск зависимости в хвостовой части объясняется сглаживающим эффектом при уменьшении числа точек, определяющих контур. Этот факт иллюстрирует рис. 6, где приведены кривизны к верхнего и нижнего исходного и интерполированных контуров.

Для проверки выполнения второго критерия были проведены расчеты распределений давления по исходному и интерполируемым контурам профиля УЯ12. На рис. 7 показаны распределения коэффициента давления на характерном для вертолетных профилей режиме = 0.8, су = 0. Здесь М — число Маха набегающего потока, су — коэффициент подъемной силы профиля. Видно идеальное совпадение эпюр давления, за исключением малой области в районе задней кромки, где сглаженность интерполированного профиля приводит к более гладкой эпюре.

Такое практически глобальное совпадение ^с» 10"31 наблюдается и на других режимах и на всем массиве рассмотренных профилей. Поэтому с достаточным основанием можно считать, что сплайн профиля £ ^а ^ с универсальным разбиением Р29 с необходимой и достаточной точностью аппроксимирует, вообще говоря, произвольный контур профиля. Для краткости далее такой

Рис. 7. Распределения коэффициента давления по исходному и интерполированному контурам (Мш = 0.8, су = 0)

Рис. 8. Приведенное отклонение ординат верхнего и нижнего контуров сглаженного профиля от ординат базового контура

Рис. 9. Кривизны верхнего и нижнего базового и сглаженного базового контуров

сплайн будем называть базовым и обозначать (basic airfoil spline). Применение такого

сплайна в задаче локальной оптимизации можно найти в [12].

Для сглаживания контуров профиля вместо базового сплайна Sb ) можно использовать

сглаживающий базовый сплайн S^m. Напомним, что все сплайны строятся в плоскости by, где

производная — (0) в носке профиля конечна. Выберем параметры сглаживания таким образом, d b

чтобы приведенное отклонение находилось в диапазоне |5sr\ < 1, как это видно на рис. 8. При этом с запасом удовлетворим принятому ранее условию необходимой точности. Кривизны верхнего и нижнего базового (S(a)) и сглаженного базового (S(aконтуров приведены на рис. 9.

Такое сглаживание не ухудшает совпадение расчетных эпюр давления (с «10_3 ). Поэтому в качестве базового сплайна оправдано использование сплайна S(aJm.

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

I=

t л( cpi - cpn )7n

Здесь cpt — целевое распределение давления, n — число расчетных точек на контуре профиля. В качестве переменных проектирования используем установленный набор генов: (xti (xi) и xui (x¡), i = 1,..., m). Поскольку носик и задняя кромка профиля фиксированы, число переменных проектирования mdv = 2(m — 2) + 2 (два дополнительных гена регулируют толщину профиля). В действительности эффективность оптимизационной процедуры можно повысить, если попеременно использовать два подпространства переменных проектирования: xTi (xi) и xUt (xi), при этом mciv = m. Вторым существенным фактором ускорения процедуры является решение оптимизационной задачи на частом разбиении, используя в качестве начального приближения решение на редком разбиении.

На основе сопоставлений результатов решения обратных задач для массивов профилей и разбиений были установлены итоговое минимальное разбиение Íjq и разбиение P15, практически идентичное базовому P29 по рассчитанным распределениям коэффициента давления. При решении обратных задач для разбиения P10 в качестве начального приближения принимался контур профиля NACA 23012.

Дальнейший сравнительный анализ проведен в окрестности трех характерных для вертолетных профилей режимах: = 0.4, cy = 1.5 (режим cymax); = 0.6, cy = 0.6 (режим

Kmax); MOT = 0.8, cy = 0 (режим Mdd). Здесь Mdd — критическое число Маха набегающего

потока, cy max — коэффициент максимальной подъемной силы профиля, Kmax — максимальное

аэродинамическое качество.

Приведены данные для профилей, построенных на основе сплайна S(aJm (разбиений P15 и базового P29), и профиля, полученного из решения обратной задачи (разбиения P10), обозначенного на рисунках, как P10 (inv). На рис. 10, 11 видно, что распределения коэффициента давления по контуру профиля разбиения P15 на всех режимах практически совпадают с распре-

Рис. 10. Распределение коэффициента давления по базовому и интерполированному контурам (М<„ = 0.4, cy= 1.5 (слева) и = 0.6, cy = 0.6 (справа))

Рис. 11. Распределения коэффициента давления по базовому, интерполированному и полученному из решения обратной задачи контурам (Мш = 0.8, су = 0)

делением коэффициента давления базового профиля. Распределения коэффициента давления по контуру профиля разбиения Р10 на режимах сутах и Ктах также близки к распределениям

коэффициента давления базового профиля. Однако на околозвуковом режиме (см. рис. 11) распределение давления по контуру профиля Рш (ту), полученному из решения обратной задачи,

лучше согласуется с базовым, чем распределение коэффициента давления для интерполированного профиля Рш. Распределения коэффициента давления по контуру Р^ (ту) не приводятся,

поскольку они практически идентичны распределениям давления по контуру Р15.

Для того чтобы оценить аэродинамические характеристики исследуемых профилей, были проведены расчеты методом СЕБ ЯЛЫБ [13]. На рис. 12, 13 видно, что на характерных режимах и величины сутах и Ктах, и зависимости су (а) и су (сх) практически совпали (Лсутах < 0.01,

ЛКтах < 0.1). На околозвуковых режимах (см. рис. 13), также наблюдается лишь незначительное отличие зависимостей ст0 (Мда) и сх ((Лсх » 0.0003, Лст » 0.0015), а величины М практически совпали. Подобные результаты (качественные и количественные) получены на всем

Рис. 12. Зависимости су (а), = 0.4 (слева) и су (сх), = 0.6 (справа) (расчет СЕБ КАКБ)

Рис. 13. Зависимости cm0 (Мш) (слева) и cx (Мш) (справа), = 0.8, cy = 0 (расчет CFD RANS)

Рис. 14. Зависимость cy (cx), = 0.76, профиль RAE2822 (расчет CFD RANS)

рассмотренном массиве профилей. В частности, на рис. 14 показана зависимость cy (cx) при

= 0.76 для трансзвукового профиля RAE2822. В этом случае AK = -0.5 — 0.5 для варианта

P10 (inv) и AK < 0.1 для варианта P15.

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

аппроксимировать сплайном S(a), построенным по 10 точкам. Такой сплайн предлагается называть аэродинамическим сплайном профиля и обозначать S(a) (aerodynamic airfoil spline). На основе сплайна S(a) аэродинамические характеристики многорежимных вертолетных профилей

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

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

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

Дополнительно можно заметить, что сплайны S(а) можно использовать для построения 3-d каркасной геометрии в CAD, для построения и быстрого перестроения расчетных сеток. Эти сплайны можно также использовать для формирования универсальной (по количеству и расположению точек по хорде) геометрической базы данных профилей. При этом любые геометрические характеристики профилей могут быть определены аналитически.

В оптимизационных задачах зависимости xt (x) и xu (x), образующие семейства профилей

типа A, или зависимости xt (x) и xt (x), образующие семейства профилей типа B, позволяют

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

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

ЛИТЕРАТУРА

1. Sobieczky H. Parametric airfoils and wings // Notes on Numerical Fluid Mechanics. 1998. V. 68, p. 71 — 88.

2. Samareh J. A. Survey of shape parameterization techniques for high-fidelity multi-disciplinary shape optimization // AIAA J. 2001. V. 39, N 5.

3. Robinson G. M. and Keane A. J. Concise Orthogonal Representation of Supercritical Airfoils // J. of Aircraft. 2001. V. 38, N 3.

4. Song W. and Keane A. A study of shape parameterisation airfoil optimisation // AIAA Paper 2004-4482.

5. Padula S. and Li W. Options for robust airfoil optimization under uncertainty // 9th AIAA Multidisciplinary Analysis and Optimization Symposium, 2002.

6. H i c k s R. M. and H e n n e P. A. Wing design by numerical optimization // J. of Aircraft. 1978. V. 15, p. 407—412.

7. S a m a r e h J. A. Aerodynamic shape optimization based on free-form deformation // AIAA Paper 2004-4630.

8. Brenda M. Kulfan and John E. Bussoletti. «Fundamental» parametric geometry representations for aircraft component shapes // AIAA Paper 2006-6948, 11th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2006.

9. N i k o l s k y A. A. Some aspects of helicopter airfoil design // Twenty first European Ro-torcraft Forum. 1995. V. 2, N 17.

10. Dadone L. U. Advanced airfoils for helicopters rotor application // US Patent N 4, 341, 795, 1982.

11. Reinch C. H. Smoothing by spline functions // Numerische Mathematik 1967. 10, p. 177 —183.

12. Никольский А. А. Оптимизация передних кромок вертолетных профилей // Ученые записки ЦАГИ. 2008. Т. XXXIX, № 4.

13. Morrison J. H. A compressible Navier — Stokes solver with two-equation and Reynolds stress turbulence closure models // NASA CR-4440, May 1992.

Рукопись поступила 13/II2013 г.