Известия РАН. Теория и системы управления, 2021, № 1, стр. 169-176

Введение. Моделирование устройств преобразования потока (ветроэнергетических и гидроэнергетических установок) как управляемых динамических систем представляет собой перспективное направление. В качестве средств реализации управления могут выступать различные элементы системы. Широкий спектр работ описывает управление посредством изменения установочного угла лопастей [1–3]. Вероятно, это направление появилось первым в области регулирования работы ветроустановок. Известны исследования по управлению частотой вращения турбины и ротора генератора установки при помощи дифференциальной планетарной передачи [4–6]. Помимо этого возможно управление внешней электрической нагрузкой в цепи генератора в зависимости от изменения скорости ветра. Подобный подход особенно актуален для малогабаритных установок [7]. Последние характеризуются не только небольшими размерами, но и подключением генератора к локальной внешней цепи, нагрузку в которой можно (и весьма целесообразно) регулировать с целью увеличения мощности [8].

Построение замкнутой математической модели малогабаритной горизонтально-осевой ветроустановки на базе относительно простой модели аэродинамического воздействия [9] позволило получить ряд новых результатов по описанию эффекта гистерезиса [7] и построить стратегии регулирования мощности [5, 10].

В работе [11] модель ветроэнергетической системы, предложенная в [9], была модифицирована для описания динамики ветроустановки оригинальной конструкции на базе кривошипно-шатунного механизма. Полученные результаты позволили разработать и спроектировать новый тип устройства преобразования энергии потока [12]. Более того, была построена динамическая модель ветромобиля, использующего для перемещения такую ветроустановку [13]. Макет ветромобиля был собран и протестирован в НИИ механики МГУ.

Для ветроэнергетического устройства указанного нового типа актуально исследование возможностей управления внешней электрической нагрузкой с целью максимизации мощности, отбираемой у потока на рабочем режиме. Подобная задача для других ветроэнергетических систем обсуждалась в [10, 14, 15]. В настоящей работе она рассматривается для кривошипно-шатунной установки. При этом динамическая модель системы усложнена по сравнению с [11] тем, что учитывается сухое трение в ползуне механизма.

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

Для формирования рабочих режимов с заданными свойствами (определенным значением мощности, отбираемой у потока) используется итерационный численно-аналитический метод [14], который представляет собой модификацию метода Андронова–Понтрягина [16, 17], примененного для решения подобной задачи в [11]. Модифицированный подход к исследованию установившихся движений позволил (по сравнению с [11]) рассмотреть значительно более широкую область значений параметров модели. Отметим, что указанный метод предназначен для поиска и формирования режимов авторотации с заданными свойствами. Для решения аналогичных задач применительно к режимам автоколебаний разработаны модификации метода, представленные в [18, 19].

1. Постановка задачи. Рассматривается модель устройства преобразования энергии потока, в основу конструкции которого положен кривошипно-шатунный механизм (рис. 1). Другие ветроустановки, в которых используются шарнирные механизмы, рассмотрены, например, в работах [20, 21]. Кривошип OA длины r может вращаться вокруг вертикальной оси O. Шатун AB имеет длину l (l > r). В точке B находится ползун, который может двигаться вдоль горизонтальной направляющей, обозначенной на рис. 1 как ось Oy.

На шатуне установлена прямоугольная лопасть, соединенная с ним жестко. Предполагается, что хорда лопасти направлена по шатуну. Лопасть имеет характерную площадь S и взаимодействует с потоком среды плотности $\rho $, имеющим на бесконечности скорость V, направленную вдоль оси Ох, которая ортогональна оси Оy. Твердое тело “стержень AB + лопасть” характеризуется центральным моментом инерции J₂ и массой m. Центр масс этого тела совпадает с точкой А, и предполагается, что компонента нормальной реакции, компенсирующая вес тела, приложена в точке А.

На валу кривошипа установлен ротор генератора установки, J₁ – момент инерции твердого тела “кривошип OA + ротор генератора” относительно оси O. Центр масс этого тела расположен на оси О.

Система имеет одну степень свободы. В качестве обобщенной координаты выберем угол φ между осью Ох и кривошипом ОА. Пусть $\theta $ – угол между шатуном AB и осью Ох, тогда cosθ = = $(r{\text{/}}l){\text{cos}}\varphi $ (рис. 1). Отметим, что минимальное значение $\sin \theta $ во время движения системы составляет $\sqrt {1 - {{r}^{2}}{\text{/}}{{l}^{2}}} > 0$.

Кинетическая энергия описанной системы имеет вид

Воздействие потока среды на лопасть будем моделировать как сумму силы D лобового сопротивления и боковой (подъемной) силы L. Будем считать, что точка приложения аэродинамических сил совпадает с точкой В. Значения аэродинамических сил определим на основе квазистатического подхода [8, 9], их проекции на оси Ox и Oy имеют вид

