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

УДК 533.69, 621.311.24

ОПТИМИЗАЦИЯ ГЕОМЕТРИИ ЛОПАСТИ ВЕТРОТУРБИНЫ, РАБОТАЮЩЕЙ В ТУННЕЛЕ ПОСТОЯННОГО ДИАМЕТРА

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

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

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

Sumbatyan Meglum Albertovich - Doctor of Physical and Mathematical Science, Professor, Faculty of Mathematicians, Mechanics and Computer Sciences, Southern Federal University, Mil'chakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: sumbat@math.sfedu.ru.

Mescheryakov Konstantin Igorevich - Assistant, Faculty of Mathematicians, Mechanics and Computer Sciences, Southern Federal University, Mil'chakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: m.keyran@gmail.com.

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

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

We propose a new method to find optimal geometry of the wind turbine blade, working in a tunnel of constant diameter. The blade is a thin twisted plate with a varying chord. The cross-sections may differ from each other by the angle of setting and the width. The attack angle is defined by the angle of setting, and by the angle which depends upon the wind speed. The fluid is accepted incompressible and inviscid. To solve the problem, we accept the hypothesis of plane sections. Under such hypotheses there is developed an integral equation, and a certain numerical method to solve this equation is proposed. Besides, there is constructed a genetic algorithm, which provides optimization of the blade. Some examples demonstrate the obtained numerical solutions.

Keywords: wind energy, wind turbine, boundary integral equations, genetic algorithm, optimization, wind turbine blade.

Рассмотрим вращающуюся лопасть, представ- диаметра. Жидкость считается идеальной. Принимая

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

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

атаки зависит от угла установки и угла, На рис. 1 изображен однородный набегающий

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

угловой скоростью вращения лопасти. Лопасть лопасти. Ось z направлена вдоль размаха лопасти,

работает в цилиндрическом туннеле постоянного перпендикулярно плоскости рисунка, причем ось x

32

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

* Y

~[1](s, y) = }ф[1](х, y)elsxdx =

—ТО

Фо1](У) , фХо(У) 0 м]

2 — í — фХХ( x y)dx,

S —то s

~[2](s, y) = 1ф[2]( x, y)eisxdx =

о

ф02]( у) фХ2]( у) 7 eisx [2], ,,

----Х02--\—Г фХХ( X, y)dx.

iS S о s

Здесь и далее нулевым индексом будем обозначать функцию, взятую при х = 0.

Таким образом, после применения преобразования Фурье по переменной х (1) принимает следующий вид:

— s 2~[1,2](s, y) +

52~[1,2](S,У) ____[ 1,2]

5у2

= у [ (s, y),

V[1'2](5, j) = j) - «p01,2](y)J.

Общее решение данного уравнения

ф[1'2](5,y) = С[[1'2] • ch(jy) + C2[12] • sh(jy) +

(2)

Рис. 1. Геометрия сечения и набегающего потока

Таким образом, относительная скорость потока образует угол р с плоскостью хг. Тогда

Г. юг п

= —, аао = - -а ойо "р .

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

Вывод граничного интегрального уравнения

Пусть ф( х,у) является потенциалом потока,

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

ф /и// (х,у) = их + ф(х,у)>

^^ = 0. (1)

Эх2 су2

Разобьем исходную область на две подобласти:

1) х <0, у < Я;

2) х > 0, у < Я.

Применим к (1) преобразование Фурье по переменной х отдельно в каждой из подобластей и возьмем интеграл дважды по частям:

+ ^ (У у)К

25 _к

Его справедливость проверяется прямой подстановкой в (2); С1 и С2 - произвольные постоянные.

Проведем обратные преобразования:

1 да /

ф[1,2](х,у) = ^ | ^ (с11,2] • е^у) + С2[1,2] • +

+ ^ Я ^[1'2](*,У)(зЬ(5(у -25 - я )

Из условий непроницания на твердых

поверхностях туннеля и лопасти и условий

непрерывности скорости вне

следующие граничные условия:

Эф[1'2]

лопасти имеем

5у

= о,

y=± R Vn =rnz • sin ауст,

Vn = U cos ауст + px cos ауст + py sin ауст и ~ U + px +фуЯуст, У e[-L, L],

(Ox (У) + P0y (У)ауст = ^«уст -U. У e[-L L] (3)

p0X](y)=pOX1^ p0>)=PO'M

y e (-R,-L) u (L, R).

Здесь L , R - длина лопасти и ширина цилиндра в данном сечении.

Решение, удовлетворяющее граничному условию на твердых стенках y = ±R, имеет вид

Ф[1,2]( X, У) = —f 4л

R J e~'sx y[1,2](s,У) х

- — R—то

х, ch(s(y + y)) + ch(s(2R—|y — у|)) .

sh(2sR) 1

о

Дифференцируя (4) по переменным x и у , получим

; R ж

фХ1'2](х,y) = -L I I g—y[1,2](s,У);

4л _ R—ж

'ch(s(y + y)) + ch(s(2 R—|y — y|))

sh(2 sR)

(5)

dsdy.

J1,2V

(6)

4л

—R—ж

XI ch(s(y+y))+ch(s(2R—|y _y|)) Idsdy,

sh(2 sR)

J!,2]

