Известия РАН. Теория и системы управления, 2020, № 6, стр. 3-14

МЕТОД ФОРМИРОВАНИЯ АВТОРОТАЦИЙ В УПРАВЛЯЕМОЙ МЕХАНИЧЕСКОЙ СИСТЕМЕ С ДВУМЯ СТЕПЕНЯМИ СВОБОДЫ

Л. А. Климина *

НИИ механики МГУ
Москва, Россия

* E-mail: klimina@imec.msu.ru

Поступила в редакцию 30.04.2020
После доработки 19.06.2020
Принята к публикации 27.07.2020

Полный текст (PDF)

Аннотация

Рассматривается автономная механическая система с двумя вращательными степенями свободы. Она находится под действием консервативных и неконсервативных сил. В соответствующей динамической системе присутствуют варьируемые параметры, которые можно рассматривать как коэффициенты усиления управляющих воздействий. Требуется подобрать значения этих параметров таким образом, чтобы сформировать в системе режим авторотации, обладающий определенными свойствами. Предложен метод итерационного поиска соответствующих значений параметров. Подход представляет собой модификацию метода Андронова–Понтрягина и, в отличие от последнего, не подразумевает наличия малого параметра в системе. Приведен пример применения метода к модели ветроэнергетической установки: подобрано значение параметра, обеспечивающее существование режима, который характеризуется высоким значением механической мощности, отбираемой у потока. В динамической системе этому режиму соответствует аттрактор, представляющий собой траекторию, близкую к периодической.

Введение. Рассмотрим механическую систему, текущее положение элементов которой можно описать двумя угловыми координатами. Пусть эта система находится под действием сил, среди которых есть как диссипативные, так и обеспечивающие приток энергии в систему. Предположим, что силы, действующие на систему, не зависят явно от времени. Здесь и далее под термином “авторотация” будем понимать такое притягивающее движение рассматриваемой механической системы, при котором по каждой из упомянутых угловых координат происходит вращение, не переходящее в колебания, не сопровождающееся сменой направления вращения, причем ни одна из соответствующих угловых скоростей не растет неограниченно. Это понимание согласуется с определением, введенным А.А. Андроновым [1, 2].

Необходимость описания авторотаций систем с двумя степенями свободы возникает во многих задачах механики: например, движение тела с неподвижной точкой в среде с сопротивлением [3], движение объекта по поверхности сферы при наличии трения [4], динамика аэроупругих систем [5], динамика двухроторных ветроустановок [6]. Проблема изучения авторотаций таких объектов достаточно сложна, поскольку в общем случае речь идет об отыскании периодических, а также в некотором смысле близких к периодическим (например, квазипериодических) траекторий автономной неконсервативной динамической системы четвертого порядка, описывающей изменение двух угловых координат и двух соответствующих им угловых скоростей.

Если авторотации отвечает периодическое решение, то можно использовать методы поиска периодических траекторий (их краткий обзор содержится, например, в [7]). Для задачи поиска периодического режима авторотации подходят, в частности, алгоритмы, предназначенные для неавтономных систем, поскольку в окрестности такого режима один из углов можно принять за новое время.

Если требуется описать не только периодические, но и близкие к периодическим решения, одним из эффективных методов является усреднение по двум углам. Его целесообразно применять при наличии малого параметра, назовем его a, при правых частях уравнений изменения угловых скоростей (когда угловые скорости можно считать медленными переменными, а углы – быстрыми). Этот метод использован, например, в [3]. Близкий подход применен в [4], где в качестве медленных переменных выступают первые интегралы порождающей консервативной системы, не совпадающие в этом случае с угловыми скоростями. При реализации методики усреднения может потребоваться построение различных усредненных систем “вдали” и “вблизи” от резонансных соотношений медленных переменных [8]. Подобная методика использована для описания автоколебаний связанных осцилляторов в [9].

Известны подходы, в которых начальное приближение для режима авторотации, полученное путем усреднения по двум углам, продолжается численно по параметру a, при этом можно (и целесообразно) варьировать и другие параметры системы [10].

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

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

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

1. Постановка задачи. Рассматривается автономная механическая система с двумя вращательными степенями свободы, находящаяся под действием неконсервативных сил. Пусть в безразмерных переменных динамика системы описывается уравнениями следующего вида (точкой обозначена производная по безразмерному времени $\tau $):

(1.1)
$\left\{ \begin{gathered} {{{\dot {\varphi }}}_{i}} = \frac{{\partial {{H}_{i}}({{\varphi }_{i}},{{\omega }_{i}})}}{{\partial {{\omega }_{i}}}}, \hfill \\ {{{\dot {\omega }}}_{i}} = - \frac{{\partial {{H}_{i}}({{\varphi }_{i}},{{\omega }_{i}})}}{{\partial {{\varphi }_{i}}}} + a{{Q}_{i}}({{\varphi }_{i}},{{\omega }_{1}},{{\omega }_{2}}), \hfill \\ \end{gathered} \right.\quad i = 1,\,2.$
$\begin{gathered} H({{\varphi }_{1}},{{\varphi }_{2}},{{\omega }_{1}},{{\omega }_{2}}) = \sum\limits_{i = 1}^2 {{{H}_{i}}({{\varphi }_{i}},{{\omega }_{i}})} ; \\ {{H}_{i}}({{\varphi }_{i}},{{\omega }_{i}}) = 0.5k_{i}^{{ - 1}}({{\varphi }_{i}})\omega _{i}^{2} + {{P}_{i}}({{\varphi }_{i}}). \\ \end{gathered} $

Здесь φi – обобщенные координаты (углы), ωi – соответствующие безразмерные импульсы; $H({{\varphi }_{1}},{{\varphi }_{2}},{{\omega }_{1}},{{\omega }_{2}})$ – безразмерная полная механическая энергия системы; ${{k}_{i}}({{\varphi }_{i}}) \geqslant 1$ – аналитические функции, определяемые структурой кинетической энергии системы; ${{P}_{i}}({{\varphi }_{i}}) \geqslant 0$ – аналитические функции, описывающие потенциальную энергию системы; a > 0 – параметр, характеризующий неконсервативные воздействия (не предполагается малым). Функции ${{Q}_{i}}$ – безразмерные обобщенные силы (аналитические функции). Все функции в правой части (1.1) 2π-периодические по переменным φi. Отметим, что правая часть уравнения, описывающего изменение переменной ω12), не зависит от угла φ21).