Здесь α – мгновенный угол атаки, т.е. угол между воздушной скоростью точки B и отрезком AB. Воздушная скорость точки B – вектор, имеющий в осях Oxy вид $\{ - V,u\,V\} $. Функции ${{C}_{L}}(\alpha )$, ${{C}_{D}}(\alpha )$ – безразмерные коэффициенты боковой силы и силы сопротивления лопасти. Промежуточные безразмерные переменные u и α определяются по формулам $u = {{V}^{{ - 1}}}r\dot {\varphi }{\text{sin}}(\theta + \varphi ){\text{si}}{{{\text{n}}}^{{ - 1}}}\theta $, $\alpha = - \arctan (u)$ + θ.

При численных расчетах использовались зависимости коэффициентов ${{C}_{D}}(\alpha )$, ${{C}_{L}}(\alpha )$ силы сопротивления и боковой силы от мгновенного угла атаки α, соответствующие стационарным аэродинамическим характеристикам плоской пластины удлинения 8 (т.е. прямоугольной пластины, у которой отношение длин сторон равно 8, более длинная сторона перпендикулярна AB).

На ползун со стороны направляющей действует сухое трение F_тр с коэффициентом f. Будем рассматривать такие движения системы, при которых скорость точки B может быть нулевой только в отдельные моменты времени, но не на каком-либо отрезке времени (далее будет показано, что в рамках построенной модели рабочие режимы системы удовлетворяют этому свойству). Соответственно используем модель трения скольжения, трение покоя не рассматриваем.

Электромеханический момент, описывающий воздействие статора генератора на ротор, будем считать линейным по угловой скорости $\dot {\varphi }$ ротора генератора [8]:

здесь С – коэффициент электромеханического взаимодействия, σ – внутреннее сопротивление генератора (константы, определяемые характеристиками генератора); $R$ – внешнее сопротивление в цепи генератора (будем рассматривать его как управляющее воздействие, подобно работе [15]).

В безразмерных переменных и параметрах уравнения движения системы имеют следующий вид:

Здесь z – безразмерная угловая скорость; $p = k(\varphi )z$ – безразмерный импульс; $\tau = Vt{\text{/}}r$ – безразмерное время; $Q(\varphi ,z)$ – безразмерная обобщенная сила; $\mu = {{m}^{{ - 1}}}{{V}^{{ - 2}}}rf$ – безразмерный коэффициент сухого трения, действующего на ползун.

Зависимость $N(\varphi ,z)$ – безразмерная функция, характеризующая проекцию нормальной реакции в ползуне на ось Ox. Параметр c – безразмерный коэффициент полезной нагрузки. Далее будем считать параметр с коэффициентом усиления управляющего воздействия (технически управление осуществляется путем изменения величины $R$ внешнего сопротивления).

Выражение для величины реакции $N(\varphi ,z)$ может быть получено выписыванием второго закона Ньютона для всей системы, для кривошипа, а также закона изменения кинетического момента всей системы относительно точки O:

Выражение $N(\varphi ,z)$ используется для замыкания системы (1.1). Приведенная формула для $N(\varphi ,z)$ является в определенном смысле рекурсивной: требуется информация о знаке величины $(N(\varphi ,z){\text{cos}}\varphi )$, чтобы вычислить ${{N}_{2}}(\varphi ,z)$. При численной реализации определяем два различных значения $N(\varphi ,z)$, соответствующих тому или иному знаку $\left( {N(\varphi ,z)\cos \varphi } \right)$, затем проверяем знаки полученных величин $N(\varphi ,z)$, чтобы установить, какой из двух рассмотренных знаков произведения $\left( {N(\varphi ,z)\cos \varphi } \right)$ корректен. При относительно небольших значениях безразмерного коэффициента трения μ знаменатель ${{N}_{2}}(\varphi ,z)$ функции $N(\varphi ,z)$ не обращается в нуль ни при каких φ и z. При этом знак слагаемого, содержащего множитель $\mu $, не влияет на знак величины $N(\varphi ,z)$ и соответственно указанная выше простая проверка достаточна для однозначного определения знака выражения $\left( {N(\varphi ,z)\cos \varphi } \right)$.

Правая часть (1.1) – функция 2π-периодическая по φ. Орбитально устойчивая 2π-периодическая по φ траектория системы (1.1) соответствует рабочему режиму устройства. Пусть функция $\tilde {z}(\varphi )$ описывает зависимость переменной $z$ от угла $\varphi $ на 2π-периодической траектории.

Для практических приложений представляют интерес значения средней безразмерной угловой скорости и механической мощности на рабочих режимах. Они зависят от коэффициента $c$ полезной нагрузки.