_ 1

R ж

(0,y) = +-1 I I sф[)1'2](y)X

4л — R—ж

( ch(s(y + y)) + ch(s(2 R— | y — y |)) [ sh(2 sR)

i R ж

фУ1'2^,y)= ±-L I |Ф0Х2](У) X

4л —R—ж

f sh(s(y + y)) — sh(s(2R— | y - y| ))sgn(y - y)

dsdy,

(7)

(8)

dsdy.

^ sh(2sR)

Тогда граничное условие (3) можно представить в виде

Ф0.

0X2I(у) + U — ® zaóñó = +

R

j-x i ф01хх2](у)x

4л — R

+ у)) - Я-1 у -у |))8ВР(у - у) К , (9)

-И ¡¡Ь(2яЯ) J ^

Введем оператор К следующим образом:

(К/)(у)=^ х 4л

х 4/(у) у + у)) - 8^(2Я- | у - у | )>^(у - у)

-Я -у вЬОЯ) '

Тогда соотношение (9) принимает вид (I ± «устК)р|[1;2](у) = -и , у е [-Ь,Ь].

Оставляя лишь главный член асимптотики в разложении (I ±аустК)-1 в ряд Неймана, получим

у)-

(i ± «уСТк) 1 (®zayCT — и)

í azaуст —U ,

т.е.

граничное условие имеет вид

Ф^СУ) = ®zaуст - U, y e [-L,L],

Ф0Х (У) = фОХ (У), У е (—R-L) ^ (L, R).

Следовательно,

ф0Х(у) = ф02Х](y), У еьRR].

(10) (11)

что

Из (2) и (9) следует,

Y[1](s, У) + v[2](s, у) = Цф01]( У) -Фо2](У)) •

Используя равенство (11) в (8), получаем

ф[1](0, у) = —ф[2](0, у) .

1 Я у

«ругч*,у) = I I е-"у)х

у 4л -я-у

х Г 8Ь(д(у + у)) - 8Ь(д(2Я~1 у - у 1 )^п(у - у) V

Из граничных условий следует

рОХ (у) + Р0у (у)ауст = рО*1 (у) + Р02 (у)ауст, х = 0

Подставив в (5) и (6) х = 0, используя соотношение и проведя преобразования, учитывающие четность и нечетность подынтегральных выражений относительно 5 , из формул (4)-(6) получим

р[,2](0,у) = ±± Чур^х

Пусть

g( У) = ф[2](0, у)—ф[1](0, у)

тогда

ф[2](У) = . Из уравнения (7)

фХ2](0, у) =

= ± J fsg(y)ch(s(^ + У)) + ch(s(2^-|У-У|))dyds.

8л

—ж—R

sh(2 sR)

Подставим сюда граничное условие (10). Заметив, что g(у) = 0 при |у| > Ь, получим

iL ж

- Ig(y) I

ch(s(y + y)) + ch(s(2 R—|y — y|))

8л

sh(2 sR)

x dyds = a zayст — U, У e[—L, L].

Для вычисления внутреннего воспользуемся известной формулой [1]

J s ch(as) ( л ( naУ 2

i--—1 ds = I — II cos— I , b >| a |,

0 sh(bs) ^ 2b Д 2b J

что позволяет преобразовать данный интеграл

» ch(s(y + y)) + ch(s(2R-|y-y|))^ =

0 _ sh(2 sR)

(12) интеграла

2

8R 2

1

1

2 л(У + y)

cos

2 л(2R— | y — y|)

4R

cos

4R

1

sin2 Л(У — y)

+ -

1

2 Л(y + y)

cos

8R

sin _____

4R 4R

Тогда интегральное уравнение (12) принимает вид

64R

I g (y)

1

1

. 2 л(у — y)

sin2 —^—— 4R

2 Л(у + y)

cos

4R

dy =

= ю zауст - U,

У e[-L, L] • (13)

Проинтегрировав по частям интеграл в (13), получим

1 L

T6R1g(y)

л(У — у) Л(У + у) ctg (У — — tg. (y

4R

4R

dy = a zayст — U

У e[-L, L] •

Согласно уравнению Бернулли

p[2] - p[1]_ (f[2]-V[1]) =

P 2

(VX2]+fX1])(FX2] -FX1])+(vy2]+№ -)

(14)

