CC BY
14КОМПЛЕКСНОЕ МОДЕЛИРОВАНИЕ ВЭУ С ВЕРТИКАЛЬНОЙ ОСЬЮ ВРАЩЕНИЯ
INTEGRATED MODELING OF VERTICAL-AXIS WIND TURBINES
Статья поступила в редакцию 04.01.10. Ред. рег. № 686 The article has entered in publishing office 04.01.10. Ed. reg. No. 686
УДК 621.458.5
КОМПЛЕКСНАЯ ПРОГРАММНО-МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ВЕТРОЭНЕРГЕТИЧЕСКОЙ УСТАНОВКИ
О. В. Матвеенко
OOO «ГРЦ-Вертикаль» 456300 Челябинская обл., г. Миасс, Тургоякское шоссе, д. 1 Тел.: +79123171805, факс: (351) 2647694; e-mail: src-vertical@gmail.com
Заключение совета рецензентов: 15.01.10 Заключение совета экспертов: 20.01.10 Принято к публикации: 25.01.10
В работе рассматривается программно-математическая модель вертикально-осевой ветроэнергетической установки с ротором типа Н-Дарье, разрабатываемой ООО «ГРЦ-Вертикаль», и обсуждаются некоторые результаты численного моделирования.
Ключевые слова: комплексная программно-математическая модель, вертикально-осевая ветроэнергетическая установка, численное моделирование.
COMPLEX PROGRAMMING MATH MODEL OF WIND POWER UNIT
O.V. Matveyenko
"SRC-Vertical", Ltd. 1 Turgoyaksky road, Miass, Chelyabinsk reg., 456300, Russia Tel.: (912) 317-1805, fax: (351) 264-7694; e-mail: src-vertical@mail.ru
Referred: 15.01.10 Expertise: 20.01.10 Accepted: 25.01.10
The article shows the programming mathematic model of vertical axis wind turbine with Darrieus rotor developed by SRC-Vertical company, with discussion of some results of numeric modeling.
Keywords: complex programming math model, vertical axis wind turbine, numerical modeling.
Введение
Главной конструктивной особенностью вертикально-осевых ветроэнергетических установок (ВЭУ) является ротор, преобразующий энергию ветра в энергию вращения вала. В работе в качестве примера рассматривается ротор ВЭУ-3 разработки ООО «ГРЦ-Вертикаль», внешний вид которого показан на рис. 1. Он устанавливается на мачте, удерживаемой в вертикальном положении тросовыми оттяжками. В ступице ротора расположен электрогенератор. Ротор ВЭУ-3 обладает способностью к самостоятельной раскрутке при скорости ветра ~ 4 м/с. Установленные на спицах ротора аэродинамические щитки, поворачиваемые центробежным регулятором поперек потока, предотвращают раскрутку ротора до опасных скоростей при любом ветре. Возможность самопроизвольной остановки ротора устраняет цифровая система управления ВЭУ, которая ограничивает нагрузку электрогенератора в зависимости от скорости ветра и скорости вращения.
Рис. 1. Ротор ВЭУ Fig. 1. Wind turbine
International Scientific Journal for Alternative Energy and Ecology № 5 (85) 2010
© Scientific Technical Centre «TATA», 2010
Работа ВЭУ сопровождается большим числом процессов различной физической природы. Это явления случайного изменения скорости ветра, аэродинамического воздействия ветра на ротор ВЭУ, колебаний и перемещений элементов конструкции, протекания электрических токов в цепях генератора и нагрузки, прохождение информационных потоков в цифровой системе управления. Эти процессы оказывают значительное влияние друг на друга и совместно определяют характер работы ВЭУ и эффективность преобразования энергии ветра в электрическую энергию. Общепризнано, что эффективным и наиболее дешевым способом исследований подобных комплексных объектов является численное моделирование с помощью ЭВМ. Ниже представлен инструмент такого моделирования - программно-математическая модель ВЭУ и некоторые результаты моделирования. Разработка модели ВЭУ проводилась для решения таких задач, как оценка сил и моментов, действующих на элементы конструкции ВЭУ, определение производительности ВЭУ при действии переменного ветра, отработка алгоритмов цифровой системы управления, оценка работы различных устройств, ограничивающих скорость вра-
щения ротора, подбор наилучших значений жесткости мачты и упругости оттяжек, определение возможности балансировки ротора в полевых условиях, оценка уровня вибрационного воздействия ВЭУ на фундамент.
Модель ВЭУ
Программно-математическая модель ВЭУ реализована в программной среде визуального программирования У1$81т и имеет блочную структуру, показанную на рис. 2.
На рис. 2 используются следующие обозначения: V - скорость ветра; п - последовательность импульсов напряжения; Ма - аэродинамический момент; М, Е - момент и сила; Мэ - электромагнитный момент; I, и - ток и напряжение; V - угловая скорость; к - коэффициент ШИМ модуляции; г - модулированный сигнал.
Каждый блок схемы, изображенной на рис. 2, моделирует отдельную подсистему ВЭУ. Это позволяет формировать модели ВЭУ разной структуры, подключая одни и отключая другие блоки, в зависимости от задач моделирования.
Рис. 2. Блок-схема ВЭУ Fig. 2. Block diagram of wind power unit
Модель ВЭУ содержит модель переменного ветра, модель аэродинамики ротора ВЭУ, модель центробежного ограничителя скорости вращения ротора, модель электрогенератора, модель упругой конструкции мачты ВЭУ, модель цифровой системы управления и ее алгоритмов. В качестве математической модели используются дифференциальные уравнения, описывающие механические и электрические явления, сопровождающие работу ВЭУ. Эти уравнения интегрируются по времени. Полученные реализации движения ВЭУ используются для последующего анализа.
Модель ветра
Модель ветра обеспечивает характеристики ветра, задаваемые нормативными документами [1]. Скорость ветра представляется в виде V(t) = V0(t) + + dV(t). Систематическая составляющая V0(t) описывает среднюю скорость ветра и задается в программе как функция времени. Так задаются нормированные порывы и буревой ветер. Статистические характеристики динамической составляющей ветра dV(t) задаются функцией спектральной плотности, например в виде [1]
65
Международный научный журнал «Альтернативная энергетика и экология» № 5 (85) 2010 © Научно-технический центр «TATA», 2010
N ^
ние приближенного равенства Бу(т) = ^ |(/'ю)| . В
1
этом случае для моделирования величины дУ({) необходимо использовать N статистически независимых источников белого шума Ы, действующих на N формирующих фильтров, выходные сигналы которых суммируются. Удовлетворительное приближение рекомендуемой спектральной плотности было получено при N = 3. На рис. 3 для иллюстрации приведена зависимость скорости ветра, полученная описанным выше способом (верхняя кривая). На этом же рисунке приведена зависимость скорости ветра, полученная измерением с помощью анемометра (нижняя кривая). Обе зависимости соответствуют одной средней скорости У0 = 4 м/с. Для получения разброса расчетной скорости ветра, равного разбросу измеренной скорости, пришлось проводить расчет при увеличенном коэффициенте шероховатости. Это объясняется тем, что измерения проводились не на ровной открытой местности, а на опушке леса вблизи зданий.
10 8
0 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000 3250 3500
Время t, с
Рис. 3. Скорость ветра на интервале 1 час Fig. 3. Wind speed on an interval 1 hour
* (")=3 * ш 2П
2 f / т \2 Л
, , ю L 1 +
2n V0
-4/3
(1)
где У0 - средняя скорость ветра; =6к0К02 - дисперсия; Ь - масштаб продольной компоненты атмосферной турбулентности; к0 - коэффициент шероховатости; ю - круговая частота.
Для моделирования динамической составляющей скорости ветра ёУ в зависимости от времени используется метод формирующего фильтра [2]. Используемые в [1] выражения для функций спектральной плотности не являются дробно-рациональными функциями, и метод точного построения формирующего фильтра, приведенный в [1], не может быть использован. Поэтому в обсуждаемой модели используется приближенное решение, когда спектральная плотность аппроксимируется суммой дробно-рациональных функций, для каждой из которых может быть найден свой формирующий фильтр. Практически это сделано с помощью специальной процедуры подбора нескольких частотных характеристик Ж&ю) (к = обеспечивающих выполне-
Модель аэродинамических характеристик ротора ВЭУ
В модели силового воздействия ветра на ротор ВЭУ сила лобового сопротивления Ех и подъемная сила Е7, действующие на каждую лопасть, вычисляются по формулам
¥х = Сх (а)д£ ; (2)
Еу = Сг (а)д£ , (3)
где а - угол атаки; Сх (а) и С7 (а) - аэродинамические коэффициенты; д = р¥2/2 - скоростной напор; р - плотность воздуха; £ - площадь лопасти. Угол
66
атаки для каждой лопасти вычисляется индивидуально, исходя из ее положения на роторе, угла поворота ротора вокруг вертикальной оси ф, угловой скорости ротора ю и мгновенного значения скорости ветра V. При этом учитываются возмущения набегающего на лопасть потока, вносимые другими лопастями. Возмущения учитываются в виде поправок к скорости набегающего потока - индуктивных скоростей. Аэродинамические коэффициенты Сх (а) и С7 (а) вычисляются в программе по табличным зависимостям. Расстояние от передней точки хорды лопасти до точки приложения аэродинамической силы вычисляется по формуле Хп = Св (а)Ъ, где Ъ -хорда; Св (а) - коэффициент центра давления, зада-
International Scientific Journal for Alternative Energy and Ecology № 5 (85) 2010
© Scientific Technical Centre «TATA», 2010
ваемый табличной зависимостью от угла атаки. Аналогичным образом определяются аэродинамические силы, действующие на щитки центробежного ограничителя скорости вращения ротора и точки их приложения. Силы и моменты, действующие на отдельные лопасти, используются для вычисления всех суммарных сил и моментов, действующих на ротор, и вычисления аэродинамической мощности ротора в каждый текущий момент времени.
В программе используются два варианта модели вращения ротора.
В первом варианте модели определяется пространственное вращение ротора, рассматриваемого как твердое тело. Угловое движение ротора определяется интегрированием динамических уравнений Эйлера, записанных в связанных с ротором осях
J х dю/dt = -юх (J хю) + M, (4)
и кинематических соотношений
dу/dt = (юх cosф-ю7 sin9)/cos 6 ; (5) d6/dt = юх sin ф + ю7 cos ф ; (6)
dф/dt = юг + (ю7 sinф-юх cosф^ш6/cos6, (7)
где ю - вектор угловой скорости; J - тензор инерции ротора; у, 6, ф - углы ориентации осей ротора относительно связанной с фундаментом, неподвижной системы координат; M - вектор момента внешних сил, определяемый суммой аэродинамического момента и момента нагрузки со стороны электрогенератора.
Во втором, упрощенном варианте моделируется вращение только вокруг вертикальной оси.
При возрастании угловой скорости ротора до предельной разрешенной величины аэродинамические щитки, размещенные на его спицах, поворачиваются вокруг осей спиц под действием центробежных сил. Это создает тормозящий момент, который вычисляется по известным аэродинамическим характеристикам щитков. В качестве грузов центробежного регулятора выступают сами щитки. В модели предусмотрен и альтернативный вариант, когда для ограничения скорости вращения ротора вместо специальных щитков используются поворотные лопасти ротора. Ниже описана модель альтернативного варианта как менее громоздкая.
Движение системы «лопасть, груз центробежного регулятора, пружина-амортизатор» описывается дифференциальным уравнением вида
Jd2у/dt2 = МТр + MПр + MГр + MA + MУ1 + MУ2, (8)
где J - суммарный момент инерции всей системы, приведенный к центру вращения лопасти с прикрепленным грузом; МТр - момент трения от амортизатора; МПр - момент пружины; МГр - момент центробежных сил, действующих на груз; МА - аэродинамический момент, вращающий лопасть; МУ1 и МУ2 -
моменты, действующие со стороны ограничителей поворота лопасти. Все моменты вычисляются относительно центра вращения лопасти.
Модель генератора
Генератор является синхронной явнополюсной электрической машиной. Магнитный поток в зазоре создается постоянными магнитами. Трехфазная обмотка статора соединена в звезду с заземленным центром. К каждой обмотке фазы подсоединена омическая нагрузка. Постоянные магниты, расположенные на роторе, создают при его вращении переменный магнитный поток (потокосцепление ¥), который индуцирует в фазовых обмотках переменное напряжение:
Евр =-ЭТ/&. (9)
Форма зависимости потокосцепления от угла поворота ротора определяется геометрией обмотки и постоянных магнитов и задается в программе табличной зависимостью. Токи в фазах генератора определяются интегрированием дифференциальных уравнений, выражающих закон Ома для каждой фазы:
0 = 11Я1 + d/dt(Ь' + М12' + М13/3) + &; (10)
0 = 12Я2 + d|dt(Ь212 + М2111 + М23/3) + 2|dt; (11)
0 = /3К3 + d|dt(Ь' + М32/2 + М31/1) + &¥3/dt, (12)
где Ь1 = Ь2 = Ь3 - индуктивность фазы; К1 = Я2 = К3 -сумма активного сопротивления фазы и нагрузки фазы; /ь '2, г'з - мгновенные значения токов в фазах. Поскольку обмотки одинаковы и развернуты относительно друг друга на угол 120 градусов, можно с достаточной точностью положить
¿1 = ¿2 = ¿3 = Ь; (13)
М12 = М13 = М21 = М23 = М31 = М32 = -Ь/2. (14)
Вращение ротора практически не изменяет собственных и взаимных индуктивностей фазовых обмоток, поэтому
дЬк/ дt = 0 (к = 1, 2, 3); (15)
дМы/ дt = 0 (к = 1, 2, 3; п = 1, 2, 3). (16)
Мгновенное значение мощности одной фазы электрогенератора определяется по формуле
Р = ^К . (17)
Полная мощность генератора определяется суммой Р = Р, + Р 2 + Р 3. (18)
Электромагнитный момент на валу генератора выражается формулой
Мм = Р>г . (19)
Международный научный журнал «Альтернативная энергетика и экология» № 5 (85) 2010 © Научно-технический центр «TATA», 2010
Рассмотрим второй, упрощенный вариант модели электрогенератора.
Амплитудное значение ЭДС вращения Евр определяется формулой
Евр = Аю г, A = const.
(20)
Ток фазы i, в пренебрежении индуктивностью обмоток, определяется как
i = Евр/(R + R),
(21)
где Я1, Я„ - активные сопротивления фазной обмотки и нагрузки фазы.
Мощность генератора находится по формуле
P = KSE2J( R + Rn).
(22)
где К - коэффициент, учитывающий осреднение.
Электромагнитный момент (момент нагрузки) на валу генератора выражается как
М = Рг/юг . (23)
Модель системы управления
Одной из основных задач системы управления является регулирование снимаемой с электрогенератора мощности с целью получения наибольшей средней производительности ВЭУ. Снимаемая мощность PC определяется исходя из текущих значений скорости ветра и скорости вращения ротора. Команды, вырабатываемые СУ, передаются координатору мощности, который регулирует снимаемую с генератора мощность PC с помощью ШИМ модуляции. Работа цифровой СУ характеризуется запаздыванием, возникающим из-за того, что ее команды вырабатываются в дискретные моменты времени. Дополнительное запаздывание возникает из-за необходимости обработки показаний датчиков, полученных на некотором временном интервале. Запаздывание команд СУ и погрешность их выработки моделируются с помощью встроенных средств VisSim; «линий задержки», генераторов шумов, квантователей. Модели алгоритмов СУ учитывают реальную вычислительную производительность применяемых процессоров. Текущее состояние ВЭУ система управления всегда оценивает с некоторой погрешностью и запаздыванием. Поэтому при моделировании работы СУ к точным параметрам состояния, получаемым численным моделированием, применяется процедура задержки и добавляется случайная составляющая, задаваемая с помощью датчиков случайных чисел. Если запаздывание сигнала, величина и закон распределения его погрешности неизвестны, то модель СУ содержит алгоритмы измерений, обработки и оценки текущего состояния ВЭУ, которые моделируют реальную погрешность и запаздывание.
Модель мачты
Мачта может быть представлена двумя вариантами моделей.
Первый вариант - жесткий стержень, шарнирно закрепленный на фундаменте и удерживаемый в вертикальном положении оттяжками. Движение мачты в каждой из двух взаимно перпендикулярных плоскостей описывается уравнением вида Jd2 9 / Ж2 = -М+ М ,
где J - момент инерции мачты с ротором; Мр - момент от оттяжек; - момент сил тяжести.
Второй вариант - последовательность упруго состыкованных друг с другом жестких элементов. Каждый элемент имеет массу, тензор инерции и шесть степеней свободы. Каждый элемент имеет жесткий каркас, состоящий из граней прямой треугольной призмы с равносторонним треугольником в основаниях. Элементы стыкуются основаниями. Вершины оснований соседних элементов попарно соединены упругими связями, стремящимися сблизить эти вершины. В недеформированном состоянии верхнее основание нижней призмы совпадает с нижним основанием верхней призмы. На рис. 4 изображены каркасы двух элементов в случае осевого растяжения мачты. Упругие связи показаны пунктиром. Каркасы условные и могут пересекаться при сжатии элементов. Силы воздействия одних элементов ВЭУ на другие вычисляются как функции их относительных перемещений и скоростей. Подбором упругости связей и размеров оснований можно моделировать все упругие свойства мачты. К самому верхнему элементу крепится ротор ВЭУ. Это обеспечивается введением стыковочных точек на валу ротора и продольной оси верхнего элемента (рис. 4). Самый нижний элемент стыкуется с фундаментом. Такая модель мачты позволяет исследовать совместные упругие колебания системы «ротор-мачта-оттяжки» при воздействии на нее аэродинамических, инерционных, упругих сил и моментов.
Рис. 4. Модель мачты Fig. 4. Mast model
International Scientific Journal for Alternative Energy and Ecology № 5 (85) 2010
© Scientific Technical Centre «TATA», 2010
Некоторые результаты моделирования
Для демонстрации возможностей описанной программно-математической модели на рис. 5-9 приведены некоторые результаты численного моделирования работы ВЭУ. Моделирование проводилось в предположении, что в роли нагрузки выступает накопитель энергии, готовый постоянно принимать любую вырабатываемую мощность. Этот режим работы ВЭУ представляет большой интерес. Действительно, он, во-первых, исключает работу вхолостую, когда ресурс ВЭУ вырабатывается без пользы, а во-вторых, исключает простои ВЭУ, увеличивающие сроки окупаемости.
Система управления задает величину снимаемой с генератора мощности Рс исходя из текущих значе-
ний скорости ветра V и скорости вращения ротора ю по алгоритму, который обеспечивает среднюю производительность ВЭУ, близкую к максимальной. Расчетная зависимость скорости ветра с учетом случайных пульсаций показана на рис. 5 (нижняя кривая). Реализовавшиеся в этих условиях скорость вращения ротора и мощность, снимаемая с генератора, показаны на рис. 5 и 6.
Из рис. 5 (верхняя кривая) видно, что скорость вращения ротора временами попадает в диапазон ю = 16-18. В эти моменты в работу включается центробежный ограничитель скорости вращения. В рассматриваемом варианте модели ВЭУ центробежный ограничитель разворачивает лопасти верхней секции ротора в тормозящее положение. Зависимость угла поворота от времени показана на рис. 7.
200 400 600 300 1000 1200
1600
20 00
2400
2600
3200 3600
Время t, с
Рис. 5. Зависимости скорости вращения ротора (верхняя кривая) и скорости ветра (нижняя кривая) от времени Fig. 5. Rotor speed (the top curve) and wind speed (the bottom curve) depending on time
5000
f 100 0
0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600
Время t, с
Рис. 6. Мощность генератора в зависимости от времени Fig. 6. Power of wind turbine depending on time
Время t, с
Рис. 7. Угол отклонения лопасти в зависимости от времени Fig. 7. Blade corner of turn depending on time
Международный научный журнал «Альтернативная энергетика и экология» № 5 (85) 2010 © Научно-технический центр «TATA», 2010
40 00
_ 30 00
li
0 250 500 750 1000 1250 1500 17 50 2000 2250 2500 2750 3000 3250 3500
Время t, с
Рис. 8. Зависимость аэродинамической силы FX от времени Fig. 8. Aerodynamic force FX depending on time
500 25 0
-250 -500 -750 -1000
250 500 750 1000 12S0 1 500 1750 2000 2250 2500 2750 3000 3250 3500
Время t, с
Рис. 9. Зависимость аэродинамической силы FY от времени Fig. 9. Aerodynamic force FY depending on time
На рис. 8, 9 приведены зависимости аэродинамических сил, действующих на ротор ВЭУ в направлении ветра ¥х и в поперечном направлении ¥Т.
Заключение
Для корректного получения представленных результатов потребовалось совместное моделирование переменного ветра, аэродинамического воздействия ветра на ротор, работы системы управления, закона управления отбором мощности, работы центробежного ограничителя скорости вращения.
Возможности представленной программно-математической модели не ограничиваются приведенным примером. Существующий набор блоков модели обеспечивает достижение всех целей моделирования, перечисленных выше.
Список литературы
1. Попов Н.А. «Рекомендации по уточненному динамическому расчету зданий и сооружений на действие пульсационной составляющей ветровой нагрузки». Госстрой России, Государственное унитарное предприятие «Центральный научно-исследовательский институт строительных конструкций имени В.А. Кучеренко ООО Еврософт». Утверждены Научно-техническим Советом ЦНИИСК 15 декабря 1999 г. М., 2000.
2. Астапов Ю.М., Медведев В.С. Статистическая теория систем автоматического регулирования и управления. М.: Наука, Главная редакция физ.-мат. литературы, 1982.
Г'-": — TATA — LXJ
International Scientific Journal for Alternative Energy and Ecology № 5 (85) 2010
© Scientific Technical Centre «TATA», 2010
