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

2015 Математика и механика № 1(33)

УДК 533, 629.7 DOI 10.17223/19988621/33/7

С.А. Орлов, С.В. Пейгин, К.А. Степанов, С.В. Тимченко

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

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

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

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

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

1 Работа выполнена при поддержке ФЦП «Кадры», соглашение .№ 14.B37.21.0733.

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

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

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

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

• Аэродинамические ограничения: заданный постоянный коэффициент подъёмной силы Сь и максимально допустимый момент тангажа См.

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

- относительная толщина секции крыла(//с),- ;

- радиус кривизны передней кромки секции крыла (Яь)г;

- угол задней кромки секции крыла (дт),;

- локальные толщины секции крыла (уЛ)р .

В этих ограничениях / = 1,...,Л^ - число секций по размаху крыла, а ] = 1,., ЛЦ/) - число ограничений на локальную толщину в секции номер /. Общее число ограничений ЛСц на крыло - 20-30.

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

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

• точный расчет сопротивления очень сложен для реальных конфигураций;

• нет общего решения проблемы глобального геометрического представления аэродинамических поверхностей;

• оптимальный поиск происходит в пространстве высокой размерности;

• необходим эффективный учет большого числа нелинейных ограничений;

• решение задачи требует огромнейшего объема вычислений.

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

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

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

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

Поэтому в качестве целевой функции в работе используется следующий функционал:

