CC BY
110
УДК 519.6:531.13+533.6.011:620.91
ВЕРИФИКАЦИЯ КОМПЛЕКСА МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ АЭРОДИНАМИКИ И ДИНАМИКИ ДВИЖЕНИЯ РОТОРА САВОНИУСА
©2013 Д. А. Сизов1, Ю. П. Онушкин1, О. А. Краснова1, О. Т. Джанибеков2
1Сызранский филиал Самарского государственного технического университета 2Казанский национальный исследовательский технический университет
им. А. Н. Туполева
Представлены основные этапы верификации комплекса математических моделей аэродинамики и динамики движения ротора Савониуса, основанных на методе дискретных вихрей. Результаты расчётов приводятся в сравнении с данными численных и натурных экспериментов, проведённых другими авторами.
Математическая модель, ротор Савониуса, метод дискретных вихрей, динамика, коэффициент использования энергии ветра.
Ветроэнергетика является одним из приоритетных направлений развития отечественной энергетики. Ветроэнергетические установки (ВЭУ) можно разделить на две основные группы: с горизонтальной и вертикальной осями вращения ротора. Для небольших хозяйств и бытовых нужд выгодно использовать вертикально-осевые установки, в частности ВЭУ с ротором Савониуса, т. к. они имеют простую конструкцию, не требуют ориентации на ветер и обладают большим начальным моментом [1]. В связи с большим распространением таких ВЭУ на первый план выходят следующие задачи: задача создания каталога типоразмеров ВЭУ с заданной геометрией роторов, позволяющего определить зависимость выходных (мощностных) характеристик от скорости ветра, диаметра ротора, удлинения лопастей; задача оптимизации формы лопастей с целью увеличения коэффициента использования энергии ветра и многие другие. Самым распространённым на данный момент методом их решения является продувка опытных образцов роторов в аэродинамических трубах. Однако существенно более дешёвым и гибким средством является математическое моделирование.
В данной статье рассматривается верификация созданного авторами про-
граммного продукта, представляющего собой комплекс математических моделей аэродинамики и динамики движения ротора Савониуса, подробно описанного в работах [2, 3]. В его основе лежит метод дискретных вихрей [4, 5], позволяющий смоделировать сложный нелинейный нестационарный процесс взаимодействия лопастей ротора ВЭУ с воздушным потоком. С помощью данного метода можно провести так называемый виртуальный эксперимент, вычислить скорость потока в произвольной точке пространства при обтекании объекта, получить безразмерные коэффициенты давлений на его поверхности, мгновенную картину вихревых структур и многое другое. Для этого необходимо представить объект (ротор ВЭУ) как дискретную систему панелей прямоугольной формы. Далее, переходя от безразмерных коэффициентов давлений к размерным нагрузкам, можно вычислить крутящий момент на роторе в произвольный момент времени и смоделировать его движение в динамике. Меняя при этом геометрию ротора, скорость набегающего потока и момент сопротивления генератора, можно получить широкий спектр результатов, анализ которых поможет решению поставленных выше задач.
Прежде чем проводить виртуальный эксперимент, необходимо провести поэтапную верификацию созданных моделей и убедиться в достаточной точности вычисления:
а)координат дискретных вихрей следа за обтекаемым объектом;
б) скоростей потока в произвольной точке пространства вблизи объекта;
в) безразмерных коэффициентов давлений на панелях разбиения поверхности объекта;
г) крутящего момента на вращающемся роторе, а также его угловых ускорения, скорости и координаты;
д) осевого момента инерции ротора;
е) энергетических характеристик нагруженного ротора (в первую очередь, коэффициента использования энергии ветра).
Каждый элемент представленного перечня соответствует отдельному этапу верификации.
Проверку правильности выстраивания вихревого следа производим на примере обтекания неподвижной пластины. На рис. 1 рассчитанный авторами вихревой след за квадратной пластиной, расположенной перпендикулярно набегающему с постоянной скоростью потоку, показан в сравнении с результатами, которые приводят С. М. Белоцерковский и М. И. Ништ
[4].
Рис. 1. Рассчитанный вихревой след за неподвижной плоской пластиной в различные моменты безразмерного времени т: т = 1 (слева), т = 2,5 (справа); а — авторы, б — работа [4]
На рис. 1, а особенно хорошо заметна указанная в работе [4] тенденция сворачивания вихревых пелен в спиралевидные жгуты. При этом те части пелены, что сходят с угловых точек, движутся быстрее, чем части, образующиеся вблизи середин кромок пластины.
Далее убедимся в правильности вычисления скоростей потока в произвольных точках пространства. На рис. 2 представлено рассчитанное поле скоростей потока в плоскости, перпендикулярной пластине. Хорошо заметны характерные зоны, которые можно наблюдать в простейшем опыте на движение пластины в гидролотке: зона торможения набегающего потока перед пластиной, зона обратного течения за ней, зоны увеличения скорости потока вблизи кромок.
Смоделировано обтекание кругового цилиндра (тонкой пластины, свёрнутой в трубу с малым зазором между смежными кромками) большого удлинения (А = 200) потоком, имеющим единичную скорость и направленным вдоль оси абсцисс (рис. 3). Из теории идеальной жидкости известно, что в точках (А и В) пересечения окружности основания цилиндра с осью абсцисс поток полностью тормозится, а в точках (С и D) пересечения окружности с осью ординат его скорость принимает максимальное значение, не зависящее от радиуса цилиндра и равное удвоенной скорости на бесконечности [6]. Результаты расчёта поля скоростей в среднем сечении цилиндра демонстрируют совпадение с данными теории.
Рис. 2. Поле скоростей потока при обтекании прямоугольной пластины (а = 90°)
Рис. 3. Поле скоростей в среднем сечении цилиндра (тонкостенной трубы)
150
Рис. 4. Расчётные и экспериментальные данные по аэродинамическим нагрузкам на прямоугольной пластине: а — авторы, б — работа [4] (1 — линейная теория, 2 — нелинейная теория, точки — эксперимент)
Обратим внимание на то, что рассчитанная скорость потока в точках, лежащих внутри трубы, стремится к нулю. Это означает, что в данном случае замена сплошного цилиндра на тонкостенную незамкнутую трубу вполне допустима и не противоречит физике процесса.
Правильность расчёта безразмерных коэффициентов давления проверим для случаев обтекания плоской пластины, цилиндра и профиля.
На рис. 4 показаны результаты произведённых расчётов безразмерных коэффициентов давления Ср при обтекании неподвижной прямоугольной пластины с удлинением X = 0,25 под углом атаки а = 20° в сравнении с данными, которые приведены в [4].
Далее рассмотрим распределение безразмерных давлений на поверхности цилиндра (тонкостенной незамкнутой трубы), определяемых формулой
р=а (1)
Ч¥
где р — размерное давление (в паскалях),
р и¥
2
— динамическое давление
набегающего потока, рв — плотность воздуха, — скорость потока. Сравним теоретические значения безразмерного давления, рассчитанные с помощью теории идеальной жидкости, с результатами, полученными посредством численного моделирования (рис. 5).
а
б
Рис. 5. Распределение безразмерного давления на поверхности круглого цилиндра (среднее сечение)
Рис. 6. Визуализация обтекания крыла (а = 10°)
Рис. 7. Зависимость коэффициента подъемной силы крыта от угла атаки
Смоделировано обтекание крыла бесконечного удлинения (в плоской постановке задачи) с профилем NACA-23011 под углом атаки а потоком, направленным вдоль оси абсцисс (рис. 6).
Полученная по результатам моделирования зависимость
коэффициента подъёмной силы крыла от угла атаки (в градусах) показана на рис. 7 в сравнении с экспериментальными данными [8]. Так как периметр профиля был разбит на конечное число панелей, коэффициент подъёмной силы
определялся по формуле
СУ = Х Ср}г С08(П ,У),
(2)
г=1
где п — число панелей разбиения; 1 — безразмерная длина панели с номером г;
(П, у) — угол между внешней нормалью
к панели с номером г и осью ординат.
Расчётные значения коэффициента получились завышенными, так как моделировалось крыло бесконечного удлинения (не учитывался скос потока, присущий крылу конечного удлинения), в то время как экспериментальные данные получены для крыла с удлинением X = 5.
Следующий этап - проверка правильности вычисления крутящего момента и кинематических параметров ротора -требует моделирования обтекания не статичных, как ранее, а движущихся тел с наперёд не известным законом движения. Простейшим примером такого движения является вращение флюгера — тонкой плоской прямоугольной пластины, которая может свободно поворачиваться во-
круг оси, проходящей вдоль одной из её сторон, под действием набегающего потока. При этом пластина будет совершать затухающие колебания, так как при образовании вихревого следа происходят потери энергии.
Рассмотрим свободное вращение квадратной пластины из дюралюминия (плотность р = 2650 кг/м ) со стороной а = 1 м и толщиной Ь = 3 мм при постоянной скорости набегающего потока = 10 м/с и начальной угловой координате ф = 90°. На каждом расчётном шаге вычислялся крутящий момент, после чего производилось двукратное интегрирование дифференциального уравнения вращательного движения. Результаты интегрирования представлены на рис.8.
С увеличением числа панелей разбиения точность вычисления осевого момента инерции увеличивается. Для рассмотренной выше тонкой пластины полагаем достаточным разбиение 30 х 30 панелей: в этом случае полученное численным методом значение осевого момента инерции равно 2,52 кгм , что отличается менее чем на 5% от истинного значения 2,65 кг м , определенного по формуле 1
(4)
J =— та2, г 3
где т — масса пластины, а — длина стороны.
Реализованный алгоритм численного нахождения момента инерции особенно полезен при рассмотрении роторов более сложной геометрии. Примером может служить ротор Савониуса с профилем лопасти, представляющим собой полуокружность радиусом г (рис. 9).
Рис. 8. Зависимость угловой координаты флюгера от времени
Из рис.8 следует, что принятая модель пластины демонстрирует затухающие колебания пластины, поскольку очевидно её стремление к самоустановке параллельно набегающему потоку.
Составление дифференциального уравнения динамики вращательного движения и, как следствие, вычисление угловых скорости и координаты ротора невозможно без знания величины осевого момента инерции ротора относительно оси вращения. В рассматриваемом программном продукте осевой момент инерции ротора произвольной формы находится по стандартной формуле
п
Л =Х тк\ (3)
к=1
где п — число панелей разбиения; тк — масса панели с номером к; Ик — расстояние от центра тяжести панели с номером к до оси вращения.
Рис. 9. Двухлопастной ротор Савониуса с разбиением на панели
Аналитическое выражение для осевого момента инерции рассматриваемого двухлопастного ротора имеет вид
1 2
J, =- МЯ2,
2
(5)
где М — масса ротора; Я = 2г — радиус ротора.
Значение осевого момента инерции, рассчитанное по формуле (5) для ротора указанной конфигурации, изготовленного из листа толщиной Ь = 1 мм (сплав Д16Т, плотность 2775 кг/м ), высотой к = 1 ми радиусом Я = 0,3 м, равно 0,11768 кг м2. Найденная численным интегрированием
величина (при 40 панелях разбиения вдоль оси вращения и 100 панелях вдоль профиля) составляет 0,11357 кгм2. Погрешность не превышает 5%.
Таким образом, задавая в программном продукте произвольные конфигурации роторов, можно получать значения осевого момента инерции с достаточной точностью, не прибегая к зачастую громоздкому выводу аналитических зависимостей для его определения.
Заключительный этап верификации - определение и контроль правильности вычисления мощностных характеристик ВЭУ, важнейшей из которых является коэффициент использования энергии ветра £ (часто выражается в процентах):
X = ' (6)
р аИЯО ¥
где N — мощность ВЭУ; И — высота ротора.
Как правило, в литературе (например в [1]) коэффициент использования энергии ветра для конкретного типа ВЭУ приводится в зависимости от быстроходности её ротора Ш, определяемой по формуле
_ ШЯ ш =
различных значениях диаметра ротора и скорости набегающего потока (табл. 1).
О ¥
(7)
где а — угловая скорость ротора.
Коэффициент использования
энергии ветра и быстроходность зависят от величины момента сопротивления генератора Мс, которым нагружается ротор. Сравнение результатов
проведённых численных экспериментов для ротора рассмотренной выше конфигурации (удлинение лопастей X = И/2Я=1,67) при постоянной скорости набегающего потока, равной 10 м/с, и при различных значениях Мс (0; 0,75; 1,5; 2,25; 3,0; 3,75; 4,5; 5,0; 5,5 Н м) с данными натурного эксперимента в аэродинамической трубе КНИТУ-КАИ (В. В. Жерехов и др.) приведено на рис.10.
Используя максимальный
коэффициент использования энергии ветра, можно определить мощность ВЭУ с лопастями данного удлинения при
Рис. 10. Зависимость коэффициента использования энергии ветра от быстроходности для двухлопастного ротора Савониуса
По результатам верификации предлагаемого комплекса математических моделей можно сделать вывод, что программный продукт может быть использован для решения широкого спектра задач, стоящих перед ветроэнергетикой:
1) создание каталога типоразмеров ВЭУ с заданной геометрией роторов, позволяющего определить зависимость выходных (мощностных) характеристик от скорости ветра, диаметра ротора, удлинения лопастей;
2) определение оптимальной с точки зрения эффективности использования энергии ветра геометрии лопастей роторов (имеющих форму дуги окружности, части лемнискаты, дуги эллипса и т. п.) для уже существующих конструкций ВЭУ с вертикальной осью вращения;
3) разработка новых конструктивных схем ВЭУ;
4) моделирование динамики выхода ВЭУ на рабочий режим;
5) моделирование поведения ВЭУ при непредвиденных обстоятельствах (например, при внезапном увеличении скорости ветра выше запланированного значения);
6) расчет нагрузки на элементы ВЭУ (например, на опоры валов).
Таблица 1. Мощность ВЭУ (в ваттах) с двухлопастным ротором Савониуса для двух различных значений удлинений лопастей
UM, м/с D, м 5 7,5 10 12,5 15
X = 1,67
0,3 1,73 5,85 13,86 27,07 46,78
0,6 6,93 23,39 55,44 108,28 187,11
0,9 15,59 52,63 124,74 243,63 421,00
1,2 27,72 93,56 221,76 433,13 748,44
X = 3,33
0,3 3,92 13,21 31,32 61,17 105,71
0,6 15,66 52,85 125,28 244,69 422,82
0,9 35,24 118,92 281,88 550,55 951,35
1,2 62,64 211,41 501,12 978,75 1691,28
Библиографический список
1. Ветроэнергетика [Текст] / под ред. Д. де Рензо; пер. с англ.; под ред. Я. И. Шефтера. - М.: Энергоатомиздат, 1982. - 272 с.
2. Математическая модель аэродинамики и динамики движения ротора Са-вониуса на основе метода дискретных вихрей (плоская задача) [Текст] / С. А. Михайлов, Д. А. Сизов, Ю. П. Онушкин [и др.] // Вестн. Казанского государственного технического университета им. А. Н. Туполева. - 2011. -№3. - С. 5-12.
3. Онушкин, Ю. П. Численное моделирование аэродинамики и динамики движения роторов вертикальноосевых ветроустановок с помощью метода дискретных вихрей (плоская и пространственная задачи) / Ю. П. Онушкин, Д. А. Сизов // Вестн. СамГТУ. Приложение к журналу. - 2011. - №3. - С. 19-26.
4. Белоцерковский, С.М. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью [Текст]/ С.М. Белоцерковский, М.И. Ништ. - М.: Наука, 1978. - 352 с.
5. Нелинейная теория крыла и ее приложения [Текст] / Т. О. Аубакиров, А. И. Желанников, С. М. Белоцерковский [и др.] - Алматы: Гылым, 1997. - 448 с.
6. Аржаников, Н. С. Аэродинамика [Текст] / Н. С. Аржаников, В. Н. Мальцев.
- М.: Государственное издательство оборонной промышленности, 1956. - 483 с.
7. Шлихтинг, Г. Теория пограничного слоя [Текст] / Г. Шлихтинг; пер. с нем.; под ред. Л. Г. Лойцянского. - М.: Наука, 1974. - 217 с.
8. Атлас аэродинамических характеристик профилей крыльев. [Текст] / Б.А. Ушаков, П.П. Красильщиков, А.К. Волков [и др.]. - М.: БНТ НКАП при ЦАГИ, 1940.
- 339 с.
VERIFICATION OF A COMPLEX OF MATHEMATICAL MODELS OF AERODYNAMICS AND DYNAMICS OF SAVONIUS ROTOR MOTION
© 2013 D. A. Sizov1, Yu. P. Onushkin1, O. A. Krasnova1, O. T. Dzhanibekov2
1Syzran branch of Samara State Technical University Kazan National Research Technical University
The paper presents the main stages of verification of a complex of mathematical models of aerodynamics and dynamics of a Savonius rotor motion based on the discrete vortex method. The results of calculations are compared with the data of experiments conducted by others authors.
Mathematical model, Savonius rotor, discrete vortex method, dynamics, wind power efficiency..
Информация об авторах
Сизов Дмитрий Александрович, старший преподаватель кафедры технической механики, Сызранский филиал Самарского государственного технического университета. E-mail: sizov.syzran@gmail.com. Область научных интересов: математическое моделирование методом дискретных вихрей.
Онушкин Юрий Петрович, кандидат технических наук, доцент, доцент кафедры технической механики, Сызранский филиал Самарского государственного технического университета. Область научных интересов: математическое моделирование методом дискретных вихрей.
Краснова Оксана Александровна, преподаватель кафедры технической механики, Сызранский филиал Самарского государственного технического университета. Email: ok-kras@ya.ru. Область научных интересов: математическое моделирование методом дискретных вихрей.
Джанибеков Олег Тофикович, начальник отдела, Казанский национальный исследовательский технический университет им. А.Н. Туполева. E-mail: oleg.djanibekov@yandex.ru. Область научных интересов: математическое моделирование методом дискретных вихрей.
Sizov Dmitry Alexandrovich, senior lecturer of the technical mechanics department, Syzran branch of Samara State Technical University. E-mail: sizov.syzran@gmail.com. Area of research: mathematical modeling using the discrete vortex method.
Onushkin Yury Petrovich, candidate of technical science, associate professor of the technical mechanics department, Syzran branch of Samara State Technical University. Area of research: mathematical modeling using the discrete vortex method.
Krasnova Oksana Alexandrovna, lecturer of the technical mechanics department, Syzran branch of Samara State Technical University. E-mail: ok-kras@ya.ru. Area of research: mathematical modeling using the discrete vortex method.
Dzhanibekov Oleg Tofikovich, head of department, Kazan National Research Technical University. E-mail: oleg.dj anib ekov@yandex. ru. Area of research: mathematical modeling using the discrete vortex method.