Среднее (по углу за период) значение коэффициента c_p механической мощности на рабочем режиме определим следующим образом:

где $P$ – механическая мощность турбины, ${{P}_{0}} = 0.5\rho {{S}_{0}}{{V}^{3}}$ – мощность потока, проходящего через осевое сечение ветроприемного механизма, S₀ – площадь осевого сечения.

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

2. Рабочие режимы системы. Будем рассматривать траектории, расположенные в области z > 0, поскольку из [11] следует, что соответствующие режимы намного предпочтительнее с точки зрения отбора мощности у потока, чем режимы, при которых кривошип вращается “по часовой” (z < 0). С целью построения 2π-периодических по φ траекторий системы (1.1) в зависимости от параметров модели воспользуемся методом [14].

В соответствии с этим подходом рассмотрим порождающую гамильтонову систему: при ε = 0 система (1.1) имеет гамильтониан ${{H}_{0}}(\varphi ,p) = 0.5k(\varphi ){{p}^{2}}$. При достаточно малых значениях ε асимптотические оценки характеристик рабочих режимов моделируемой установки получены в [11] в предположении отсутствия трения (f = 0) путем применения метода Андронова–Понтрягина к системе (1.1) как к системе, близкой к гамильтоновой. В настоящей работе используем траектории порождающей консервативной системы с гамильтонианом H₀ в качестве нулевой итерации ${{\tilde {z}}_{0}}(\varphi )$ для построения $2\pi $-периодических траекторий, последующие приближения ${{\tilde {z}}_{n}}(\varphi )$ определим по формулам, предложенным в статье [14]. Будем считать, что итерационный процесс сошелся, если максимум (по φ) модуля разности ${{\tilde {z}}_{n}}(\varphi )$ и ${{\tilde {z}}_{{n - 1}}}(\varphi )$ меньше некоторого заданного числа. Если при этом ${{\tilde {z}}_{n}}(\varphi ) > 0$ при всех φ, то используем функцию ${{\tilde {z}}_{n}}(\varphi )$ как итоговое приближение для 2π-периодической траектории $\tilde {z}(\varphi )$.

Рассмотрим случай, соответствующий малогабаритной ветроустановке с конкретными значениями геометрических и инерционных параметров: ${{J}_{1}} = 0.00375$ кг ⋅ м², ${{J}_{2}} = 0.0375$ кг ⋅ м², r = 0.1 м, $l = 0.5$ м, $S = 0.003$ м², $m = 1.225$ кг, ${{S}_{0}} = 0.015$ м². Пусть величина скорости ветра $V = 7$ м/с; плотность воздуха $\rho = 1.125$ кг/м³. При этих значениях параметров получаем $\varepsilon = 0.5J_{1}^{{ - 1}}\rho S{{r}^{3}}$ = 0.0005.

Помимо этого рассмотрим значения параметров модели, соответствующие работе устройства в потоке воды: $\rho = 1000$ кг/м³, $V = 0.5$ м/с (остальные параметры те же, что и в первом случае). При этом значение безразмерного параметра ε = 0.4. Отметим, что конструкция установки позволяет при использовании системы в воде полностью изолировать ползун от воздействия жидкости, поместив его в защитный кожух, что является характерным (одним из основных) преимуществом установок шарнирного типа, в которых рабочий элемент совершает колебания, а не ротации. Таким образом, при рассмотрении системы как гидроэнергетической по-прежнему можно считать, что в ползуне действует именно сухое трение.

Во всех приведенных случаях знаменатель ${{N}_{2}}(\varphi ,z)$ выражения (1.2) для реакции $N(\varphi ,z)$ не обращается в нуль ни при каких значениях безразмерного коэффициента $\mu $ сухого трения.

На рис. 2, 3 представлены бифуркационные диаграммы средней безразмерной угловой скорости 〈z〉 кривошипа и среднего коэффициента c_p механической мощности устройства на режимах авторотации, отвечающих 2π-периодическим решениям системы (1.1), для случая $\varepsilon = 0.0005$ (рис. 2), а также для случая ε = 0.4 (рис. 3). Они построены при различных значениях коэффициента сухого трения f ($f = 0,0.1,0.25,0.5$). На этих диаграммах сплошные линии отвечают притягивающим периодическим решениям, пунктирные – отталкивающим. Достаточное условие орбитальной устойчивости/неустойчивости периодической траектории: отрицательное/положительное значение угла наклона касательной к бифуркационной кривой в точке бифуркационной диаграммы, соответствующей данной траектории. Это условие – следствие того, что параметр c в системе (1.1) отвечает за поворот векторного поля [17].

