Разработка и особенности математической модели ветротурбины Дарье

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

Разработка и особенности математической модели ветротурбины Дарье

Д.В. Беляков

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

Ключевые слова - ветротурбина Дарье, режим авторотации, осреднение

I. Введение

Первые ветряные электростанции, которые могли быть эффективно использованы для получения электроэнергии, появились в конце 19 века, хотя человечество использовало энергию ветра уже сотнями лет. Их можно разделить на две большие группы: с горизонтальной и вертикальной осью вращения. Несмотря на то, что самыми распространенными являются на сегодняшний день горизонтально расположенные ветряные роторы в форме несущего винта, вертикальный ветряк остается самой первой моделью, которая использовалась когда-либо для производства энергии. Недостаток всем известных горизонтальных моделей состоит в том, что постоянно должны подстраиваться под часто меняющий направления ветер. Это означает, что самые лучшие горизонтальные ветряки большую часть времени отводят на идеальное выравнивание по ветру, тем самым они работают не в полную мощность. В 1931 году французский конструктор Жорж Дарье предложил свой вариант ротора, который имеет от двух и более плоских лопастей. Эта конструкция не зависит от направления ветра. Математическая модель ветротурбины Дарье будет рассмотрена в настоящей статье.

Рис. 1 Ветротурбина Дарье

Статья получена 4 августа 2015.

Д. В. Беляков, к.т.н., доцент, МАТИ-РГТУ им. К.Э. Циолковского, г. Москва

II. Постановка задачи

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

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

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

сопротивления SA, SB, направленные против скоростей центров давления VA, Vb , пластинок относительно среды

и подъемные силы PA, Рв, направленные им

ортогонально. При этом величины аэродинамических сил равны:

| SA |= s(a + S)VA2 = 0.5 росх (а + S)Va2 , | Р

| SB |= s(e + S)V2 = 0.5 росх (в + S)V2, |

Где а, в - углы атаки между векторами окружных скоростей ( Vao , Vbo ) точек A, B и векторами Va , Vb , p, s - аэродинамические функции углов атаки, сх , су -безразмерные аэродинамические функции, р - плотность воздуха, О - площадь одной пластинки. Зависимости Сх (а) и су (а) определены из продувок прямоугольных

пластинок с заданным удлинением в аэродинамической трубе и являются экспериментальными данными. Их типичный вид (для удлинения Л = 8) представлен на рисунке 3.

В качестве обобщенных координат, определяющих положение тела, введем угол в отклонения стержня АВ от оси x. Для описания распределения скоростей точек

1

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

нашего тела, зададим величину вектора скорости ветра V и абсолютную угловую скорость стержня o.

В системе координат, связанной с вектором V (см. рис. 6), для векторов VAO ,VA ,VB имеем представление:

VAO = {-rosin 0, rocos 0},

VA = {V - rosin 0, rocos 0},

VB = {V + rosin 0, - rocos 0}

Следовательно, силы сопротивления можно определить как:

SA = s(a + S)VA {rosin 0- V, - rocos 0}, SB

А подъемные силы будут иметь вид:

PA = p(a + S)VA {rocos 0,

типичный вид (для удлинения Я = 8) представлен на рисунке 3.

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

V и абсолютную угловую скорость стержня o.

В системе координат, связанной с вектором V (см. рис. 2), для векторов VAO ,VA ,VB имеем представление: VAO = {-rosin 0, rocos 0}, VA = {V - rosin 0, rocos 0}

= s(P + S)^-V -osfh 0 ocosocos 0}

Следовательно, силы сопротивления можно определить как:

SA = s(a+S)VA{rmsin0~V, - rocos 0},

rosin 0—v}, pB = p{p + -Х-0o8ffl0

Рис. 2 Рассматриваемое тело Где а, в - углы атаки между векторами окружных скоростей (VAO, VBO) точек A, B и векторами V, VB, p, s - аэродинамические функции углов атаки, Сх , су -безразмерные аэродинамические функции, р - плотность воздуха, 7 - площадь одной пластинки. Зависимости cx (а) и су (а) определены из продувок прямоугольных

пластинок с заданным удлинением в аэродинамической трубе и являются экспериментальными данными [6]. Их

а

Рис. 3 Аэродинамические функции: а) подъемной силы; б )силы сопротивления

подъемные силы будут иметь вид:

PA = p(a+S)VA{rocos0, rosin0-V},