. (15)

x

L

ж

s

+

2

л

л

+

L

Обозначим:

+ у]1] )= 2П + ф02]

- )=ф02х]-ф0],

Бз.= У2]+уу1 )=ф0у +<$,

*4.4у2] - уу1] )=ф02у-ф0].

Для нахождения разности давлений достаточно определить величины ; S4 находится

непосредственно численным решением уравнения (14); Б2 можно получить путем преобразования граничных условий (3)

ф02*] (у) - ф101 (у) = -ф0у (у)ауст - ф0у (у)ауст) > отсюда Б2 = -Б4ауст.

Подставляя в уравнение (8) и учитывая, что

f sh(as) , n na f—^^ ds = —tg—. 0 sh(bs) 2b 2b

b >| a | [1], получим

S3 = 9^(0, У) + ф[2](0, y)

.[2]/

- tg

n(2R- | y - -1 - -)

4R

i«- +

-¡У+укг)-«(ф^ - ]*.

Здесь учтено, что вне лопасти ф^Су) - ф^Су) = 0 . Значения Бз находятся численно после численного решения уравнения (14).

Для нахождения Б1 преобразуем граничное условие (3)

фсх](у) + ф0х](у) = 2(югауСт - и)--фу]( у)ауст +ф02у]( у)ауст

=2(югауст-и)-БзауСТ.

Следовательно, из уравнения Бернулли (15) имеем p[2] _ p[1] = р S1S2 + S3S4 =

= Р-

(2(югауст - С/)- Sз«уст)(- S4ауст)+ SзS.

2

Метод численного решения полученного ГИУ и результаты

Интегральное уравнение (14) является сингулярным; для его решения можно использовать численную схему, данную в [2]. Проведем дискретизацию (14) на каждом сечении лопасти. Тогда для каждого сечения получим следующую систему линейных уравнений:

тк^у(у (4Я(у - у;) ]-(у;+у) ]]х

х (у;+1 - у ; ) = югауст - и, ] = 0. п,

где п - количество точек разбиения; у =2ЬИп; у^ = 2Ь(] + 0,5)/п ; Ь - ширина лопасти в данном сечении (рис. 2).

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

Fj = Fj cos p, Ffl = Fj sinp.

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

H

M = J Fj zdz , где H - длина лопасти.

о

Рис. 2. Расположение уг- (х) и уj (о) на лопасти

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

Данная система решается с помощью библиотек линейной алгебры.

Проведена оптимизация геометрии лопасти с помощью генетического алгоритма. Геометрия лопасти представлена в виде

= а0 (1 + а1аз г + а2 (1 -а1 )г) (1 + а 2 г )(1 + аз 2)

= Ь0(1 + Ьфз г + Ь2(1 - Ь1 )г) ( ) (1 + Ь2 г )(1 + Ьз г) ,

где г - безразмерная величина, представляющая собой расстояние от оси вращения до текущего сечения лопасти; а( г) - функция, задающая угол установки в

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

Ген потомка = <

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

Стратегия мутации взята из [3].

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

Пример. Рассчитывалась лопасть размахом 3,8 м, работающая в туннеле диаметром 4 м. Ширина хорды ограничена до 0,4 м. Скорость набегающего потока -7,8 м/с, угловая скорость вращения - 20 рад/с. В результате работы генетического алгоритма получены оптимальные параметры а, и Ь, (табл. 1).

Таблица 1

Коэффициенты для оптимальной лопасти

i ai b,

0 0,04 0,249

1 0,46 0,330

2 0,39 0,620

3 0,21 0,315

Соответствующие значения угла установки и ширины лопасти приведены в табл. 2.

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

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

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

Вращающий момент составил 592,24 Нм.

Таблица 2

Геометрия лопасти

Z, м b, м

0 2,31 0,249

0,38 2,07 0,21

0,76 1,89 0,19

1,14 1,73 0,17

1,52 1,60 0,15

1,90 1,49 0,14

2,28 1,39 0,13

2,66 1,30 0,12

3,04 1,23 0,11

3,42 1,16 0,10

3,80 1,10 0,10

Литература

1. Прудников А.П., Брычков Ю.А., Маричев О.И. Интегралы

и ряды: Элементарные функции. М., 1981. 298 с.

2. Белоцерковский С.М., Лифанов И.К. Численные методы

в сингулярных интегральных уравнениях и их применение в аэродинамике, теории упругости, электродинамике. М., 1985. 78 с.

3. Muhlenbein H., Schlierkamp-Voosen D. Predictive Models

for the Breeder Genetic Algorithm - I. Continuous Parameter Optimization // EVOLUTIONARY COMPUTATION. 1993. Vol. 1. P. 25-49.

9 апреля 2014 г.

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

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

ген родителя 2, а > 1,25