0,1 + [(//с)* - (/ /с)], если (//с)* > (Г /с);

0,2 + [Л* -Яь],

Q =

0,3 + [0* -0T] 0,5,

CD

если если если если

rl > rl ;

/ (t) < У (t);

где у" (/), у (/) - соответственно контуры верхней и нижней поверхностей крыла, задаваемые при помощи линий Безье, опорные точки которых собственно и необходимо найти в процессе оптимизации.

1. Метод решения

Уравнения Навье - Стокса вязкой сжимаемой жидкости могут быть записаны в следующей форме:

qt + div C = div V, (1)

где тензор C = f g, h) содержит конвективные члены, тензор V = (r, s, t) содержит вязкие члены, q = (р, pu, pv, pw, E)T , p - плотность, (u, v, w) - вектор скорости, E -энергия, t - время, f g, h - невязкие (конвективные) потоки и r, s, t - вязкие потоки.

Компоненты невязких потоков имеют вид

f(q) = uq + p(0, 1, 0, 0, u)T,

g(q) = vq + p(0,0, 1, 0, v)T, h(q) = wq + p(0, 0, 0, 1,w)T.

Компоненты вязких потоков имеют вид

r(q) = Д(0, Tu, T21, Тз1, 0:)T

S(q) = Д(0, T12, T22, T32, C2)T

t(q) = д(0, T13, T23, T33, C3)T, где составляющие тензора вязких напряжений задаются следующим образом:

Tn = (4/3)u - (2/3) vy - (2/3)wz,

T21 = T12 = uy + vx,

Т31 = Т13 = Uz + wx, Т22 = (4/3) Vy - (2/3)Ux - (2/3) Wz,

T32 = T23 = Uz + Wy, T33 = (4/3)wz - (2/3)«x - (2/3) Vy,

O = u T11 + v T12 + W T13 + (c2)x / ((y - 1)Pr),

02 = u T21 + V T22 + W T23 + (c2)y / ((y - 1)Pr),

03 = u T31 + V T32 + W T33 + (c2)z / ((y - 1)Pr).

В этих формулах ц - значение вязкости, y - отношение теплоёмкостей газа, Pr -число Прандтля,

p = (y - 1)[E - 0,5р(м2 + v2 + w2], c2 = YP/P, H = (E + p)/p.

2.1. Численный метод

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

Рассмотрим структурированную расчётную сетку, состоящую из ячеек, которые образуют структуру типа (i, j, k). Интегрируя по каждой ячейке, мы получаем систему обыкновенных дифференциальных уравнений, к которой применяется процедура интегрирования по времени.

Запишем аппроксимацию для ячейки с индексами i, j, k:

(Qy^yA + [C ■ (Sn)]i+0,5j,k - [C ■ (Sw)] i-0,5j,k +

+ [C ■ (Sw)]i,j+0:5,k - [C ■ (Sn)]Ii;-0,5,k +

+ [C ■ (Sw)]ij,k+0,5 - [C ■ (Sw)]i,j,k-0,5 = (2)

= [V ■ (Sn)]1+0,5j',k - [V ■ (Sn)]i-0,5j',k

+ [V ■ (Sn)],,j+0..5,k - [V ■ (Sw)]ij-0,5,k

+ [V ■ (Sw)] ij ,k+0,5 - [V ■ (Sw)] ij,k—0,5 , где Qj - объём ячейки, qy,k - среднее значение по ячейке и S - площадь грани ячейки. Дробные индексы указывают, с какой стороны ячейки берётся поток в квадратных скобках.

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

В пространственной аппроксимации потоков конвективные потоки на гранях ячеек интерполируются по данным в центрах ячеек посредством двух интерполяционных операторов: характеристического оператора первого порядка и ENO (Essentially non-Oscillatory - Существенно не-Осцилляционного) оператора высокого порядка. Шаблон характеристической схемы первого порядка, фактически применяемой при релаксации, состоит из одной точки, выбранной по знаку соответствующего собственного числа.

Метод ENO (предложенный А. Хартеном и С. Ошером и затем упрощённый С.-В. Шу и С. Ошером, [4, 5] применяется посредством выбора интерполяционного шаблона по локальным характеристикам и гладкости потоков и может изме-

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

Мы также используем метод поправки на дефект, в котором целевая дискретизация уравнений на самой тонкой сетке отличается от дискретизации, применяемой при многосеточной релаксации. Для релаксации используется аппроксимаци-онный оператор первого порядка точности, а поправка на дефект вычисляется посредством аппроксимации £N0 высокого порядка.

Таким образом, мы применяем схему £N0 только при вычислении поправки на дефект, а большая часть работы осуществляется посредством относительно дешёвой в вычислительном отношении характеристической схемы первого порядка. Вязкие члены аппроксимируются напрямую. Для интегрирования по времени используется схема Рунге-Кутта третьего порядка, сохраняющая полную вариацию (ТУБ), разработанная С.-В. Шу и С. Ошером (см. [5]). Более подробное описание численного метода дано в [6] и [7].

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

В этом разделе приводим результаты одноточечной и многоточечной оптимизации крыла самолета Богтег-728 [8]. Главная точка оптимизации была выбрана при Сь = 0,5 и М = 0,78. Дополнительные точки оптимизации были выбраны при М = 0,8 (С = 0,5) и при М = 0,2 (С = 1,19). Геометрические ограничения были наложены на толщину крыла, радиус кривизны передней кромки секции крыла и угол задней кромки секции крыла. Дополнительное (аэродинамическое) ограничение накладывалось на момент тангажа.

Оптимизационные условия и ограничения приведены в таблице. В общей сложности были выполнены три оптимизации: две одноточечные оптимизации и одна трехточечная. Число секций крыла было равно трем - корневая, средняя и концевая секции. Одноточечные оптимизации отличаются величиной С М. В первой из оптимизаций момент подъемной силы был свободным, в то время как вторая оптимизация (а также многоточечная) сохраняли момент подъемной силы на уровне оригинального крыла. В дальнейшем эти три варианта будут обозначаться соответственно СазеБ728_1, СазеБ728_2 и СазеБ728_3.

Крыло Богше-728, ограничения и условия оптимизации

с; М С* ^м (Г/С)! (Г/С)2 (Г/С)3

Са8еБ728 1

0,500 0,78 1,0 -0,130 0,142 0,119 0,104

Са8еБ728 2

0,500 0,78 1,0 -0,072 0,142 0,119 0,104

Са8еБ728 3

0,500 0,78 0,75 -0,072 0,142 0,119 0,104

0,500 0,80 0,23 -0,082 0,142 0,119 0,104

1,190 0,20 0,02 -0,012 0,142 0,119 0,104

Главная точка дизайна лежит в высоком трансзвуковом диапазоне, где в потоке около оригинального крыла развивается сильный скачек. Одноточечная оптимизация (Са5еБ728_1) устраняет этот скачек. Соответствующие продольные распределения давления приведены на рис. 1 и 2, соответственно.

Cp

-0,8

0,8

1,6

д ü Д Д Л

& & ó,

li

ь

0

0,2

0,4

0,6

0,8

Z/C

0,2

X/C

-0,2

-0,4

Рис. 1. Оригинальное крыло Dornier-728, M = 0,78; CL = 0,5. Распределение давления по верхней и нижней поверхностям крыла

0

0

Cp

-0,8

0,8

1,6

0

0,2

0,4

0,6

0,8

Z/C

0,2

X/C

0,2

0,4

Рис. 2. Оптимизированное крыло Богтег-728 (Са8еБ728_1), М = 0,78; С = 0,5. Распределение давления по верхней и нижней поверхностям крыла

0

В условиях такого проектирования, было достигнуто снижение на 32,6 каунта (от начальных 213,7 каунта). Минимально возможное значение сопротивления при этом значении коэффициента подъемной силы составляет около 180 каунтов, что практически идентично сопротивлению оптимизированного крыла.

В этом оптимизационном случае достигалось значение момента тангажа, достаточно далекое от исходного крыла. Интересно оценить влияние этого ограничения на результаты оптимизации. Более строгое ограничение на момента тангажа (Caí = -0,07), использованное в варианте CaseD728_2, выливалось в снижение

уменьшения сопротивления: в этом случае полное сопротивление оптимизированного крыла равно 189,7 каунтов по сравнению с 181,1 каунтами в случае Са8еБ728_1.

Характеристики нештатной оптимальной формы Са5еБ728_1 вполне удовлетворительны с точки зрения подъемной силы и сопротивления, но она имеет слишком высокое значение отрицательного момента тангажа. С другой стороны, профиль Са5еБ728_2 имеет подходящее значение момента тангажа в связи с введением соответствующих ограничений, но его взлетные характеристики ухудшились по сравнению с оригинальной крылом. Для обеспечения как хорошего качества проектирования, так и хороших характеристик крыла вне главной точки проектирования проведена многоточечная оптимизация (Са5еБ728_3).

Результаты проведенной оптимизации обладают следующими характеристиками: в главной точке оптимизации (Сь = 0,5 и М = 0,78) общее сопротивление оптимизированного крыла составило 192,8 каунтов (значение, близкое к Са5еБ728_2). Соответствующее поперечное распределение давления в средней секции крыла показано на рис. 3.

Ср

-0,8

0,8

1,6

д А

й Л Д

— Л

г I? 17 V

(

V

НС

0,2

-0,2

-0,4

0

0,2

0,4

0,6

0,8

XIС

Рис. 3. Оптимизированное крыло Богтег-728 (Са8еБ728_3), М = 0,78; Сь = 0,5. Распределение давления по верхней и нижней поверхностям крыла

0

0

Для вторичной точки оптимизации при большом значении числа Маха сопротивление составляет 212,6 каунта (по сравнению с 244,1 каунта у начального крыла и 216,7 каунта у формы Са5еБ728_2). При этом интенсивность первоначального скачка уменьшается за счет многоточечной оптимизации.

В целом характеристики оптимизированных крыльев в сравнении с исходным представлены на рис. 4 - 6. Поляры (подъемная сила)/(сила сопротивления) при М = 0,78 (рис. 4) и при М = 0,80 (рис. 5) показывают, что оптимизация позволила улучшить аэродинамические характеристики при высоких коэффициентах подъемной силы (С > 0,4) без снижения производительности при более низких Сь.

Что касается зависимости роста коэффициента сопротивления от числа Маха в крейсерском режиме Сь = 0,5 (рис. 6), применение многоточечной оптимизации привело к сдвигу точки кризиса волнового сопротивления в сторону больших значений числа Маха (от оригинального значения М = 0,76 до 0,80).

Сь 0,5 0,4 0,3 0,2

0,1

................................ .................................

1Л;ГСГЧ1ЬК ГЛИ \fVII4ld - СаБе 13723 1 ..................

СагеР72В 3 ........

50

100

150

200

Со

Рис. 4. Поляры (подъемная сила)/(сила сопротивления) при M = 0,78

Сь

0,1

................................ ..................

У

у □ОЙШЕР? 728 \ЛЛККЗ --

/ Са£е0728_3 ........

50

100

150

200

Со

Рис. 5. Поляры (подъемная сила)/(сила сопротивления) при M = 0,80

Со

300

250

200

150

ОСЖЫ1ЕР 728 -Сазе0728_3 - -------

//

/

0,6

0,65

0,7

0,75

0,8

М

Рис. 6. Зависимость роста коэффициента сопротивления от числа Маха

Ожидалось, что включение точки взлета в многоточечную оптимизацию должно, по крайней мере, сохранить крылу высокую подъемную силу при малых числах Маха. Действительно, анализ результатов подтвердили это предположение. В самом деле, подъем кривой оптимизированного крыла немного лучше, чем у оригинала.

ЛИТЕРАТУРА

1. Holland J.H. Adaptation in natural and artificial systems. Ann Arbor: The University of Michigan Press, 1975.

2. Michalewicz Z. Genetic algorithms + data structures = evolution programs. New York: Springer-Verlag, 1992. Artificial Intelligence.

3. Пейгин С.В., Periaux J., Тимченко С.В. Применение генетических алгоритмов для оптимизации формы тела по тепловому потоку // Математическое моделирование. 1998. Т. 10. № 9. C. 111-122.

4. Harten A., Osher S. Uniformly high-order accurate non-oscillatory schemes // I. SIAM Journal of Numerical Analysis. 1987. V. 24. P. 279.

5. Shu C.-W. and Osher S. Efficient implementation of essentially non-oscillatory shock capturing schemes 1 // Journal of Computational Physics. 1989. V. 83. No. 1. 1989. P. 3

6. Epstein B. Averbuch A. and Yavneh I. An accurate ENO driven multigrid method applied to 3D turbulent transonic flows // Journal of Computational Physics. 2001. V. 168. P. 316-328.

7. Epstein B., Peigin S.V. Implementation of WENO (Weighted Essentially Non-oscillatory) Approach to Navier-Stokes Computations // International Journal of CFD. 2004. V. 18. No. 3.

8. Казаков В.Ю. , Пейгин С.В., Тимченко С.В. Оптимизация траектории входа в атмосферу земли по интегральному тепловому потоку // ЖПМТФ. 2000. Т. 41. № 4. С. 112-121.

Статья поступила 20.01.2015 г.

Orlov S.A., Peigin S.V., Stepanov K.A., Timchenko S.V. EFFECTIVE IMPLEMENTATION OF NONLINEAR CONSTRAINTS IN OPTIMIZATION OF THREE-DIMENSIONAL TRANSONIC WINGS. DOI 10.17223/19988621/33/7

A new approach to effective implementation of non-linear constraints for optimization problems solving using genetic algorithms (GA) is proposed. The feature of the approach is to change the traditional strategy in which the search path can only pass through permissible (satisfying the constraints) point by admitting routes passing through both valid and invalid points. The basic idea of this approach is that the information from the "forbidden" (i.e., not satisfying the constraints) areas can be very important, and the path to the optimal point that runs through these areas can appear to be much shorter. The method was applied to the problem of multiobjective optimization of aerodynamic shapes depending on different geometrical and aerodynamic constraints. The results showed that the method maintains high reliability with preservation of the traditional GA computational costs (in the calculation of the objective function on the basis of the full Navier-Stokes equations) at an acceptable level.

Keywords: Conditional optimization, genetic algorithms, nonlinear constraints, Navier - Stokes equations, three-dimensional transonic wings.

ORLOV Sergey Aleksandrovich (Candidate of physical and mathematical science, Assoc. prof of Tomsk State University, Tomsk, Russian Federation)

PEIGIN Sergey Vladimirovich (Doctor of physical and mathematical science, Professor, Tomsk State University, Tomsk, Russian Federation)

STEPANOV Kiril Aleksandrovich (Postgraduate, Tomsk State University, Tomsk, Russian Federation)

TIMCHENKO Sergey Viktorovich (Doctor of physical and mathematical science, Tomsk State University, Tomsk, Russian Federation)

E-mail: tsv@ftf.tsu.ru

REFERENCES

1. Holland J.H. Adaptation in natural and artificial systems. Ann Arbor, The University of Michigan Press, 1975.

2. Michalewicz Z. Genetic algorithms + data structures = evolution programs. New York, Springer-Verlag, 1992. Artificial Intelligence.

3. Peygin S.V., Periaux J., Timchenko S.V. Primenenie geneticheskikh algoritmov dlya optimizatsii formy tela po teplovomu potoku. Matematicheskoe modelirovanie, 1998, vol. 10, no. 9, pp. 111-122. (in Russian)

4. Harten A., Osher S. Uniformly high-order accurate non-oscillatory schemes. I. SIAM Journal of Numerical Analysis, 1987, vol. 24, pp. 279.

5. Shu C.-W. and Osher S. Efficient implementation of essentially non-oscillatory shock capturing schemes 1. Journal of Computational Physics, 1989, vol. 83, no. 1, pp. 3

6. Epstein B. Averbuch A. and Yavneh I. An accurate ENO driven multigrid method applied to 3D turbulent transonic flows. Journal of Computational Physics, 2001, vol. 168, pp. 316-328.

7. Epstein B., Peigin S.V. Implementation of WENO (Weighted Essentially Non-oscillatory) Approach to Navier-Stokes Computations. International Journal of CFD, 2004, vol. 18, no. 3.

8. Kazakov V.Yu. , Peygin S.V., Timchenko S.V. Optimizatsiya traektorii vkhoda v atmosferu zemli po integral'nomu teplovomu potoku. Prikladnaya mekhanika i tekhnicheskaya fizika, 2000, vol. 41, no. 4, pp. 112-121. (in Russian)