PB = p(e + S)VB{-rocos0, - V - rosin 0}

Вычислим момент аэродинамических сил, приложенных к пластинкам относительно центра масс.

M0 = r(V2(p(a+S)sina-s (a+S)cosa) + VB2(p(P + S)sinв-s(P + S)cosp)ez

J0=rtyAPa+S sina-s (a+S)oosO)+

+V?

Соотношения, определяющие величины VA,VB, a,в , через V ,0,o, имеют вид:

2

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

VA sin a = -V cos в, VA cos a = ro-V sin в A A (2) VB sinв = Vcose, VB cosв = ro+Vsine

После интегрирования системы уравнений (1) -(2), мы можем окончательно определить положение тела.

Таким образом, уравнения (1)-(2) представляют замкнутую систему, полностью описывающую движение рассматриваемого тела в рамках принятых предположений.

III. ПРОСТЕЙШИЕ СТАЦИОНАРНЫЕ РЕЖИМЫ

Будем искать простейшие установившиеся режимы движения, при которых у тела отсутствует вращение. При этом 0 = 0. Правая часть уравнения изменения

кинетического момента системы (1.1.1) обращается в нуль, так что

[p(a+8) sin a- s(a+8)cosa+

+р(в + 8)sin в - s(e + 8) cos в] = 0

Действительно, аэродинамические функции s(a), p(a) имеют период П, а углы атаки при поступательном движении связаны соотношением: в = П + a.

Проводя простейшие преобразования, получим:

Рис. 4 Первое положение равновесия

Рис. 5 Второе положение равновесия

[ p(a + 8)sina- s(a + 8)cosa + p(e + 8)sin в- s(в + 8)cos в] =

= [ p(a + 8)sin a - s(a + 8)cos a + p(n + a+8) sin(n + a) - s(n + a + 8) cos(n + a)] =

= [ p(a + 8)sina- s(a + 8)cosa- p(a+ 8)sina + s(a+ 8)cosa] = 0

Таким образом, в рассматриваемом случае, суммарный момент аэродинамических сил равен нулю для любого положения тела e = e*

При отсутствии вращения из кинематических соотношений (1.1.2) следует: tga = ctge, откуда

IV. СТАЦИОНАРНЫЙ РЕЖИМ АВТОРОТАЦИИ

Рассмотрим движение тела с высокой угловой скоростью. При этом окружные скорости пластинок

VAo, VBO велики по сравнению с абсолютной

скоростью центра масс V , поэтому углы атаки a, в

„ П a+e = —. 2

ctg (a) =

p

cy

= k

Следует отметить два наиболее характерных положения равновесия (см. рис.4, 5):

близки к нулю. Из уравнения изменения кинетического момента (1.1.1) видно, что его правая часть достаточно мала. Предположим, что существует стационарный режим, при быстром вращении ветротурбины. Итак, пусть турбина раскручена до высокой угловой скорости, т.е.

I VAO 1>>1 V l,a,в ^ 0.

Аэродинамические функции p, s разложим в ряд,

учитывая, что a, в ^ 0

