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

УДК 533.6

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

© 2013 г. А.А. Бондарчук, К.И. Мещеряков, М.А. Сумбатян

Бондарчук Алексей Алексеевич - кандидат физико-математических наук, старший преподаватель, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, e-mail: melchior@list.ru.

Мещеряков Константин Игоревич - аспирант, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, e-mail: m.keyran@gmail.com.

Сумбатян Межлум Альбертович - доктор физико-математических наук, профессор, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, e-mail: sumbatma@mail.ru.

Bondarchuk Aleksey Alekseevich - Candidate of Physical and Mathematical Science, Senior Lecturer, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, Russia, 344090, e-mail: melchior@list.ru.

Mescheryakov Konstantin Igorevich - Post-Graduate Student, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: m.keyran@gmail.com. Sumbatyan Mezhlum Albertovich - Doctor of Physical and Mathematical Science, Professor, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia, email: sumbatma@mail.ru.

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

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

The paper presents a method of finding the optimal geometry of the rotating blade of a wind turbine. Every normal section is a plane airfoil of the same class. The sections may differ from each other by their sizes and by angles of incidence. The attack angle for every normal section is defined by the incidence angle as well as by the angle which depends upon the wind speed. The problem is solved with the following assumptions: the fluid is nonviscous; at every section the problem is plane i.e. 2D. Under these hypotheses there is constructed a genetic algorithm which can optimize the blade. Some results of calculations are demonstrated.

Keywords: aerodynamics, wind energy, optimization, genetic algorithms, blade, profile.

Постановка задачи и основные гипотезы

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

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

этом ось вращения параллельна оси х .

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

Однородный набегающий на лопасти поток

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

вниз и равна юг, таким образом, относительная скорость потока образует угол р с плоскостью хг . При этом угол атаки а ад меньше по сравнению со случаем, когда вращательное движение отсутствует. Угол установки зависит от угла закрутки у(г),

который определяет, насколько закручено сечение на расстоянии г от оси вращения по отношению к углу установки вблизи оси:

2 2 ч2 п

