CC BY
29Статья поступила в редакцию 23.04.12. Ред. рег. № 1288
The article has entered in publishing office 23.04.12. Ed. reg. No. 1288
УДК 621.548.4
ОПТИМИЗАЦИЯ ГЕОМЕТРИИ ЛОПАСТИ ДЛЯ МАЛЫХ ВЕТРОЭНЕРГЕТИЧЕСКИХ УСТАНОВОК
Ю.Н. Жесткое, М.А. Сумбатян, А.А. Бондарчук
ООО «Персональные энергосистемы» 344116 Ростов-на-Дону, ул. 2-я Володарского, д. 76/23а, оф.213 Тел. (863)2189024, факс (863)2189024, e-mail: es-energy@aaanet.ru
Заключение совета рецензентов: 27.04.12 Заключение совета экспертов: 02.05.12 Принято к публикации: 04.05.12
В статье разработана методика оптимизационного выбора геометрии лопасти малых ветроэнергетических установок. Метод основан на максимизации функционала, которым выражается в виде некоторых квадратур величина ветровой мощности ВЭУ. Используется сочетание инженерной теории ветряка Г.Х. Сабинина и точный расчет на основе конечно-элементного пакета АНСИС. В рамках этой методики удается проанализировать геометрии двух различных лопастей: на основе крылового симметричного профиля NACA-0015 и специального профиля FX-63-137. В результате предложенного метода оптимизации геометрии выбран профиль, причем полученное оптимальное изменение хорды и угла установки (крутки) вдоль размаха впервые позволило достичь значения КИЭВ = 0,52 для 10-киловаттных установок.
Ключевые слова: ветроустановка, лопасть, оптимальная геометрия, профиль, ветровая энергия.
OPTIMIZATION OF BLADE GEOMETRY FOR SMALL WIND TURBINES Yu.N. Zhestkov, M.A. Sumbatyan, A.A. Bondarchuk
Limited Company "Personal Energy Systems" 76/23a Volodarskogo str., Rostov-on-Don, Russia Tel. (863)2189024, fax (863)2189024, e-mail: es-energy@aaanet.ru
Referred: 27.04.12 Expertise: 02.05.12 Accepted: 04.05.12
We propose a new method to optimize blade geometry for small wind turbines. This is based on a maximization of the functional which expresses wind energy in terms of some integrals. We combine an engineering Sabinin theory with exact aerodynamic FEM ANSYS CFX calculations. In frames of the proposed methodology we analyze the geometry of two different blades, made on the basis of symmetric wing airfoil NACA-0015 and a special airfoil FX-63-137. As a result, the proposed optimization technique allows us to firstly achieve the wind energy coefficient = 0.52 for 10 kWt wind turbines, by an optimal distribution of the chord and twist angle over the blade length.
Keywords: wind turbine, blade, optimal geometry, airfoil, wind energy.
Введение
Классическая аэродинамика не дает методики выбора оптимальной геометрии профиля лопасти ВЭУ прямым аналитическим методом. В связи с этим созданы каталоги профилей, имеющих те или иные преимущества по таким характеристикам, как максимальное аэродинамическое качество и число Рейнольдса Яе, ему соответствующее, профильное сопротивление, критический угол атаки и характер срыва, положение фокуса профиля и т.д. Обширные каталоги с результатами продувок были созданы в
NACA [1], в ЦАГИ [2], а также в европейских научных организациях. Речь идет о профилях, разработанных как для летательных аппаратов, так и для ВЭУ. Отличие российских условий применения ВЭУ - преобладание умеренных среднегодовых скоростей ветра 3-4 м/с, что делает актуальным поиск оптимальной конфигурации для условий с малыми скоростями ветра.
В данной работе приведена методика поиска оптимальной аэродинамической конфигурации лопасти ветряной турбины. Проведен численный эксперимент на установке мощностью 10 кВт.
Международный научный журнал «Альтернативная энергетика и экология» № 04 (108) 2012 © Научно-технический центр «TATA», 2012
Описание геометрии турбины и предлагаемый метод расчета
Бетц показал [3], что максимально возможный коэффициент использования энергии ветра (КИЭВ) равен 0,593. На сегодня значение 0,5 превышено только для мегаваттных установок. Для исследуемой ВЭУ выбран ротор горизонтально-осевого типа диаметром 10 метров, установленный за башней, ориен-
тация на ветер - свободная. Расчетная скорость ветра по выходной мощности V = 7,7 м/с. Расчетная частота вращения ротора п = 90 об/мин. Число лопастей I = 2. В работе развивается новый метод проектирования геометрии лопасти, который состоит в сочетании теории Г.Х. Сабинина [4] и методов оптимизации по выходной мощности ветряка. Корректность такого подхода оценивается путем точного 3Б расчета на основе МКЭ. Геометрия лопастей показана на рис. 1.
Рис. 1. Предварительный внешний вид геометрии лопастей Fig. 1. Preliminary draft of blades geometry
Используя схему аэродинамических сил и скоростей из теории Г.Х. Сабинина [4], имеем
W = V (V — и1 )2 + (-юг - Mj)2 =
= (V — = U;
sin р
z„ = ctg ß =
ЮГ + Mj
V-и
(1)
C pb 2 Cypb 2 F, = -^-W ; F, = -^-W2. x 2 y 2
(2)
Ry = Fy, sin ß - Fx, cos ß =
= ^W2 2
C„
C ctg ß
>+ctg2 ß л/1 + ctg2 ß
pbw C
271
+ z~
Kl-VZu ),
(3)
где Ж - относительная скорость набегающего потока; и1 и и1 - компоненты изменения скорости потока за счет взаимодействия с турбиной; 2и - относительная быстроходность в данном сечении г; в - угол между вектором Ж и плоскостью вращения.
Силы, действующие на лопасть в каждом сечении г, раскладываются на составляющие вдоль набегающего потока (сила сопротивления) и ортогонально к нему (подъемная сила):
где ц = Сх/Су - величина, обратная аэродинамическому качеству профиля.
Тогда крутящий момент определяется по формуле (здесь £ < г < Ь - изменение радиальной координаты вдоль размаха лопасти)
М =[ Кгйг =Р[ Ь (г )Ж 2(г )С (г) 1 -Ц(Г) Ыг . (4)
< 2 * V1 + г2(г)
Не все величины, входящие в (4), являются независимыми. Так, имеют место соотношения связи [4]:
ibCy = 8пг-----Ц=
y (1 + e)(1 - e)2 (
где р - массовая плотность жидкости; Ь = Ь(г) - длина хорды профиля; безразмерные величины: Сх -коэффициент сопротивления; Су - коэффициент подъемной силы.
Основной интерес представляет момент результирующей силы Яу, действующей на хорду, относительно оси вращения. Его величина равна
1 ^ 1 +(4«А2 ){(1 - е)/(1 + е)} 2и = 2-2(1—е)-' (5)
где I = 2 - число лопастей; е = - коэффициент торможения потока, в каждом сечении свой.
International Scientific Journal for Alternative Energy and Ecology № 04 (108) 2012
© Scientific Technical Centre «TATA», 2012
Ю.Н. Жесткое, М.А. Сумбатян, А.А. Бондарчук. Оптимизация геометрии лопасти для малых ветроэнергетических установок
Заметим, что коэффициент подъемной силы для данного профиля является известной функцией угла атаки, который также переменен по размаху: а = в -ф, Cy = Cy(a) (где ф - угол крутки).
Введем понятие «прямой» и «обратной» задачи. В прямой задаче предполагается, что известны физические параметры, т.е. скорость ветра V, угловая скорость вращения ю, радиус лопасти L, а также известна геометрия лопасти. Нужно рассчитать мощность, создаваемую аэродинамическим моментом: Ывк = Мю. Ясно, что в прямой задаче известны все геометрические параметры - типа данных в табл. 1. По формулам теории Г.Х. Сабинина полный момент (4), а с ним и мощность, легко вычисляются простым суммированием (интегрированием) вкладов всех сечений по размаху лопасти. Заметим, что если выбрано N сечений, то в прямой задаче аэродинамическая мощность является функцией 2N величин Ь, ф,-, i = 1, ..., N; Е = Е(Ь, ф,). На практике мы принимали N = 10.
Таблица 1
Углы крутки и длины хорд лопасти для профилей NACA-0015 / FX-63-137
Table 1
Twist angle and chord length distribution for airfoils NACA-0015 / FX-63-137
Обратная задача поиска оптимальной конфигурации формулируется следующим образом. Полный набор параметров Ь, ф,-, г = 1, ..., N заранее не известен. Требуется найти их все из условия максимума мощности на ветроколесе = Мю. Обычно физические параметры (скорость ветра, угловая скорость и диаметр ротора) считаются известными, однако в различных постановках к 2N геометрических неизвестных можно добавлять также и какие-то из физических параметров. В любом случае математически обратная задача состоит в поиске максимума функции Е = Е(Ь,-, ф,) при условии, что каждый из входящих параметров задан в своем диапазоне изменения.
Ограничения на значения проектных параметров могут быть разными. Например, можно требовать монотонного убывания хорды b с ростом номера сечения , а можно и не требовать. Так, для первой из двух рассмотренных ниже лопастей это условие выполняется, а для второй - нет. То же самое - по углу крутки. Для последнего можно требовать, чтобы в каждом сечении угол крутки был положительным, а можно и не требовать. Иногда для реализации более экономной технологии изготовления желательно выдержать трапециевидную форму лопасти в плане, тогда закон изменения параметров b должен быть линейным по . При любом подходе максимизация функции Е = Е(Ь,, ф,) происходит в выбранной области изменения проектных переменных.
Для поиска оптимального значения целевой функции мы использовали метод глобального случайного поиска [5]. Эффективность алгоритма основана на следующих двух его особенностях: 1) случайный выбор совокупности искомых параметров в окрестности точек с наибольшими значениями функции происходит более часто, чем в окрестности точек с меньшими значениями; 2) области, в которых выбираются значения случайного поиска набора искомых параметров, постепенно стягиваются к малым окрестностям точек с наибольшими значениями целевой функции. В итоге решение задачи максимизации мощности определяет оптимальные значения всех величин, стоящих под интегралом, как функции от радиуса r.
Построенная таким методом оптимальная геометрия далее тестируется расчетами по МКЭ.
Сравнение с расчетами по МКЭ
Геометрия расчетной области - усеченный конус, ось которого совпадает с осью вращения лопасти. Радиус меньшего основания конуса, расположенного вверх по течению - 15 м, радиус большего основания, расположенного вниз по течению - 20 м, расстояние между основаниями - 40 м. На меньшем основании ставится условие типа Inlet (воздух поступает в область с заданной скоростью через эту поверхность), на большем основании и на боковой поверхности ставится условие типа Outlet (воздух выходит из области через эти поверхности). Построение расчетной сетки выполнялось в универсальном сеточном пакете ANSYS ICEM CFD. Количество элементов в МКЭ - от 1,2 млн до 4,5 млн. На всех конструкционных элементах ВЭУ ставится граничное условие прилипания.
В МКЭ-пакете ANSYS CFX проводится нестационарный расчет при заданной скорости потока воздуха и частоты вращения. Вводится вращающаяся система координат, связанная с лопастью ветроустановки (дополнительные инерционные силы, появляющиеся во вращающейся системе координат, учитываются автоматически). В такой системе лопасть покоится, обтекатель приводного вала вращается, но остается в исходном объеме.
r bb м фь град
0,075 0,560/1,105 28,45/37,37
0,10 0,540/1,175 23,47/33,19
0,20 0,450/1,025 11,56/20,51
0,30 0,386/0,811 6,29/13,46
0,40 0,330/0,658 3,27/9,32
0,50 0,276/0,551 1,14/6,72
0,60 0,236/0,473 -0,24/4,98
0,70 0,198/0,416 -1,49/3,76
0,80 0,172/0,372 -2,26/2,89
0,90 0,150/0,337 -2,84/2,25
1,00 0,130/0,309 -3,22/1,74
Международный научный журнал «Альтернативная энергетика и экология» № 04 (108) 2012 © Научно-технический центр «TATA», 2012
Использовалась модель сжимаемого теплопроводного вязкого газа. В ANSYS CFX это реализуется выбором модели газа AirldealGas в сочетании с выбором модели теплопроводности TotalEnergy. Для учета турбулентности вместо классической к-г модели, которая недостаточно хорошо описывает пограничные слои вблизи стенок, использовалась более точная SST-модель турбулентности Ментера (Shear Stress Transport).
В начальный момент времени ветровой поток движется равномерно и прямолинейно, кроме узкого диска вблизи лопасти, который движется вместе с ней по закону вращения твердого тела. Также в качестве начальных условий могут использоваться результаты расчетов на одном из предыдущих шагов.
Первым шагом при проведении численного расчета является построение геометрии лопасти. Для этого требуется в нескольких сечениях задать форму аэродинамического профиля, его размеры и угол крутки. После этого на полученный каркас лопасти натягивается B-spline поверхность, строится расчетная сетка и загружается в пресолвер ANSYS CFX.
В численных экспериментах были использованы два профиля, показанных на рис. 2. Продольная ось лопасти с профилем NACA-0015 находится на расстоянии 30%, а с профилем FX-63-137 - на расстоянии 31% длины хорды лопасти от передней кромки.
Рис. 2. Общий вид крылового профиля NACA-0015 и профиля FX-63-137 Fig. 2. A sketch of the wing airfoil NACA-0015 and airfoil FX-63-137
Лопасть развернута вогнутой частью к набегающему потоку воздуха. Профиль FX-63-137 изогнут существенно сильнее, чем NACA-0015, и имеет вытянутую острую кромку, что требует более аккуратного построения МКЭ-сетки.
Анализ полученных результатов
Основные результаты расчета приводятся в табл. 1 и 2. Крыловой профиль NACA-0015 был выбран из тех соображений, что на нем отрабатывалась принятая методика расчетов, которые носили тестовый характер. С аэродинамической точки зрения спроектированная лопасть на основе этого профиля не является удачной. В частности, начиная с седьмого сегмента вращающий момент при использовании профиля становится отрицательным, также имеются отрицательные углы крутки. Однако на этой лопасти нам удалось показать соответствие между разными подходами при проектировании.
Таблица 2
Распределение момента вдоль лопасти для профилей NACA-0015/FX-63-137
Table 2
Moment distribution along the blade for airfoils NACA-0015/FX-63-137
Сегменты по длине лопасти, м Вращающий момент, Нм
0,375-0,5 2,209/5,450
0,5-1,0 25,345/47,795
1,0-1,5 28,990/88,903
1,5-2,0 23,388/124,542
2,0-2,5 20,957/150,863
2,5-3,0 4,547/169,505
3,0-3,5 -8,477/193,413
3,5-4,0 -2,998/216,293
4,0-4,5 -10,595/228,650
4,5-5,0 -24,668/128,973
Суммарная величина 58,698/1354,386
Далее мы перешли к проектированию принципиально новой лопасти с высокими аэродинамическими характеристиками на основе профиля FX63-137. Для нее также подтвердилось соответствие между теорией Г.Х. Сабинина и точным расчетом по МКЭ. Полученное оптимальное распределение параметров b,, ф,, i = 1, ..., N вдоль лопасти приведено в табл. 1 и 2.
В результате предложенного метода оптимизации геометрии выбран профиль FX-63-137, причем полученное оптимальное изменение хорды и угла установки (крутки) вдоль размаха впервые позволило достичь значения КИЭВ = 0,52 для 10-киловаттных установок.
Работа выполнена при поддержке ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2013 годы», ГК № 16.516.11.6106.
Список литературы
1. Jacobs E.N., Ward K.E., Pmkerton R.M. The characteristics of 78 related airfoil sections from tests in the variable-density wind tunnel // NACA Report No. 460. Washington, D. C., 1933.
2. Кравец А.С. Характеристики авиационных профилей. М.-Л.: ГИОП, 1939.
3. Glauert H. Airplane propellers. In: Durand WF, editor. Aerodynamic theory. New York: Dover Publications, 1963.
4. Сабинин Г.Х. Теория и аэродинамический расчет ветряных двигателей // Труды НИИ Промышленности, № 482. ЦАГИ. 1931. вып. 104.
5. Жиглявский А. А. Математическая теория глобального случайного поиска. Л.: Изд-во ЛГУ, 1985.
Г'-": — TATA — LXJ
International Scientific Journal for Alternative Energy and Ecology № 04 (108) 2012
© Scientific Technical Centre «TATA», 2012