s(a) = 2po(cX8 + c8a), р(a) = 2pa(cy8 + c'y8a),

s

,где:

S8 = s(8), s8 = s'(8), p8 = p(8), p8 = P (8),

сх8 = сх8

Левые части кинематических соотношений (2) линеаризуем при a, в ^ 0

3

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

VAa = -V cos в VA = ra-V sine

vBe = v cose

VB = ra+ V sine

Пользуясь этими соотношениями, сделаем в уравнениях (1) , переход от переменных a,j3,VAVB к переменным

e,V ,о

После вычислений получим уравнения:

JO = pjr{-cxS(ra)2 + (c'yS cos2 в- cxS sin2 в+

+(cyS-c'xs )sinecose)V2}

Выделив в правых частях уравнения (3.1.1) периодическую часть, среднее значение которой равно нулю, и стационарную часть, получим уравнение:

1

Ja =pjr{-cxS (a)2 + - V 2(c'yS- cxS) +

1

+ 2 V2[(c'y, + cxS)cos 2в + (c yS-c'xS )sin 2в)] }

Проведем осреднение правых частей уравнения (6) на отрезке [0,2п], оставив в ней только стационарные слагаемые и учтем, что среднее значение Y на обороте равно нулю:

1

Je = pjr{-cxS

Осредненное уравнение (3.1.3) имеет стационарное решение О0:

о =-

V

r

cyS С xS

2с„

xS В случае малых значений S

справедливы соотношения:

cyS = c'y0S c'yS = c'x0 , cxS = cx0 + cx2S' . cxS = 2cx2S  При малых значениях S будем иметь:

/ rO

V _

r^l 2

a0 =-

Режим (9) называется режимом авторотации.

Можно сделать следующие выводы:

1. Для фиксированных пластинок, величина

стационарного значения (00 зависит от параметров: S, r, V и аэродинамических коэффициентов cx0, cx2.

2. Чем меньше значения аэродинамических функций cx0, cx2, тем больше (й0. 3. Чем больше установочный

угол S пластинок, тем меньше угловая скорость ветротурбины.

V. УРАВНЕНИЯ ПЕРВОГО ПРИБЛИЖЕНИЯ РЕЖИМА АВТОРОТАЦИИ

Проведем исследование устойчивости установившегося решения (9). Введем малые отклонения: (О = (й0 + X и

линеаризуем осредненное уравнение (7). Уравнение первого приближения режима авторотации будет иметь в ид :

jX = -2sxSr2 X

Корень характеристического уравнения отрицателен и режим авторотации асимптотически устойчив.

VI. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ

ТЕЛА

Полученные ранее основные аналитические

результаты использованы для построения с помощью математического пакета MATLAB 6.5 программы, реализующей численное интегрирование уравнений движения тела (приложение программа 1). При поиске численного решения используется процедура ode45, реализующая методы Рунге-Кутта четвертого и пятого порядка с автом атич ес ким в ыбор о м ш ага :

При запуске программы вводятся числовые значения параметров : j , r , J , S : splast=0.32; a2=splast/8; a=sqrt(a2); r=8*a; r0=1.27;G=0.05; del=3;V Значения аэродинамических функций для

прямоугольных пластинок с удлинением Л = 8 вводятся следующим образом: alfa=- 0.5 * pi : ( pi* 5 )/ 180:0.5*pi;

s=[ 1 .4 1 .3 8 1 .3 5 1 .3 1 .25 1.17 1.09 0.99 0.89 0.79 0.71 0.61 0.2 0.28 0.4 0.5 0.61 0.71 0.79 0.89 0.99 1.09 1.17 1.25 1.3

p=[0 -0.17 -0.28 -0.4 -0.5 -0.6 -0.65 -0.72 -0.77 -0.81 -0.8' 0.4 0.78 0.75 0.73 0.79 0.87 0.9 0.87 0.81 0.77 0.720.65 0

После этого вычисляются стационарные значения О0,

полученные при помощи метода осреднения:

d1p0=(p(:,20)-p(:,19))/hh

d2s0=

sdel=s(:,19)+0.5*d2s0*delA2

d1sdel=d2s0*del

pdel=d1p0*del

bb=sqrt

wst=

Далее вводится вектор начальных условий:

У0 = [в(0), в(0)] и временной интервал Tspan , на

котором проводится интегрирование:

Tspan=[0 5]; Yo=[0 10]

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

function dydt =vetroturbina1(T,Y);

% функция, вычисляющая правые части уравнений dydt=zeros(2,1); % a column vector;

4

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

global v r G s p alfa nap del dima V;

tgalfa1=(-V*cos(Y(1)))/(r*Y(2)-V*sin(Y( 1)));

tgbeta1= (V*cos(Y(1)))/(r*Y(2)+V*sin(Y(1)));

alfa 1=atan(tgalfa 1);

beta1=atan(tgbeta1);

sa=spline(alfa,s,alfa1);

pa=spline(alfa,p,alfa1);

sb=spline(alfa, s ,beta 1);

pb=spline(alfa,p,beta1);

VA2=VA2+(r*Y(2))A2-2*V*r*Y(2)*sin(Y(1));

VB 2=VA2+(r*Y(2))A2+2*V*r*Y(2) *sin(Y(1)); dydt(1)=Y(2);

dydt(2)=nap*(r*(VA2*(pa*sin(alfa1)-

sa*cos(alfa1))+VB2*(pb*sin(beta1)-

-sb*cos(beta1))))/G;

При работе программы процедура ode45 постоянно обращается к файлу-функции vert1. Далее происходит графический вывод полученного численного решения: [T,Y]=ode45(@vert1,Tspan,Yo,a1); plot(Y(:,1),Y(:,2),'M'), grid on

После выхода тела на режим авторотации и окончания интегрирования уравнений движения тела происходит сглаживание полученного решения. Численные значения угловой скорости, скорости спуска и угла планирования аппроксимируются набором кубических сплайнов и проводится их осреднение на последнем обороте. При этом используются функции fnint и fnval из пакета расширения Spline Toolboox. Полученные численные результаты сравниваются с теоретическими оценками, найденными при помощи метода осреднения: sf=max(Y (:,1)) ;ob=fix(sf/(2*pi)) ;pp=spline(Y (:,1),Y(:,2));

I= fnint

wshet2=fnval(I,2*pi*(ob-1))-fnval(I,0);wshet=(wshet1-

wshet2)/

stats=[ wst ]; shet=[wshet ]

VII. ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ПРОС ТЕЙШИХ ПОЛОЖЕНИЙ РАВНОВЕСИЯ ТЕЛА

Проведем моделирование движения тела при отсутствии вращения, имеющего пластинки с удлинением Л = 8. Зададим размеры тела, установочный угол пластинок, момент инерции, начальные условия и проведем численное исследование движения тела около

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

V = 10 м/c., r = 1.6 м., J=0.01kf м2, 7 = 0,32м2,

8 = 00. Выберем начальные условия, близкие к стационарным значениям: [#(0), #(0)] = [0,3]. В

качестве результатов моделирования сформируем фазовую зависимость угловой скорости (рис. 7).

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

Проведем теперь моделирование движения тела в окрестности простейшего режима вертикального спуска (рис. 5). При создании модели, возьмем следующие числовые значения параметров: V = 10 м/c.,

r = 1.6 м., J=0.01kf  м2, 7 = 0,32 м2, 8 = 00.

Выберем начальные условия, близкие к стационарным

• п

значениям: [#(0),#(0)] = [-,1.95]. В качестве

результатов моделирования сформируем фазовую зависимость угловой скорости (см. рис. 7).

Рис. 7

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

VIII. ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ВРАЩЕНИЯ ТУРБИНЫ ВБЛИЗИ РЕЖИМА АВТОРОТАЦИИ

Проведем моделирование движения турбины на режиме авторотации. Приведем некоторые результаты численного интегрирования для тела заданных размеров, имеющего пластинки с удлинением Л=8 при различных значениях массы, момента инерции, площади пластинок, перекоса пластинок. При построении модели, используем следующие числовые значения параметров:

М = 10 м / c, J=0.1ktQm2, 7 = 0,32м2, r = 1.6 м.

5

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

Начальные условия будут иметь вид: [#(0), (9(0)] = [0,70]. В качестве результата

моделирования изобразим на рисунках 9-14 фазовые зависимости угловой скорости, скорости центра масс и траекторию центра масс при установочных углах

пластинок: 8 = {0°,1°, 2°,3°, 40,50}.

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

Рис. 10

Рис. 11

8

О

Рад.

V сек. у Аналитические значения

94.51420157266902

93.58095831281183

90.93837860168216

86.99025210302818

82.23526992959449

77.12832888080796

О

Рад.

V сек. у

Расчетные значения

93.01331674391736

87.12227849360717

68.08085059501784

43.85354503174106

25.57964325550195

Вращения нет

Таблица 1

0

0

0

1

0

2

0

3

0

4

0

5

Рис. 8

Рис. 9

Численное решение, изображенное для случая 8 = 0° позволяет сделать следующие выводы:

6

International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 9, 2015

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

Во-вторых, численное решение притягивается к

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

В-третьих, сравнивая численные и стационарные значения, мы убеждаемся в том, что стационарные значения угловой скорости, при достаточно малых значениях S отличаются от соответствующих численных значений на величину порядка малого параметра, что соответствует основным выводам из теоремы Боголюбова.

В-четвертых, увеличение момента инерции приводит к увеличению времени выхода на стационарный режим.

По виду численного решения, изображенного на

рисунках, изображенных для случаев S = 10 — 5° позволяет сделать следующие выводы:

Во-первых, при увеличении перекоса пластинок

уменьшается угловая скорость, что вытекает также из формул первого приближения метода осреднения. Во-вторых, чем больше величина перекоса пластинки, тем больше расхождение между теоретическими оценками и результатами численного интегрирования (рис 43). При S = 0 результаты численного счета достаточно близки к аналитическим оценкам. Это можно объяснить тем, что формулы режима авторотации, полученные при помощи метода осреднения, имеют смысл только при достаточно малых значениях S . В-третьих, когда величина S становится больше пяти градусов, авторотация тела практически прекращается.

IX. ЗАКЛЮЧЕНИЕ Перечислим выводы работы:

1. Разработана математическая модель движения ветротурбины.

2. Указаны некоторые положения равновесия.

3. Показано, что в широком диапазоне параметров существует стационарный режим авторотации, который асимптотически устойчив.

4. Установлена возможность управления углом отклонения пластинок скоростью вращения турбины.

5. На базе специализированной системы компьютерной математики MATLAB 6.5 разработана программа для имитационного моделирования движения тела.

сопротивляющейся среде». Издательство «МЭИ», журнал «Вестник МЭИ», выпуск № 4 2006 г., стр. 5-10.

[2] Беляков Д.В. "Исследование и особенности

математической модели движения несимметричного авторотирующего тела в квазистатической среде". Издательство «Новые технологии», журнал

"Мехатроника, Автоматизация, Управление". Выпуск № 11. 2007 г., стр. 20-24

[3] Локшин Б.Я. , Привалов В.А., Самсонов В.А. “Введение в задачу о движении точки и тела в сопротивляющейся среде”. Издательство Московского университета. 1992.

[4] Паршин Д.Е. Самсонов В.А. «Качественный анализ в задаче о движении аэродинамического маятника.»1992г. (МГУ НИИмеханики, отчет 419).

[5] Б.Я. Локшин, В.А. Привалов, В.А. Самсонов «Качественный анализ движения твердого тела в сопротивляющейся среде».Шестой Всесоюзный съезд по теоретической и прикладной механике. Ташкент, 1986 г.

[6] Табачников В.Г. «Стационарные характеристики

крыльев на малых скоростях во всем диапазоне углов атаки.» Труды ЦАГИ 1974 г. выпуск 1621

[7] Strickland J.N. , Smith T. and Sun K. ,”A vortex Model of the Darrious Turbine: An analitycal and experimental stady”, Final Report submitted to Sandia Laboratories on contract # 13-5602, 1982

[8] Г.И. Шоломович, В.П. Карликов, О.В. Зверев «Исследование по оптимизации ортогональных

энергоагрегатов для использования энергии воздушных и водных течений. Оценка влияния кривизны и неоднородности потока на на аэродинамические

характеристики профиля»,-М.,1990 г.-33 л.(МГУ.

НИИмеханики, отчет 3910)

[9] Paraschivoiu J., “Double Multiple Stremeamtube model with Recent Improvements”, Journal of Energy, vol.7 no.3

[10] В.Ф. Журавлев Д.М. Климов «Прикладные методы в теории колебаний» Издательство «Наука» 1988 г.

[11] Самсонов В.А. Беляков Д.В. «Оптимизация

движения спускающегося авторотируюшего тела». Седьмой Международный Аэрокосмический Конгресс IAC'07. Полные доклады. 27-31 августа 2012 г., Москва, Россия. Электронный вид. Регистрационный номер: ISBN 7-85312-056-7

[12] Имитационное моделирование спуска

осесимметричного авторотирующего тела в квазистатической среде. Девятая Международная научнопрактическая конференция «Современные

информационные технологии и ИТ образование» 13-15 ноября 2014 г. Сборник трудов стр. 856-864

[13] Беляков Д. В. Оптимизация спуска осесимметричного авторотирующего тела в воздушной среде //International Journal of Open Information Technologies. - 2015. - Т. 3. -№. 2. - С. 11-16.

[14] Беляков Д. В. Математическое моделирование движения осесимметричного авторотирующего тела в квазистатической среде //International Journal of Open Information Technologies. - 2015. - Т. 3. - №. 3. - С. 7-16.

Библиография

[1] Беляков Д.В., Самсонов В.А., Филиппов В.В. «Исследование движения несимметричного тела в

7

Development and characteristics of the mathematical model of the

Darrieus wind turbine

D.V. Belyakov

Abstract - This paper is devoted to the questions of the mathematical modeling Darrieus wind turbine in the quasi-static environment. We show the existence of the auto-rotations mode. It lets us use the system as a wind generator. Our paper offers a set of program for estimation the adequacy of the proposed model.

Keywords - Darrieus wind turbine, auto-rotation, averaging

8