Vfo/ =У2 + (юг) , ^Р = —, (Хйд = л/2-адйд-р,

а óño = а óño (0) + Y(z) •

(1)

Основные математические соотношения

Силы, действующие на лопасть в каждом сечении г, раскладываются на составляющие вдоль набегающего потока (сила сопротивления) и ортогонально к нему (подъемная сила):

Fx = CxPb-

Fy = Cypb-

= & (v2 + ю2 z 2

_1 Í

= v2 +Ю2 z 2 (Cy vä - Cx roz).

i+tg 2p Vi+tg 2p

N = Mb = Jb( z)yjv2ä + ю2 z2 (Cy vá - Cx bz) zdz. (4) 2 0

Еще одной немаловажной характеристикой ВЭУ является действующая на нее опрокидывающая сила Рх

Рх = Рх' 51П(аад + а опд ) + РУ С05(аад + а опд ) =

х

= Fx> cos p + Fy> sin p =

= pb t,2 2

(v2 + b2 z 21

C

Cy tgp

Vi+tg 2P л/i+tg 2p

(2)

2 ' y y 2 где b = b(z) - длина хорды профиля. Здесь введены

следующие безразмерные величины: Cx - коэффициент

сопротивления; Cy - коэффициент подъемной силы;

полные размерные силы в (2) приведены на единицу длины вдоль размаха лопасти [1, 2].

Найдем некоторые характеристики ВЭУ

Fy = Fy. sin(a^ +aдяд)-Fx. cos(a^ +a6m) =

= Fy' cos p -Fx' sinp =

Крутящий момент определяется по формуле

Ь Р Ь I 2 2 2

М = \RyZdz = — |Ь(г)д/У2 + ю г (СуУ^ - Схюz)zdz, (3)

0 2 0

где Ь - длина лопасти. Здесь для простоты принято, что корневое сечение лопасти находится вплотную к оси, так что нижний предел интегрирования вдоль хорды можно для определенности положить равным нулю.

Очевидно, что при очень малых ю подынтегральное выражение (3) положительно. Раскручивание приводит к увеличению ю, при этом второе слагаемое, содержащее ю в виде множителя, работает с отрицательным знаком, т.е. против вращения. При достижении критического значения ю интеграл (3) обращается в 0, и дальнейшее раскручивание лопасти становится невозможным. Таким образом, при вращении лопасти ВЭУ существует внутренний механизм, гарантирующий, что лопасть не сможет набирать угловую скорость бесконечно.

Важной характеристикой ВЭУ является мощность

= -у-\/ V2 +ю2 г 2 (Сх vá + С y юг).

Эта характеристика определяет механическую нагрузку на мачту ВЭУ и, соответственно, должна быть ограничена.

Расчет коэффициентов Сх и Су методом конечных элементов

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

методом конечных элементов (МКЭ) в программном пакете ANSYS CFX 11.0. Решалась двумерная задача стационарного обтекания крылового профиля равномерным потоком воздуха при заданных углах атаки. Использовалась структурированная гексаэдри-ческая сетка, построенная в универсальном сеточном пакете ANSYS ICEM CFD. Количество элементов вычислительной сетки - 59400. При расчете использовалась модель вязкой несжимаемой жидкости.

Численный расчет проводился для крылового профиля NACA 4412. Профили NACA, к которым относится четырехзначный профиль NACA 4412, строятся по определенным правилам с использованием аналитических выражений, описывающих выпуклость (кривизну) средней линии профиля и его толщину по всей длине [3]. В номере профиля содержится информация о его кривизне m (в первой цифре номера), положении максимального утолщения p (во второй цифре), максимальной толщине t (в двух последних цифрах номера). У профиля NACA 4412, согласно этим правилам, максимальная толщина составляет 12 % от его длины; средняя линия имеет выпуклость, равную 4 % длины; максимальное утолщение находится на расстоянии, равном 40 % длины от головной части (m = 0,04b, p = 0,40b, t = 0,12b).

Уравнения, описывающие форму средней линии профиля, имеют вид

m(2px-x2), 0 < x < p

P

m

(i - p)2

(i - 2p) + 2px - x

p < x < i.

Здесь х е[0,1] - безразмерная координата, ось которой направлена вдоль хорды профиля. Толщина профиля изменяется по закону

yt =-X

0,2

<(0,2969Vx -0,i260x-0,35i6x2 + 0,2843x3 -0,i0i5x4).

v

2

2

c

t

В параметрической записи уравнения верхней и нижней линий контура примут вид

Хи = x - yt sin0, Уи = yc + yt cos 0,

xL = x + yt sin0,y¿ = Ус - yt cos 0,0 = (^j • Результаты расчета приведены в табл. 1.

Таблица 1

Зависимость коэффициентов Cx и Cy от угла атаки

Углы атаки Cx Cy

-4;4 0,00852; 0,0260 0,0624; 0,762

-3;3 0,00765; 0,0211 0,139; 0,673

-2;2 0,00777; 0,0169 0,223; 0,583

-1;1 0,00878; 0,0134 0,310; 0,491

0 0,0107 0,400

5;6 0,0313; 0,0372 0,849; 0,930

7,8 0,0433; 0,0494 1,007; 1,075

9,10 0,0548; 0,0590 1,1863; 1,153

11,12 0,0586; 0,0547 1,126; 1,051

Построим аналитические аппроксимации зависимостей Сх (а ¿^) и Су (а¿¿) с точностью,

соответствующей точности расчетов по МКЭ. В простейшем случае приближения с тремя членами степенного и тригонометрического типа имеют вид

Cx (aó ) = 0,0143 + 4,46-10"5 aáó + 4,35-10-8a cv (aáó ) = 0,444 +1,58 -10

->-8 2

-3

aáó - 7,67-10-7aáó , -3,

(5)

Cx (aáó ) = 0,0324 + 0,0117sin(3,9 -10-3 aáó ) -

-0,0229отs(3,•10-3а¿^ ), Су(а¿,д ) = 0,606 + 0,442sin(3,9•Ю-3а¿,д )--0,257cos(3,9•Ю-3а¿1д ), где угол а¿¿ - в радианах. При этом приближение (5) обеспечивает точность с относительной погрешностью в 7 %, (6) - в 5. Более сложные аппроксимации обеспечивают точность с относительной погрешностью менее 1 %.

Обратная задача об оптимальной геометрии лопасти турбины ВЭУ заключается в максимизации выражения для мощности (4) [4]. Решение данной оптимизационной задачи проводится при помощи генетического алгоритма.

Выбор генов генетического алгоритма и построение целевой функции

Будем считать, что v¿ и ю постоянны; оптимизация производится по функциям а ж (г) и

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

-3,

(6)

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

м м

а ¿„б (г)= Еаггг, Ь(г) = £ ЬГ. (7)

г=0 г=0

Подставляя (7) в (1), а затем полученные соотношения в (5), попутно переводя градусы в радианы, получим выражения для коэффициентов подъемной силы и сопротивления:

Cx = Cx (z, a, ) = 0,0143 + 0,46080 x

л/2 - 2a,z - arctan i=0

v v« у у

л 1 +

(

+ 4,633200

M

л/2 - Zaizi - arctan i=0

v v« У У

-2

Cy = Cy (z, a, ) = 0,444 +16,380 x

(8)

M

л/2 - 2aizi - arctan i=0

v v « У У

1

(

- 81,64800

M

л -2

л/2 - Za,z - arctan i=0 i

v v « У У

-2

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

образом, можно разбить интеграл в (4) на сумму интегралов, вынести из каждого интеграла неизвестные коэффициенты, вычислить интегралы

(¿)сЪ на этапе построения целевой функции, что

существенно сократит время работы алгоритма.

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

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

N = ^Г+1 - Г)/2(Ь ^2 +ю2г2 х

2 г=0

X (Су (, а1 )у3 - Сх (21, а1))юг1)+

/ 2 2 2 + Ьг+ю гг'+1 (Су(гг+Ьаг+1^3 -

- Сх (^г+1, аг+1))югг+1)гг+1),

где - координаты сечений лопасти; п -количество сечений.

x

2

x

Функция N обладает низкой вычислительной стоимостью и позволяет использовать как приближения Сх и Су из (5) и (6), так и более сложные

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

Построение генетического алгоритма

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

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

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

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

ген потомка =

ген родителя 1, а < —0,25

а • ген родителя 1 + (1 — а) • ген родителя 2, ае [—0,25;1,25].

ген родителя 2,

а > 1,25

В случае если полученный ген не удовлетворяет ограничениям, число а выбирается заново. Приведенный алгоритм является смесью дискретной и промежуточной рекомбинации. Применительно к данной задаче использование этой формулы приводит к существенному сокращению количества поколений, необходимых для сходимости процесса, которое уменьшается на 25-35 % по сравнению с чистой промежуточной рекомбинацией, описанной в [6].

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

новое значение гена = старое значение гена + а • 8, где знак «+» или «-» выбирается с равной вероятностью; а = 0,5(miхgi — mingi); maхgi -максимальное допустимое значение данного гена; mingi - минимальное допустимое значение данного

т _.

гена, 8 = 2 а^ 2 1; ai = 1 с вероятностью , иначе -

i=1

0; т - параметр, взятый равным 20.

Сходимость: в случае если за 10 % от максимального числа поколений значение целевой функции для лучшей особи изменилось не более чем на 1 %, алгоритм считается сошедшимся.

Результаты вычислительного эксперимента

С применением данного генетического алгоритма проведена условная оптимизация трехлопастной ВЭУ со следующими наложенными ограничениями:

— 4° < аад <12°; 0,05 1 < Ь <0,25 1 . Опрокидывающая сила ограничена сверху значением 60 Н. Скорость набегающего потока - 8 м/с, ю =146,6 рад/с. Полученная в результате расчета результирующая мощность ветроустановки составила 2970 Вт. Значения полученных углов установки и длин хорды приведены в табл. 2.

Таблица 2

Форма лопасти, полученной из генетического алгоритма

Координаты сечения, м Длина хорды, м Угол установки, град.

0,08 0,25 22,3

0,17 0,17 6,3

0,27 0,14 1,77

0,36 0,11 0,546

0,46 0,11 0,258

0,56 0,10 0,118

0,65 0,099 0,316

0,75 0,085 0,382

0,8 0,082 0,440

Литература

1. Седов Л. И. Плоские задачи гидродинамики и

аэродинамики. М., 1996. 448 с.

2. Лойцянский Л. Г. Механика жидкости и газа. М., 1973.

848 с.

3. Jacobs E.N., Ward K.E., Pinkerton R.M. The characteristics

of 78 related airfoil sections from tests in the variable-density tunnel // NACA Report. 1933. №. 460. Р. 58.

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

4. Сумбатян М.А., Бондарчук А.А. Обратная задача об

оптимальной геометрии лопасти турбины ВЭУ // Современные проблемы механики сплошной среды: тр. XIV междунар. конф. Ростов н/Д, 2010. Т. 1. 319 c.

5. Панченко Т. В. Генетические алгоритмы. Астрахань,

2007. 87 с.

6. Muehlenbein H., Schlierkamp-Voosen D. Predictive Models for

the Breeder Genetic Algorithm: I. Continuous Parameter Optimization // Evolutionary Computation. 1993. Vol. 1. P. 25.

15 ноября 2012 г.