Рассмотрим задачи описания/формирования авторотаций механической системы, поведение которой моделируется уравнениями (1.1). Таким режимам соответствуют траектории системы, на которых переменные ${{\omega }_{i}}$ не меняют знака и являются периодическими или же в некотором смысле близкими к периодическим функциями времени. Без ограничения общности рассмотрим траектории, расположенные в области ${{G}_{0}} = \{ {{\omega }_{1}} > {{\omega }_{0}} > 0,\;{{\omega }_{2}} > {{\omega }_{0}} > 0\} $, где ω0 – некоторое положительное значение. Случаи, когда обе или одна из величин ωi меньше нуля, можно рассмотреть аналогично.

Для описания периодических траекторий можно принять один из углов ${{\varphi }_{i}}$ за новое время и применить методы исследования неавтономных систем (например, [19]).

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

Пусть в систему добавлены 2π-периодические по переменным φi управляющие воздействия вида ${{b}_{i}}{{f}_{i}}({{\varphi }_{i}},{{\omega }_{1}},{{\omega }_{2}})$, где ${{\omega }_{i}}{{f}_{i}}({{\varphi }_{i}},{{\omega }_{1}},{{\omega }_{2}}) > 0$ при ${{\omega }_{i}} \ne 0$, т.е. система (1.1) при наличии управления имеет вид