Проиллюстрировано изменение средней (по углу) безразмерной угловой скорости и среднего (по углу) коэффициента c_p механической мощности в зависимости от коэффициента c полезной нагрузки. Для всех построенных $2\pi $-периодических траекторий выполнено $\tilde {z}(\varphi ) > 0$ при всех $\varphi $, тем самым скорость точки B (ползуна) не остается нулевой ни на каком интервале времени.

Можно видеть, что при $\varepsilon = {\text{0}}{\text{.0005}}$ в случае  f = 0 есть два интервала значений параметра c, для которых существует по два притягивающих 2π-периодических режима. При  f = 0.1 и 0.25 такой интервал один. В связи с наличием этих интервалов наблюдается гистерезис характеристик рабочих режимов в зависимости от направления изменения параметра c. При  f = 0.5 для всех значений коэффициента $c$ существует только один 2π-периодический режим, который является притягивающим.

При $\varepsilon = 0.4$ в случаях  f = 0 и 0.1 есть два диапазона значений параметра c, в которых есть по два притягивающих 2π-периодических режима, а в случае  f = 0.25 такой диапазон один.

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

3. Управление коэффициентом полезной нагрузки. Рассмотрим возможность регулирования механической мощности путем изменения внешнего сопротивления в цепи генератора. Для конкретной установки все размерные параметры модели фиксированы, кроме величины V скорости потока, и значения R внешнего сопротивления, которые описывают текущие условия эксплуатации. При изменении скорости V можно регулировать величину R внешнего сопротивления таким образом, чтобы значение коэффициента c полезной нагрузки соответствовало максимуму коэффициента c_p мощности: $c = {{c}_{{opt}}}$. Значение безразмерного коэффициента $\mu $ трения зависит от величины V скорости потока. Соответственно при изменении V меняется значение коэффициента c, доставляющее максимум коэффициента c_p мощности (на рис. 2, 3 это отвечает переходу от одной диаграммы к другой). В связи с этим для построения зависимости величины R от скорости V потока требуется информация о функции ${{c}_{{opt}}}(\mu )$, а не об одном значении ${{c}_{{opt}}}$ (в отличие от модели, где не учитывалось трение). В остальном выбор значения ${{R}_{{opt}}}$, которое требуется для поддержания максимальной механической мощности, аналогичен предложенному в работе [15] для другой ветроустановки:

Функция ${{c}_{{opt}}}(\mu )$ строится по абсциссам максимумов диаграмм ${{c}_{p}}(c)$, соответствующих различным значениям $\mu $.

Поскольку значение R всегда положительно, существует предельное значение скорости V потока, при превышении которого нельзя реализовать максимум функции ${{c}_{p}}(c)$ (аналогичное свойство было отмечено для другого типа ветроустановки в работе [15]).

4. Обсуждение результатов. Получены бифуркационные диаграммы для характерных значений средней безразмерной угловой скорости и среднего коэффициента механической мощности на рабочих режимах устройства, которые соответствуют 2π-периодическим траекториям системы (1.1). Продемонстрировано, что в некоторых диапазонах значений параметров существуют два притягивающих периодических режима, области притяжения которых разделяет отталкивающий периодический режим. При таких значениях параметров имеет место гистерезис безразмерной угловой скорости и коэффициента мощности в зависимости от направления изменения коэффициента c, отвечающего за внешние условия функционирования устройства: скорость потока и нагрузку со стороны потребителей в цепи генератора.

Кроме того, из построенных диаграмм следует наличие в модели притягивающих периодических режимов, характеризующихся относительно высокими значениями безразмерной угловой скорости: $\left\langle z \right\rangle > 1$. На таких режимах скорость точки В крыла на определенных фазах движения превышает по величине скорость ветра. Режимы данного типа характерны для ветроустановок с высокими значениями выходной мощности: турбины пропеллерного типа и типа Дарье [22–25].

Продемонстрированы некоторые тенденции зависимости характеристик рабочих режимов от коэффициента сухого трения в ползуне, которое в модели описывается параметром  f. В безразмерных параметрах и переменных построены серии бифуркационных диаграмм (рис. 2, 3), соответствующих различным значениям безразмерного коэффициента трения.

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

При моделировании внешних воздействий использованы подходы [7–9]. Математическая модель представлена в форме динамической системы второго порядка. Рабочему режиму установки соответствует периодическая орбитально устойчивая траектория системы. Для отыскания периодических решений и проверки их устойчивости применен метод [14]. Проведен параметрический анализ рабочих режимов. В частности, исследовано влияние трения в ползуне механизма при функционировании в различных средах (воздух, вода). Получены оценки коэффициента мощности на рабочих режимах. Продемонстрировано наличие эффекта гистерезиса мощности в зависимости от направления изменения коэффициента полезной нагрузки.

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