(1.2)
$\left\{ \begin{gathered} {{{\dot {\varphi }}}_{i}} = \frac{{\partial {{H}_{i}}({{\varphi }_{i}},{{\omega }_{i}})}}{{\partial {{\omega }_{i}}}}, \hfill \\ {{{\dot {\omega }}}_{i}} = - \frac{{\partial {{H}_{i}}({{\varphi }_{i}},{{\omega }_{i}})}}{{\partial {{\varphi }_{i}}}} + a\left( {{{Q}_{i}}({{\varphi }_{i}},{{\omega }_{1}},{{\omega }_{2}}) - {{b}_{i}}{{f}_{i}}({{\varphi }_{i}},{{\omega }_{1}},{{\omega }_{2}})} \right), \hfill \\ \end{gathered} \right.\quad i = 1,\,2.$

Считаем, что функции ${{f}_{i}}({{\varphi }_{i}},{{\omega }_{1}},{{\omega }_{2}})$ заданы (например, ${{f}_{i}} \equiv {\text{sign}}({{\omega }_{i}})$ или ${{f}_{i}} \equiv {{\omega }_{i}}$) и выбору подлежат значения коэффициентов bi усиления управляющих воздействий.

Пусть за счет выбора коэффициентов bi требуется сформировать в системе (1.2) такую траекторию $\left( {{{\varphi }_{i}},{{\omega }_{i}}} \right)$, на которой среднее по времени при $\tau \to \infty $ значение функции ${{\omega }_{i}}(\tau )$ равно $\Omega _{i}^{*}$, т.е.

(1.3)
$\mathop {\lim }\limits_{T \to \infty } \,\frac{1}{T}\int\limits_0^T {{{\omega }_{i}}(\tau )d\tau } = \Omega _{i}^{*}.$

Обозначим такую траекторию через $\left( {{{\varphi }_{i}},{{\omega }_{i}}} \right) = \left( {{{{\hat {\varphi }}}_{i}}(\tau ),{{{\hat {\omega }}}_{i}}(\tau )} \right)$.

2. Частично усредненные системы. Опишем первый этап процедуры определения коэффициента b1 управления. Для b2 соответствующий этап полностью аналогичен. Рассмотрим вспомогательную частично усредненную систему второго порядка вида

(2.1)
$\left\{ \begin{gathered} {{{\dot {\varphi }}}_{1}} = \frac{{\partial {{H}_{1}}({{\varphi }_{1}},{{\omega }_{1}})}}{{\partial {{\omega }_{1}}}}, \hfill \\ {{{\dot {\omega }}}_{1}} = - \frac{{\partial {{H}_{1}}({{\varphi }_{1}},{{\omega }_{1}})}}{{\partial {{\varphi }_{1}}}} + a({{Q}_{1}}({{\varphi }_{1}},{{\omega }_{1}},\Omega _{2}^{*}) - {{b}_{1}}{{f}_{1}}({{\varphi }_{1}},{{\omega }_{1}},\Omega _{2}^{*})). \hfill \\ \end{gathered} \right.$

В уравнениях (2.1) величина $\Omega _{2}^{*}$ выступает как параметр системы. Условно говоря, при построении системы (2.1) на основе (1.2) переменную ${{\omega }_{2}}$ как аргумент функций ${{Q}_{1}}$ и  f1 усреднили по времени.

Замечание. При построении системы (2.1) существенно, что функции ${{H}_{1}}$, ${{Q}_{1}}$ и  f1 не зависят явно от угла ${{\varphi }_{2}}$. Дело в том, что информации о среднем по времени (в смысле (1.3)) значении $\Omega _{2}^{*}$ переменной ${{\omega }_{2}}$ не достаточно, вообще говоря, чтобы определить среднее по времени значение функции, зависящей от ${{\varphi }_{2}}$.

Будем искать значение параметра b1, обеспечивающее в системе (2.1) существование 2π-периодической по переменной ${{\varphi }_{1}}$ траектории, проходящей через точку ${{\varphi }_{1}} = 0,\,\;{{\omega }_{1}} = {{\omega }_{{10}}}$. Обозначим такую траекторию ${{\omega }_{1}} = {{\tilde {\omega }}_{1}}({{\varphi }_{1}};{{\omega }_{{10}}})$. Решать эту задачу можно с помощью различных подходов, например [19, 20]. В настоящей работе используем итерационный алгоритм, который предложен в [7], основан на модификации метода Андронова–Понтрягина [1, 2, 21] и существенно использует то, что параметр ${{b}_{1}}$ в системе (2.1) отвечает за поворот векторного поля (в смысле определения, данного в [2]). Отметим, что система (2.1) относится именно к тому классу, для которого указанный подход применим. В соответствии с этим алгоритмом определим нулевое и последующие приближения ${{y}_{n}}({{\varphi }_{1}},{{\omega }_{{10}}})$ для траектории ${{\tilde {\omega }}_{1}}({{\varphi }_{1}};{{\omega }_{{10}}})$:

${{y}_{0}}(\varphi ,{{\omega }_{{10}}}) = \sqrt {2{{k}_{1}}(\varphi )({{h}_{{10}}} - {{P}_{1}}(\varphi ))} ,\quad {{h}_{{10}}} = {{H}_{1}}(0,{{\omega }_{{10}}}),$
(2.2)
$\begin{gathered} {{b}_{{10}}}({{\omega }_{{10}}}) = {{\left( {\int\limits_0^{2\pi } {{{f}_{1}}(\varphi ,{{y}_{0}}(\varphi ,{{\omega }_{{10}}}),\Omega _{2}^{*})\,} d\varphi } \right)}^{{ - 1}}}\int\limits_0^{2\pi } {{{Q}_{1}}(\varphi ,{{y}_{0}}(\varphi ,{{\omega }_{{10}}}),\Omega _{2}^{*})\,} d\varphi ; \\ {{y}_{n}}(\varphi ,{{\omega }_{{10}}}) = \operatorname{Re} \sqrt {2{{k}_{1}}(\varphi )\left( {{{h}_{{10}}} - {{P}_{1}}(\varphi ) + a\int\limits_0^\varphi {({{Q}_{1}}(\vartheta ,{{y}_{{n - 1}}}(\vartheta ,{{\omega }_{{10}}}),\Omega _{2}^{*}) - {{b}_{{1(n - 1)}}}{{f}_{1}}(\vartheta ,{{y}_{{n - 1}}}(\vartheta ,{{\omega }_{{10}}}),\Omega _{2}^{*}))} \,d\vartheta } \right)} , \\ \end{gathered} $
${{b}_{{1n}}}({{\omega }_{{10}}}) = {{\left( {\int\limits_0^{2\pi } {{{f}_{1}}(\varphi ,{{y}_{n}}(\varphi ,{{\omega }_{0}}),\Omega _{2}^{*})\,} d\varphi } \right)}^{{ - 1}}}\int\limits_0^{2\pi } {{{Q}_{1}}(\varphi ,{{y}_{n}}(\varphi ,{{\omega }_{0}}),\Omega _{2}^{*})\,} d\varphi ,\quad n \geqslant 1.$

Значение ω10 будем варьировать в некотором интервале ${{\omega }_{{10}}} \in [{{A}_{1}},{{B}_{1}}]$, таком, чтобы внутри этого интервала можно было ожидать расположение той величины ω10, которая соответствует периодической траектории со средним по времени значением переменной ω1, равным $\Omega _{1}^{*}$. Далее будем предполагать, что итерационная процедура (2.2) сходится при каждом ${{\omega }_{{10}}} \in [{{A}_{1}},{{B}_{1}}]$. Тогда, как показано в [7], функция ${{\tilde {\omega }}_{1}}({{\varphi }_{1}};{{\omega }_{{10}}}) = \mathop {\lim }\limits_{n \to \infty } {{y}_{n}}({{\varphi }_{1}},{{\omega }_{{10}}})$ описывает 2π-периодическую по переменной φ1 траекторию, которая существует в системе (2.1) при ${{b}_{1}} = {{\tilde {b}}_{1}}({{\omega }_{{10}}}) = \mathop {\lim }\limits_{n \to \infty } {{b}_{{1n}}}({{\omega }_{{10}}})$ и проходит через точку ${{\varphi }_{1}} = 0,$ ${{\omega }_{1}} = {{\omega }_{{10}}}$.

Вычислим среднее по времени за период значение ${{\bar {\Omega }}_{1}}({{\omega }_{{10}}})$ переменной ω1 вдоль траектории ${{\tilde {\omega }}_{1}}({{\varphi }_{1}};{{\omega }_{{10}}})$ (в общем случае оно отличается от среднего по углу):

${{\bar {\Omega }}_{1}}({{\omega }_{{10}}}) = 2\pi {{\left( {\int\limits_0^{2\pi } {{{{\left( {{{{\tilde {\omega }}}_{1}}(\varphi ;{{\omega }_{{10}}})} \right)}}^{{ - 1}}}d\varphi } } \right)}^{{ - 1}}}.$

Определим то значение ${{\omega }_{{10}}} = \omega _{{10}}^{*}$, при котором ${{\bar {\Omega }}_{1}}({{\omega }_{{10}}}) = \Omega _{1}^{*}$. Оно единственно (если существует) в силу свойств поворота поля. Факт непересечения периодических траекторий, которые существуют при различных значениях параметра, отвечающего за поворот поля, справедлив для систем более широкого класса, чем (2.1) [22]. Обозначим $b_{1}^{*} = {{\tilde {b}}_{1}}(\omega _{{10}}^{*})$.

Аналогично на основе построения частично усредненной системы для переменных φ2, ω2 определим значение $b_{2}^{*}$, обеспечивающее в такой системе существование $2\pi $‑периодической по переменной φ2 траектории ${{\tilde {\omega }}_{2}}({{\varphi }_{2}};\omega _{{20}}^{*})$, вдоль которой среднее по времени значение ${{\bar {\Omega }}_{2}}({{\omega }_{{20}}})$ переменной ω2 равно $\Omega _{2}^{*}$ (предполагаем, что искомое значение существует).

Предположение 1. При ${{b}_{i}} = b_{i}^{*}$ в системе (1.2) можно ожидать существования близкой к периодической траектории $\left( {{{{\hat {\varphi }}}_{i}}(\tau ),{{{\hat {\omega }}}_{i}}(\tau )} \right)$, на которой среднее по времени при $\tau \to \infty $ значение функции ${{\hat {\omega }}_{i}}(\tau )$ равно $\Omega _{i}^{*}$, т.е. выполнено (1.3).

Если предположение 1 для конкретных $\Omega _{i}^{*}$ верно, то проекции траектории $\left( {{{{\hat {\varphi }}}_{i}}(\tau ),{{{\hat {\omega }}}_{i}}(\tau )} \right)$ на плоскости $({{\varphi }_{i}},{{\omega }_{i}})$ приближенно описываются построенными функциями ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$. До проверки предположения 1 будем называть кривую $({{\varphi }_{i}},{{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*}))$ “претендентом” на приближение для близкой к периодической траектории.

Будем устанавливать справедливость предположения 1 численно путем прямого интегрирования системы (1.2) с начальными условиями, расположенными на кривых ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$. Следует отметить, что если при этом траектория $\left( {{{{\hat {\varphi }}}_{i}}(\tau ),{{{\hat {\omega }}}_{i}}(\tau )} \right)$ существует, но является неустойчивой как в прямом, так в обратном времени, то для проверки предположения 1 не достаточно прямого интегрирования. Тогда для верификации можно использовать методы численного исследования, примененные в [10].

2. Модификации подхода. Рассмотрим важные смежные задачи, которые в общем случае более сложны, чем исходная задача выбора двух коэффициентов ${{b}_{i}}$.

Модификация 1. Пусть требуется обеспечить в системе (1.2) существование близкой к периодической траектории, характеризующейся заданным значением одной скалярной величины, и пусть при этом выбору подлежит только один параметр в законе управления. Иными словами, в системе имеет место дефицит управляющих воздействий (смысл этого термина обсуждается в работе [23]). Например, пусть требуется, чтобы среднее по времени значение функции ${{\hat {\omega }}_{1}}(\tau )$ составляло $\Omega _{1}^{*}$, и пусть при этом ${{b}_{1}} = 0$. Тогда используем предложенный выше алгоритм, варьируя $\Omega _{2}^{*}$ как параметр в некотором диапазоне. Построим зависимость $b_{1}^{*}(\Omega _{2}^{*})$. Будем искать значение $\Omega _{2}^{*}$, обеспечивающее $b_{1}^{*}(\Omega _{2}^{*}) = 0$ (возможно, таких значений несколько). Для такого $\Omega _{2}^{*}$ (если оно будет определено) вычислим значение $b_{2}^{*}$, которое и представляет собой искомый параметр закона управления.

Модификация 2. Принципиально важна следующая модификация. Предположим, требуется найти расположенные в некоторой области близкие к периодическим (например, квазипериодические) траектории системы, в которой отсутствует управление, т.е. системы (1.1). Формально расширим систему до вида (1.2). Будем варьировать значения $\Omega _{i}^{*}$ в некоторых интервалах. Построим зависимости $b_{i}^{*}(\Omega _{1}^{*},\Omega _{2}^{*})$, используя описанную выше процедуру для частично усредненных систем. Найдем те пары $(\Omega _{1}^{*},\Omega _{2}^{*})$, при которых $b_{i}^{*}(\Omega _{1}^{*},\Omega _{2}^{*}) = 0$. Каждой такой паре может соответствовать близкая к периодической траектория системы (1.1) – в смысле предположения 1. Требуется дополнительная численная проверка, аналогичная описанному выше основному случаю. Отметим, что в этом случае нулевое приближение метода совпадает с отысканием неподвижных точек системы, которая получается из уравнений (1.1) при переходе к переменным ${{\varphi }_{i}}$, Hi и усреднении по двум углам. Это легко показать, рассмотрев формулы для нулевой итерации алгоритма [7].

3. Пример применения в задаче о двойной ветротурбине Дарье. Рассмотрим механическую систему, состоящую из двух вертикально-осевых ветротурбин Дарье, расположенных в стационарном потоке ветра скорости V (рис. 1) плотности ρ. Валы турбин соосны. Лопасти турбин ориентированы так, чтобы поток поддерживал вращение турбин в противоположных направлениях. На валу первой турбины закреплен статор электрогенератора, на валу второй – ротор. Генератор подключен к локальной цепи с внешним сопротивлением R.

Рис. 1.

Схема вертикально-осевой ветроустановки, состоящей из двух турбин Дарье

Обозначим через ${{\varphi }_{1}}$ и ${{\varphi }_{2}}$ углы поворота нижней и верхней турбин. Эти углы отсчитываются от направления скорости потока против часовой стрелки и по часовой стрелке соответственно, если наблюдать с конца оси Oz вращения турбин.

Предполагается, что ${{\varphi }_{1}}$ – это угол между направлением ветра и державкой одной из лопастей нижней турбины, которая выбрана в качестве “первой” лопасти. Тогда угол между направлением ветра и державкой лопасти “с номером k” составляет ${{\varphi }_{1}} + 2\pi (k - 1){\text{/}}3$. Заменяя ${{\varphi }_{1}}$ на ${{\varphi }_{2}}$, получаем аналогичные формулы для углов положения радиальных державок верхней турбины.

Пусть J1, J2 – моменты инерции относительно оси Oz для первой и второй турбин; S – характерная площадь одной лопасти. Лопасти находятся под действием потока воздуха. Аэродинамический профиль каждой лопасти характеризуется коэффициентами ${{C}_{d}}(\alpha )$ и ${{C}_{l}}(\alpha )$ силы сопротивления и подъемной (боковой) силы соответственно; $\alpha $ – мгновенный угол атаки. При численных расчетах будем использовать для аэродинамических коэффициентов аппроксимации (рис. 2) экспериментальных данных, которые соответствуют профилю RAF34 (результаты экспериментов взяты из [24]). Считаем, что хорда профиля лопасти перпендикулярна направлению от центра давления (точки приложения аэродинамических сил) к оси вращения ротора. Пренебрегаем возможным перемещением центра давления вдоль профиля. Расстояние от центра давления до оси Oz вращения ротора обозначим через r (одинаково для всех лопастей).

Рис. 2.

Аэродинамические характеристики профиля лопасти

Взаимодействие между ротором и статором генератора описывается электромеханическим моментом. Предполагается, что этот момент является линейной функцией угловой скорости $({{\omega }_{1}} + {{\omega }_{2}})$ ротора относительно статора [25]:

$T = \frac{C}{{(\sigma + R)}}({{\omega }_{1}} + {{\omega }_{2}}).$

Здесь С – коэффициент электромеханического взаимодействия, σ – внутреннее сопротивление генератора. Момент T действует со стороны статора на ротор и направлен так, что замедляет вращение второй турбины (соединенной с ротором генератора). Такой же по величине момент действует со стороны ротора на статор и направлен так, что замедляет вращение первой турбины (соединенной со статором).

Обозначим через ${{\omega }_{1}} = r{{\dot {\varphi }}_{1}}{\text{/}}V$ быстроходность (безразмерная угловая скорость) первого ротора, ${{\omega }_{2}} = r{{\dot {\varphi }}_{2}}{\text{/}}V$ – второго (точкой обозначена производная по времени t).

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

$\begin{gathered} \left\{ {\begin{array}{*{20}{l}} {{{{\dot {\varphi }}}_{1}} = {{\omega }_{1}},} \\ {{{{\dot {\varphi }}}_{2}} = {{\omega }_{2}},} \\ {{{{\dot {\omega }}}_{1}} = a\left( {W({{\varphi }_{1}},{{\omega }_{1}}) - b({{\omega }_{1}} + {{\omega }_{2}})} \right),} \\ {{{{\dot {\omega }}}_{2}} = aB\left( {W({{\varphi }_{2}},{{\omega }_{2}}) - b({{\omega }_{1}} + {{\omega }_{2}})} \right),} \end{array}} \right. \\ a = \frac{{\rho S{{r}^{3}}}}{{2{{J}_{1}}}},\quad B = \frac{{{{J}_{1}}}}{{{{J}_{2}}}},\quad b = \frac{{2C}}{{V\rho S{{r}^{2}}(\sigma + R)}}, \\ \end{gathered} $
(3.1)
$W(\varphi ,\omega ) = \sum\limits_{k = 1}^3 {\sqrt {u_{k}^{2} + w_{k}^{2}} } \left( {{{C}_{l}}({{\alpha }_{k}}){{u}_{k}} - {{C}_{d}}({{\alpha }_{k}}){{w}_{k}}} \right).$
$\begin{gathered} {{\alpha }_{k}} = \left\{ \begin{gathered} \arctan \left( {\frac{{{{u}_{k}}(\varphi )}}{{{{w}_{k}}(\varphi ,\omega )}}} \right),\quad {\text{если}}\quad {{w}_{k}}(\varphi ,\omega ) \geqslant 0, \hfill \\ \arctan \left( {\frac{{{{u}_{k}}(\varphi )}}{{{{w}_{k}}(\varphi ,\omega )}}} \right) + \pi ,\quad {\text{если}}\quad {{w}_{k}}(\varphi ,\omega ) < 0, \hfill \\ \end{gathered} \right. \\ {{u}_{k}} = \cos \left( {\varphi + \frac{{2\pi }}{3}k} \right),\quad {{w}_{k}} = \omega + \sin \left( {\varphi + \frac{{2\pi }}{3}k} \right). \\ \end{gathered} $

Здесь точкой обозначена производная по безразмерному времени $\tau = Vt{\text{/}}r$; величины ${{u}_{k}},{{w}_{k}}$ являются промежуточными переменными (k – номер лопасти). При вычислении значений ${{u}_{k}},{{w}_{k}}$ для первой/второй турбины подставляем $\varphi = {{\varphi }_{i}},\;\omega = {{\omega }_{i}}$ соответственно. Коэффициент B (отношение моментов инерции турбин) при дальнейших расчетах считаем равным 0.5. Безразмерный параметр b – коэффициент внешней нагрузки – может варьироваться за счет изменения величины R внешнего сопротивления в цепи и выступать как коэффициент усиления управляющего воздействия (аналогичный подход к задачам управления режимами функционирования ветроустановок рассмотрен в [2628]).

Для случая $a \ll 1$ зависимость ωi на авторотациях от параметра b описана в [6] на основе анализа системы, усредненной по двум углам. В частности, установлено, что в случае $a \ll 1$ максимальное значение механической мощности установки достигается на притягивающем режиме авторотации, на котором каждая безразмерная угловая скорость близка к значению ${{\omega }_{{opt}}} \approx 10$, соответствующее значение коэффициента внешней нагрузки $b = {{b}_{{opt}}} \approx 0.11$.

3.2. Аттракторы в случае a = 1. Рассмотрим случай a = 1. Пусть требуется определить значение b* коэффициента b, такое, что существует притягивающий режим авторотации, на котором среднее по времени значение ${{\omega }_{1}}$ составляет $\Omega _{1}^{*} = {{\omega }_{{opt}}}$, а величина $\Omega _{2}^{*}$ (среднее от ${{\omega }_{2}}$) принимает значение, максимально близкое к ${{\Omega }_{{opt}}}$ (модификация 1). Пусть помимо этого требуется для случая $b = b{\kern 1pt} *$ описать все аттракторы, существующие у системы (3.1) в области ${{G}_{0}}$ (модификация 2).

При $a = 1$ построим вспомогательные системы вида (2.1) для исходных уравнений (3.1) и применим предложенный итерационный метод в модификации 1. Получаем, что оптимальное с позиций максимизации механической мощности значение параметра b остается практически тем же, что и при $a \ll 1$: $b{\kern 1pt} * \approx 0.11$. Это связано с тем, что при относительно больших значениях величин ${{\omega }_{i}}$ зависимость правой части (3.1) от углов ${{\varphi }_{i}}$ несущественна.

Для описания аттракторов, существующих в случае $b = b{\kern 1pt} *$, построим последовательные приближения для близких к периодическим траекторий системы (3.1) при b = 0.11 на основе модификации 2. Получаем, что при таком значении коэффициента b в области положительных $\Omega _{i}^{*}$ существует пять пар $(\Omega _{1}^{*},\Omega _{2}^{*})$, обеспечивающих $b_{1}^{*} = b_{2}^{*} = 0.11$. Прямое численное интегрирование системы (3.1) с начальными условиями ${{\varphi }_{i}}(0) = 0$, ${{\omega }_{i}}(0) = \omega _{{i0}}^{*}$ показывает, что две из найденных пар $(\Omega _{1}^{*},\Omega _{2}^{*})$ отвечают аттракторам системы (3.1). Об одной из них уже шла речь выше: для нее $\Omega _{i}^{*} \approx {{\omega }_{{opt}}}$, а значения ${{\hat {\omega }}_{i}}(\tau )$ на соответствующем притягивающем решении системы (3.1) с небольшими отклонениями колеблются около константы ${{\omega }_{{opt}}}$.

Другая пара: $(\Omega _{1}^{*},\Omega _{2}^{*}) = (0.43,0.49)$. Она отвечает притягивающему режиму авторотации с низкими угловыми скоростями. На рис. 3 черным цветом изображены проекции на плоскость $(\sin 3{{\varphi }_{i}},{{\omega }_{i}})$ соответствующих построенных итерационно приближенных кривых ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$. При их построении выполнено 20 итераций вида (2.2). Зеленым цветом показано приближение для аттрактора, построенное путем прямого численного интегрирования системы (3.1) методом Рунге–Кутты (с течением времени проекция траектории заполняет область, закрашенную зеленым цветом).

Рис. 3.

Приближения для аттрактора, построенные предложенным итерационным алгоритмом (черный) и прямым численным интегрированием (зеленый). Случай b = 0.11

Из рис. 3 видно, что построенные приближения ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$ позволили хорошо описать интегральное поведение близкой к периодической траектории $\left( {{{{\hat {\varphi }}}_{i}}(\tau ),{{{\hat {\omega }}}_{i}}(\tau )} \right)$.

В отношении эволюции аттракторов системы с ростом параметра a интересно отметить, что чем ближе аттрактор системы к плоскости ${{\omega }_{1}} = 0$, ${{\omega }_{2}} = 0$, тем он чувствительнее к изменению параметра a. Например, при a = 1, $b = 0.11$ для аттрактора, отвечающего $\Omega _{i}^{*} \approx {{\omega }_{{opt}}}$, значения безразмерных угловых скоростей слабо отличаются от констант, а для другого аттрактора, существующего при тех же параметрах, отклонение ${{\hat {\omega }}_{i}}(\tau )$ от констант существенно.

3.3. Случай b = 0.5. Для сравнения рассмотрим эволюцию с ростом a аттрактора, который существует в случае b = 0.5 при относительно небольших значениях ${{\omega }_{i}}$. При $a \ll 1$ этому аттрактору соответствует неподвижная точка $P = (\Omega _{1}^{*},\Omega _{2}^{*}) = (0.13,0.13)$ усредненной системы нулевого приближения. Применяя метод последовательных итерационных приближений (модификация 2), получаем, что при a = 0.1 аттрактору отвечает пара $(\Omega _{1}^{*},\Omega _{2}^{*}) = (0.11,0.15)$. На рис. 4 черным цветом изображены проекции соответствующих построенных итерационно приближенных кривых ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$ на плоскости $(\sin 3{{\varphi }_{i}},{{\omega }_{i}})$. Зеленым цветом показано приближение для аттрактора, построенное путем прямого численного интегрирования системы (3.1) методом Рунге–Кутты (аналогично рис. 3). Относительные колебания ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$ около средних значений $\Omega _{i}^{*}$ очень существенны, как видно из рисунка.

Рис. 4.

Приближения для аттрактора, построенные предложенным итерационным алгоритмом (черный) и прямым численным интегрированием (зеленый). Случай b = 0.5

Интересно отметить, что для случая b = 0.5 в усредненной системе нулевого приближения есть только одна неподвижная точка, расположенная в области G0 (указанная выше точка P). Иными словами, при достаточно малых a существует только одна пара $(\Omega _{1}^{*},\Omega _{2}^{*})$, $\Omega _{i}^{*} > 0$, удовлетворяющая условию $b_{i}^{*}(\Omega _{1}^{*},\Omega _{2}^{*}) = 0.5$. При увеличении параметра a точка P существенно смещается, и помимо этого, в область ${{G}_{0}}$ через ее границу “приходят” дополнительные пары $(\Omega _{1}^{*},\Omega _{2}^{*})$, обеспечивающие $b_{i}^{*} = 0.5$.

Здесь и далее дополнительными будем называть пары $(\Omega _{1}^{*},\Omega _{2}^{*})$, существующие при конечных значениях a, которые при продолжении по параметру a в направлении уменьшения этого параметра (вплоть до $a \ll 1$) исчезают (сливаясь друг с другом) или покидают область ${{G}_{0}}$.

В частности, при a = 0.1 существует дополнительная пара $(\Omega _{1}^{*},\Omega _{2}^{*}) = (0.10,0.16)$. Прямое численное интегрирование системы с начальными условиями в окрестности кривых ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$, соответствующих данной дополнительной паре, показывает, что вблизи этих кривых проходит граница областей притяжения различных аттракторов (как в прямом времени, так и в обратном времени). В прямом времени при начальных условиях $(0,0,\omega _{{10}}^{*},\omega _{{20}}^{*})$ система выходит на притягивающий режим, описанный выше на рис. 4. Соответствующий переходный процесс показан зелеными кривыми на рис. 5. При начальных условиях $(0,0,\omega _{{10}}^{*} - 0.005,\omega _{{20}}^{*} + 0.005)$ первая турбина переходит на режим колебаний (с очень небольшой амплитудой) около положения, в котором ${{\varphi }_{1}} = 0.11$, а вторая турбина при этом выходит на режим ротации. Траектория с указанными начальными условиями изображена красными кривыми на рис. 5. В обратном времени при начальных условиях $(0,0,\omega _{{10}}^{*},\omega _{{20}}^{*})$ обе угловые скорости уходят в плюс бесконечность, а из точки $(0,0,\omega _{{10}}^{*} - 0.005,\omega _{{20}}^{*} + 0.005)$ – в минус бесконечность.

Рис. 5.

“Претендент” на приближение для неустойчивой близкой к периодической траектории, построенный предложенным итерационным алгоритмом (черный); траектории системы (3.1), полученные прямым численным интегрированием с начальными условиями $(0,0,\omega _{{10}}^{*},\omega _{{20}}^{*})$ (зеленый) и $(0,0,\omega _{{10}}^{*} - 0.005,\omega _{{20}}^{*} + 0.005)$ (красный). Случай b = 0.5

Таким образом, можно предположить, что кривые ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$ приближенно описывают близкую к периодической траекторию, которая неустойчива и в прямом, и в обратном времени (траектория седлового типа). Для уточнения этого предположения можно использовать численный алгоритм, примененный в [10]. Назовем полученную пару ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$ “претендентом” на приближение для неустойчивой близкой к периодической траектории.

При построении кривых ${{\tilde {\omega }}_{i}}({{\varphi }_{i}};\omega _{{i0}}^{*})$ в случае b = 0.5, a = 0.1 было выполнено по 30 итерационных приближений вида (2.2) для каждого рассмотренного значения $\omega _{{i0}}^{*}$.

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

Соответствующая динамическая система (1.2) имеет четвертый порядок. Рассматривается случай, когда авторотации отвечает близкая к периодической, но не периодическая траектория. Предложен метод выбора коэффициентов усиления управляющих воздействий для формирования близкой к периодической траектории с заданными средними по времени значениями угловых скоростей, отвечающих двум угловым координатам. Метод основан на построении двух вспомогательных систем второго порядка и применении к каждой из них итерационного алгоритма [7]. Предложенный подход определяет набор коэффициентов управления, а также “претендента” на приближение для близкой к периодической траектории в фазовом пространстве. Требуется дополнительная прямая численная проверка существования близкой к периодической траектории в окрестности найденного “претендента”.

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

Рассмотрено применение метода в задаче о динамике двойной ветротурбины Дарье. Предложенный метод позволил оценить влияние параметров модели на характеристики рабочих режимов устройства. В частности, описаны некоторые элементы эволюции аттракторов при изменении параметров, а также построен “претендент” на приближение для близкой к периодической траектории неустойчивой и в прямом, и в обратном времени.

Отметим, что на основе методов [29, 30] возможно создание модификации подхода, предложенного в настоящей работе, которая позволит искать приближения для близких к периодическим траекторий автоколебательного типа.

Список литературы

  1. Andronov A.A. Les Cycles Limites de Poincaré et la Théorie des Oscillations Auto-Entretenues // Comptes Rendus de l’Academie des Sci. 1929. V. 189. P. 559–561.

  2. Андронов А.А., Леонтович Е.А., Гордон И.И., Майер А.Г. Теория бифуркаций динамических систем на плоскости. М.: Наука, 1967.

  3. Акуленко Л.Д., Лещенко Д.Д., Черноусько Ф.Л. Быстрое движение вокруг неподвижной точки тяжелого твердого тела в сопротивляющейся среде // Изв. РАН. МТТ. 1982. № 3. С. 5–13.

  4. Шалимова Е.С. Стационарные и периодические режимы в задаче о движении тяжелой точки по вращающейся сфере при наличии вязкого трения // Вестн. МГУ. Сер. 1. Математика, механика. 2014. № 4. С. 43–50.

  5. Selyutskiy Yu.D. Dynamics of an Elastically Mounted Wing in Flow // AIP Conference Proceedings. 2018. V. 2046. P. 020085-1–020085-6.

  6. Досаев М.З., Климина Л.А., Локшин Б.Я., Селюцкий Ю.Д., Шалимова Е.С. О режимах авторотации двухроторной ветротурбины Дарье // Изв. РАН. МТТ. 2021. № 1.

  7. Климина Л.А. Метод построения периодических решений в управляемой динамической системе с цилиндрическим фазовым пространством // Изв. РАН. ТиСУ. 2020. № 2. С. 5–16.

  8. Волосов В.М. О методе усреднения // Докл. АН ССР. 1961. Т. 137. № 1. С. 21–24.

  9. Кондрашов Р.Е., Морозов А.Д. О глобальном поведении решений системы двух уравнений Дюффинга–Ван дер Поля // Нелинейная динамика. 2011. Т. 7. № 3. С. 437–449.

  10. Kamiyama K., Komuro M., Endo T. Bifurcation of Quasi-Periodic Oscillations in Mutually Coupled Hard-Type Oscillators: Demonstration of Unstable Quasi-Periodic Orbits // Intern. J. Bifurcation and Chaos. 2012. V. 22. № 6. P. 1230022.

  11. Schilder F., Osinga H.M., Vogt W. Continuation of Quasi-Periodic Invariant Tori // SIAM J. Applied Dynamical Systems. 2005. V. 4. № 3. P. 459–488.

  12. Parker T.S., Chua L.O. Practical Numerical Algorithms for Chaotic Systems. N.Y.: Springer-Verlag, 1989.

  13. Kaas-Petersen Chr. Computation of Quasiperiodic Solutions of Forced Dissipative Systems // J. Comput. Phys. 1985. V. 58. № 3. P. 395–408.

  14. Боголюбов Н.Н., Митропольский Ю.А., Самойленко А.М. Метод ускоренной сходимости в нелинейной динамике. Киев: Наук. думка, 1969.

  15. Bush J., Gameiro M., Harker S., Kokubu H., Mischaikow K., Obayashi I., Pilarczyk P. Combinatorial-Topological Framework for the Analysis of Global Dynamics // Chaos: An Interdisciplinary J. Nonlinear Science. 2012. V. 22. № 4. P. 047508.

  16. Kamiyama K., Komuro M., Endo T. Algorithms for Obtaining a Saddle Torus Between two Attractors // Intern. J. Bifurcation and Chaos. 2013. V. 23. № 9. P. 1330032.

  17. Zhou B., Thouverez F., Lenoir D. A Variable-Coefficient Harmonic Balance Method for the Prediction of Quasi-Periodic Response in Nonlinear Systems // Mechanical Systems and Signal Processing. 2015. V. 64. P. 233–244.

  18. Chen G., Dunne J.F. A Fast Continuation Scheme for Accurate Tracing of Nonlinear Oscillator Frequency Response Functions // J. Sound and Vibration. 2016. V. 385. P. 284–299.

  19. Самойленко А.М. Численно-аналитический метод исследования периодических систем обыкновенных дифференциальных уравнений. I // Укр. мат. журн. 1965. Т. 17. № 4. С. 82–93.

  20. Гринь А.А., Рудевич С.В. Признак Дюлака–Черкаса для установления точного числа предельных циклов автономных систем на цилиндре // Дифференц. уравнения. 2019. Т. 55. № 3. С. 328–336.

  21. Понтрягин Л.С. О динамических системах, близких к гамильтоновым // ЖЭТФ. 1934. Т. 4. № 9. С. 234–238.

  22. Черкас Л.А., Гринь А.А., Булгаков В.И. Конструктивные методы исследования предельных циклов автономных систем второго порядка (численно-алгебраический подход). Гродно: ГрГУ, 2013.

  23. Formalskii A.M. Stabilization and Motion Control of Unstable Objects. Berlin/Boston: Walter de Gruyter, 2015.

  24. Кравец А.С. Характеристики авиационных профилей. М.: Гос. изд. оборонной промышленности, 1939.

  25. Досаев М.З., Самсонов В.А., Селюцкий Ю.Д., Вен-Лон Лю, Чин-Хуэй Линь. Бифуркации режимов функционирования малых ветроэлектростанций и оптимизация их характеристик // Изв. РАН. МТТ. 2009. № 2. С. 59–66.

  26. Климина Л.А., Шалимова Е.С. Двухпропеллерная ветроэнергетическая установка с дифференциальной планетарной передачей // Мехатроника, автоматизация, управление. 2017. Т. 18. № 10. С. 679–684.

  27. Ишханян М.В., Климина Л.А. Ветротурбина класса “Савониус-Магнус” с коническими лопастями: динамика и управление // Изв. РАН. ТиСУ. 2020. № 4. С. 160–168.

  28. Досаев М.З., Климина Л.А., Селюцкий Ю.Д., Шалимова Е.С. Двойная ветротурбина, работающая на основе эффекта Магнуса // Вестн. МГУ. Сер. 1. Математика, механика. 2020. № 4. С. 65–69.

  29. Климина Л.А. Метод поиска периодических траекторий центрально-симметричных динамических систем на плоскости // Дифференц. уравнения. 2019. Т. 55. № 2. С. 159–168.

  30. Климина Л.А., Селюцкий Ю.Д. Метод построения периодических решений в управляемой динамической системе второго порядка // Изв. РАН. ТиСУ. 2019. № 4. С. 3–15.

Дополнительные материалы отсутствуют.