Разработка алгоритма расчета дисбаланса винтов вертолетов с соосной схемой
В работе представлен алгоритм определения дисбаланса винтов вертолета с соосной схемой вращения по данным вибрационных измерений. Для решения задачи разработана математическая модель динамики системы с шестью степенями свободы, учитывающая поступательные и угловые перемещения конструкции, противоположное направление вращения роторов, гироскопические эффекты и возмущающие воздействия, обусловленные дисбалансом. Для получения обобщённых координат по экспериментальным данным используются траектории трёх неколлинеарно расположенных трёхосевых вибродатчиков. На основе взаимного движения точек восстанавливаются поступательная и вращательная составляющие движения конструкции. Полученные траектории преобразуются в обобщённые координаты динамической модели и используются для расчёта вибрационного отклика системы. Предложенный подход позволяет связать измеряемые вибрационные параметры с характеристиками дисбаланса верхнего и нижнего роторов и использовать математическую модель для разработки алгоритма их определения.
дисбаланс ротора, соосные винты, вибрационный отклик, многотельная динамика, математическое моделирование
1 Введение
Дисбаланс винтов вертолётов является одной из ключевых причин возникновения вибраций в системе приводов и оказывает существенное влияние на надежность, ресурс и эксплуатационные характеристики летательного аппарата [1–3]. Даже небольшое отклонение распределения массы относительно оси вращения приводит к появлению периодических центробежных сил, амплитуда которых возрастает пропорционально квадрату угловой скорости вращения. Возникающие вибрации передаются на элементы конструкции, увеличивая динамические нагрузки и ускоряя износ узлов.
Особую сложность представляет балансировка вертолётов с соосной схемой несущей системы, в которой два ротора расположены на одной оси и вращаются в противоположных направлениях. Вследствие такого взаимодействия вибрационный отклик определяется совместным воздействием дисбалансов обоих роторов, их фазовых соотношений, взаимного расположения и динамических характеристик конструкции [4–9]. В отличие от однороторных систем, измеренный вибрационный сигнал не позволяет непосредственно разделить вклад каждого ротора, что существенно усложняет задачу идентификации параметров дисбаланса.
Для контроля технического состояния вертолётов применяются системы вибродиагностики, использующие акселерометры и методы спектрального анализа для выделения гармонических составляющих, связанных с вращением роторов и различными источниками возмущений [10, 11]. Подобные методы позволяют обнаруживать наличие дисбаланса и оценивать его величину, однако для соосной системы дополнительной задачей является разделение совместного вибрационного отклика на составляющие, обусловленные верхним и нижним роторами.
Существующие методы балансировки, как правило, основаны на последовательном внесении корректирующих воздействий и анализе изменения вибрационного отклика [3, 12–14]. Для несущих винтов вертолётов эта задача решается как совместная балансировка и выравнивание лопастей по траектории [15–19], для двухроторных систем — методом коэффициентов влияния [20, 21]. При этом точность определения корректирующих параметров непосредственно зависит от количества и качества проведенных измерений. Для соосных роторов число независимых параметров возрастает за счёт необходимости определения характеристик дисбаланса каждого ротора, включая его величину и фазовое положение. Поэтому применение стандартных процедур, требующих нескольких последовательных итераций балансировки, может приводить к значительному увеличению продолжительности испытаний.
Проведение натурных испытаний вертолетной техники связано с высокой стоимостью эксплуатации и обслуживания летательного аппарата, расходом ресурса агрегатов, необходимостью привлечения испытательного персонала и обеспечением специальных требований безопасности. Кроме того, каждый дополнительный режим работы и каждое изменение конфигурации несущей системы требуют повторного проведения измерений. Поэтому актуальной является задача разработки методов идентификации дисбаланса, обеспечивающих требуемую точность при минимальном количестве экспериментальных прогонов.
С этой точки зрения перспективным является применение методов параметрической идентификации, основанных на совместном анализе вибрационного отклика системы и математической модели её динамики [22, 23]. Использование физически обоснованной модели позволяет связать измеряемые перемещения и вибрации с параметрами дисбаланса и определить их значения непосредственно по экспериментальным данным. При этом постановка задачи оптимизации позволяет одновременно учитывать вклад нескольких параметров и использовать информацию, содержащуюся во всем измерительном интервале, вместо последовательного изменения отдельных параметров.
Для соосной несущей системы особенно важным является учёт противоположных направлений вращения роторов, их фазовых характеристик и гироскопического взаимодействия [24, 25]. При этом использование нескольких пространственно разнесенных точек измерения позволяет восстановить не только локальный вибрационный отклик, но и пространственное движение конструкции. Такая информация может быть использована для построения более информативной модели идентификации параметров дисбаланса [26].
Таким образом, актуальной научно-технической задачей является разработка метода определения параметров дисбаланса соосных роторов по измеренному пространственному вибрационному отклику, обеспечивающего высокую точность оценки при минимальном количестве экспериментальных прогонов. Решение данной задачи позволяет снизить стоимость и продолжительность испытаний, а также создать основу для последующей автоматизации процедур диагностики и балансировки несущих систем вертолётов.
2 Определение модели
Положение системы описывается относительно двух систем координат (Рис. 1): неподвижной \(OXYZ\) и связанной с несущей конструкцией \(O'X'Y'Z'\), начало которой совпадает с центром вращения \(O'\) (центром масс конструкции). Описание пространственного движения системы с использованием обобщённых координат соответствует общепринятому подходу к моделированию многотельных механических систем [27]. Ось вращения обоих роторов совпадает с осью \(Z'\) связанной системы координат. Координаты вектора \(x\), \(y\), \(z\) задают положение точки \(O'\) в неподвижной системе координат, а углы \(\varphi\), \(\theta\), \(\psi\) — ориентацию связанной системы относительно неподвижной. Обобщённые координаты описывают перемещения центра вращения \(q_t\) и малые угловые отклонения системы \(q_r\):
\[ q = \begin{bmatrix} q_t & q_r \end{bmatrix}^{\top},\qquad q_t = \begin{bmatrix} x & y & z \end{bmatrix}^{\top},\qquad q_r = \begin{bmatrix} \varphi & \theta & \psi \end{bmatrix}^{\top}. \tag{1}\]
где \(\varphi\) — угол крена (угол поворота относительно оси \(X\)), \(\theta\) — угол тангажа (угол поворота относительно оси \(Y\)), \(\psi\) — угол рыскания (угол поворота относительно оси \(Z\)).
Система состоит из несущей конструкции (масса \(m_{\text{str}}\), тензор инерции \(J_{\text{str}}^0\)), верхнего и нижнего роторов (масса \(m_r\), тензоры инерции \(J_r^0\), начальные углы установки \(\beta_{0,\text{up}}\), \(\beta_{0,\text{down}}\)) и \(N_s\) упругих опор.
Конструкция опирается на \(N_s \ge 3\) упруго-вязких опор, расположенных в точках \(P_i\), \(i = 1,\ldots,N_s\). Для устойчивого положения несущей конструкции по всем шести степеням свободы необходимо не менее трёх опор, не лежащих на одной прямой. В недеформированном состоянии положение \(i\)-й опоры относительно центра вращения задаётся радиус-вектором \(r_i^0\). Опоры считаются безынерционными, а их жёсткость и демпфирование — диагональными по осям \(X\), \(Y\), \(Z\):
\[ k = \begin{bmatrix} k_{x,1} & k_{y,1} & k_{z,1} & \ldots & k_{x,N_s} & k_{y,N_s} & k_{z,N_s} \end{bmatrix}, \tag{2}\]
\[ c = \begin{bmatrix} c_{x,1} & c_{y,1} & c_{z,1} & \ldots & c_{x,N_s} & c_{y,N_s} & c_{z,N_s} \end{bmatrix}. \tag{3}\]
Дисбаланс каждого из винтов задаётся несимметричным распределением массы относительно оси вращения. Для описания дисбаланса вводятся масса неуравновешенного элемента \(m_i\), расстояние \(e_i\) от оси вращения до его центра масс и фаза дисбаланса \(\alpha_{0,i}\). Вектор параметров дисбаланса системы имеет вид
\[ p = \begin{bmatrix} m_{\text{up}} & e_{\text{up}} & \alpha_{0,\text{up}} & m_{\text{down}} & e_{\text{down}} & \alpha_{0,\text{down}} \end{bmatrix}, \tag{4}\]
где \(m_i\) — величина дисбалансирующей массы \(i\)-го винта, \(e_i\) — расстояние от оси вращения до центра масс дисбалансирующей массы, \(\alpha_{0,i}\) — фаза дисбаланса.
Частота вращения обоих роторов у вертолётов с соосной схемой синхронизированы и равны по модулю, а направления их вращения являются противоположными. В связи с этим текущие фазы дисбалансов изменяются с одинаковой по модулю, но противоположной по знаку угловой скоростью:
\[ \alpha(t) = \int\limits_0^t \omega(\tau)\,d\tau, \tag{5}\]
\[ \alpha_i(t) = \alpha_{0,i} + s_i\,\alpha(t), \tag{6}\]
где \(\omega(t)\) — общая для обоих роторов угловая скорость вращения по модулю в момент времени \(t\), а \(s_i\) — коэффициент, определяющий направление вращения ротора (\(s_i = -1\) при вращении \(i\)-го ротора по часовой стрелке и \(s_i = 1\) в противоположном случае).
Таким образом, положение дисбалансирующей массы относительно оси вращения определяется радиус-вектором
\[ \rho_i(t) = \begin{bmatrix} e_i \cos\alpha_i(t) \\ e_i \sin\alpha_i(t) \\ 0 \end{bmatrix}. \tag{7}\]
Динамика системы описывается линейной матричной моделью второго порядка [28]:
\[ M(\alpha(t))\,\ddot q + G(p,\omega)\,\dot q + \dot M(\alpha(t),\omega(t))\,\dot q = Q_{\text{dis}}\big(q,p,\alpha(t),\omega(t)\big) + Q_{\text{sup}}(q,\dot q), \tag{8}\]
где \(M(\alpha(t))\) — матрица масс системы (см. п. 2.1), \(\dot M(\alpha(t),\omega(t))\) — производная матрицы масс по времени (см. п. 2.1.1), \(G(p,\omega(t))\) — гироскопическая матрица (см. п. 2.2), \(Q_{\text{dis}}(q,p,\alpha(t),\omega(t))\) — вектор обобщённых сил от дисбалансов роторов (см. п. 2.4), \(Q_{\text{sup}}(q,\dot q)\) — вектор обобщённых сил от упругих опор (см. п. 2.3).
2.1 Матрица масс M
Для рассматриваемой системы, имеющей поступательные и вращательные степени свободы, матрица масс имеет блочную структуру
\[ M(\alpha) = \begin{bmatrix} m I_3 & 0 \\ 0 & J(\alpha(t)) \end{bmatrix},\qquad m = m_{\text{str}} + 2m_r, \tag{9}\]
где \(m\) — суммарная масса системы без учёта дисбалансирующих масс, \(I_3\) — единичная матрица третьего порядка, \(J(\alpha)\) — тензор инерции системы относительно точки \(O'\).
Тензор инерции системы формируется из инерционных характеристик несущей конструкции и двух роторов. Каждое тело (конструкция, верхний ротор, нижний ротор) \(k \in \{\text{str}, \text{up}, \text{down}\}\) задаётся массой \(m_k\), тензором инерции \(J_k\) относительно собственного центра масс и радиус-вектором центра масс \(r_k\) в системе координат сборки. Для роторов \(r_{\text{rot},i} = p_i + c_i\), где \(p_i\) — положение системы координат ротора, а \(c_i\) — положение его центра масс в этой системе.
Поскольку при вращении роторов их ориентация изменяется относительно конструкции, исходные тензоры инерции роторов преобразуются в текущую ориентацию:
\[ J_i(\alpha) = R_i(\alpha)\,J_r^0\,R_i^{\top}(\alpha), \qquad i \in \{\text{up}, \text{down}\}, \tag{10}\]
где \(R_i(\alpha)\) — матрица поворота вокруг оси вращения \(Z\):
\[ R_i(\alpha) = R_z\big(\beta_{0,i} + s_i\alpha\big),\qquad R_z(\beta) = \begin{bmatrix} \cos\beta & -\sin\beta & 0 \\ \sin\beta & \cos\beta & 0 \\ 0 & 0 & 1 \end{bmatrix}, \tag{11}\]
где \(\beta_{0,i}\) — начальный угол установки \(i\)-го ротора, определяется относительно опорного углового положения, соответствующего моменту формирования импульса датчиком частоты вращения. В используемой измерительной системе применяются импульсные датчики с одним импульсом за оборот, что позволяет фиксировать момент прохождения ротора через заданное угловое положение и определять фазу.
Поворот происходит вокруг оси ротора, проходящей через его центр масс, поэтому положение центра масс ротора при вращении не меняется, а изменяется только его тензор инерции. Для роторов с тройной симметрией лопастей тензор \(J_i\) не зависит от угла поворота.
Для учёта различного положения центров масс тел относительно точки, к которой приводится тензор, используется теорема Гюйгенса–Штейнера. Перенос тензора инерции тела из его центра масс \(C\) в точку \(O\) определяется выражением
\[ J_O = J_C + m\left( d^{\top} d\, I_3 - d\,d^{\top} \right),\qquad d = r_O - r_C. \tag{12}\]
Формирование тензора \(J(\alpha)\) выполняется в два этапа.
Сборка. Определяются суммарная масса и общий центр масс системы:
\[ m = \sum_k m_k,\qquad r_c = \frac{1}{m}\sum_k m_k r_k, \tag{13}\]
а тензоры инерции всех тел переносятся в общий центр масс и суммируются:
\[ J_c(\alpha) = \sum_k \left( J_k + m_k \left( d_k^{\top} d_k\, I_3 - d_k d_k^{\top} \right) \right),\qquad d_k = r_k - r_c, \tag{14}\]
где \(J_{\text{str}}\) — постоянный тензор инерции несущей конструкции, \(J_{\text{up}} = J_{\text{up}}(\alpha)\), \(J_{\text{down}} = J_{\text{down}}(\alpha)\).
Перенос. Полученный тензор переносится из общего центра масс \(r_c\) в центр вращения \(q_t\):
\[ J(\alpha) = J_c(\alpha) + m\left( d_0^{\top} d_0\, I_3 - d_0 d_0^{\top} \right),\qquad d_0 = r_c - q_t. \tag{15}\]
Таким образом, матрица инерции \(J(\alpha)\) изменяется во времени вследствие изменения ориентации роторов. В рассматриваемой модели дисбалансирующие массы не включаются в состав тел, формирующих матрицу инерции (они не участвуют в определении \(m\), \(r_c\) и \(J\)), а их влияние учитывается отдельно в векторе обобщённых сил.
2.1.1 Производная матрицы масс \(\dot M\)
Матрица масс зависит от времени только через фазу вращения роторов \(\alpha(t)\), поэтому её производная по времени определяется как
\[ \dot M(\alpha,\omega) = \frac{\partial M}{\partial \alpha}\,\dot\alpha = \omega(t)\,\frac{\partial M}{\partial \alpha}. \tag{16}\]
Поступательный блок \(mI_3\) постоянен, так как масса системы без дисбалансирующих масс не меняется. Поэтому производная имеет ненулевой только вращательный блок:
\[ \dot M(\alpha,\omega) = \begin{bmatrix} 0_{3\times3} & 0_{3\times3} \\ 0_{3\times3} & \dot J(\alpha,\omega) \end{bmatrix},\qquad \dot J(\alpha,\omega) = \omega\,\frac{\partial J}{\partial \alpha}. \tag{17}\]
Здесь \(\dot M\) имеет размер \(6\times6\), соответствующий вектору обобщённых координат \(q\). Тензор \(J(\alpha)\) вычисляется методом центральных конечных разностей. Для малого шага по времени \(\Delta t\):
\[ \dot J(\alpha,\omega) \approx \frac{J\big(\alpha + \omega(t)\Delta t\big) - J\big(\alpha - \omega(t)\Delta t\big)}{2\Delta t}. \tag{18}\]
Для роторов с тройной симметрией лопастей тензоры \(J_{\text{up}}\) и \(J_{\text{down}}\) не зависят от угла поворота. В этом случае \(J(\alpha) = \text{const}\) и \(\dot J = 0\), а \(\dot M\) можно не учитывать.
В рамках рассматриваемой модели принимается, что инерционные характеристики несущей конструкции и роторов постоянны, а изменение результирующего тензора инерции обусловлено только изменением ориентации роторов относительно связанной системы координат. Деформации лопастей, изменение их конфигурации и перераспределение масс в процессе вращения не учитываются.
Дальнейшее развитие модели может предусматривать учёт переменных инерционных характеристик, зависящих от дополнительных параметров системы. К ним относятся масса и пространственное распределение дополнительной нагрузки, деформации лопастей под действием центробежных сил, а также изменение их конфигурации при изменении угла установки или угла атаки. В общем случае такие эффекты могут быть представлены зависимостью тензора инерции от вектора параметров и текущего состояния системы:
\[ J = J(\alpha, q, \dot q, \lambda), \tag{19}\]
где \(\lambda\) — совокупность дополнительных параметров, характеризующих конфигурацию и инерционные свойства системы. Введение такой зависимости позволит расширить модель для описания динамики системы с учётом изменения конфигурации её элементов и перераспределения массы.
2.2 Гироскопическая матрица G
Гироскопические эффекты вызываются вращением двух контрвращающихся роторов с дисбалансами. Роторы вращаются вокруг оси \(Z'\) корпуса, поэтому их суммарный вектор кинетического момента направлен вдоль этой оси:
\[ G(p,\omega) = \begin{bmatrix} 0_{3\times3} & 0_{3\times3} \\ 0_{3\times3} & G_{\text{rot}} \end{bmatrix},\qquad G_{\text{rot}} = \begin{bmatrix} 0 & H_{\text{sum}} & 0 \\ -H_{\text{sum}} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}. \tag{20}\]
Кинетический момент \(i\)-го ротора относительно оси вращения равен \(H_i = s_i I_i^z \omega\). Осевой момент инерции ротора с дисбалансом, представленным сосредоточенной массой \(m_i\) на расстоянии \(e_i\) от оси, по теореме Гюйгенса–Штейнера равен:
\[ I_i^z = J_{r}^{z} + m_i e_i^2,\qquad i \in \{\text{up}, \text{down}\}, \tag{21}\]
где \(J_{r}^{z}\) — момент инерции ротора без дисбалансирующей массы относительно оси вращения. Суммарный кинетический момент двух роторов из (21):
\[ H_{\text{sum}} = s_{\text{up}}H_{\text{up}} + s_{\text{down}}H_{\text{down}} = \omega\left( s_{\text{up}} m_{\text{up}} e_{\text{up}}^2 + s_{\text{down}} m_{\text{down}} e_{\text{down}}^2 \right). \tag{22}\]
При равных дисбалансах (20):
\[ H_{\text{sum}} = 0,\qquad G = 0_{6\times6}. \tag{23}\]
Для неидентичных роторов гироскопический момент также будет определяться разностью инерций роторов: \(J_{\text{up}}^z - J_{\text{down}}^z\).
2.3 Упругие опоры. Матрицы демпфирования C и жёсткости K
Несущая система опирается на \(N_s\) упругих опор. Каждая \(i\)-я опора является упруго-вязкой связью с диагональными матрицами жёсткости и демпфирования:
\[ K_i = \operatorname{diag}\left( k_{x,i}, k_{y,i}, k_{z,i} \right),\qquad C_i = \operatorname{diag}\left( c_{x,i}, c_{y,i}, c_{z,i} \right),\qquad i = 1,\ldots,N_s. \tag{24}\]
Положение \(i\)-й опоры в недеформированном состоянии задаётся радиус-вектором \(r_i^0 = r_{S,i} - q_t\), проведённым из центра вращения \(q_t\) в точку крепления опоры \(r_{S,i}\). При повороте системы на углы \(\varphi\), \(\theta\), \(\psi\) он принимает вид:
\[ r_i = R(\varphi,\theta,\psi)\,r_i^0, \tag{25}\]
где \(R\) — матрица поворота корпуса, задаётся последовательностью элементарных поворотов вокруг осей \(X\), \(Y\), \(Z\) неподвижной системы координат:
\[ R(\varphi,\theta,\psi) = R_z(\psi)\,R_y(\theta)\,R_x(\varphi), \tag{26}\]
\[ R_x(\varphi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\varphi & -\sin\varphi \\ 0 & \sin\varphi & \cos\varphi \end{bmatrix},\qquad R_y(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix},\qquad R_z(\psi) = \begin{bmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{bmatrix}. \]
Перемещение и скорость точки крепления опоры равны
\[ u_i = q_t + r_i - r_i^0,\qquad v_i = \dot q_t + \dot q_r \times r_i, \tag{27}\]
где \(\dot q_r\) принимаются за вектор угловой скорости корпуса (допустимо при малых углах). Сила, действующая на систему со стороны \(i\)-й опоры, и её момент относительно центра вращения:
\[ F_i = -K_i u_i - C_i v_i,\qquad M_i = r_i \times F_i. \tag{28}\]
Вклад всех опор в вектор обобщённых сил равен
\[ Q_{\text{sup}} = \sum_{i=1}^{N_s} \begin{bmatrix} F_i \\ r_i \times F_i \end{bmatrix}. \tag{29}\]
2.4 Вектор обобщённых сил от дисбаланса \(Q_{\text{dis}}\)
Основным возмущающим воздействием являются центробежные силы, возникающие из-за дисбаланса верхнего и нижнего роторов [13, 29]. Дисбалансирующая масса \(m_i\), расположенная на расстоянии \(e_i\) от оси вращения, при угловой скорости \(\omega\) создаёт силу, модуль которой равен
\[ F_i(t) = m_i e_i \omega(t)^2,\qquad i \in \{\text{up}, \text{down}\}. \tag{30}\]
Для характеристики дисбаланса каждого ротора введём статический момент дисбаланса
\[ U_i = m_i e_i. \tag{31}\]
Сила направлена вдоль радиус-вектора \(\rho_i(t)\) (от оси вращения к дисбалансирующей массе). В системе координат, связанной с ротором, её направление задаётся единичным вектором (7):
\[ n_i^{\text{loc}}(t) = \begin{bmatrix} \cos\alpha_i(t) \\ \sin\alpha_i(t) \\ 0 \end{bmatrix}. \tag{32}\]
Системы координат роторов ориентированы вместе с корпусом, поэтому при его повороте на углы \(\varphi\), \(\theta\), \(\psi\) направление силы в глобальной системе координат получается умножением на матрицу поворота \(R = R(\varphi,\theta,\psi)\) (см. (26) и (32)):
\[ F_i(t) = F_i(t)\,R\,n_i^{\text{loc}}(t) = m_i e_i \omega^2(t)\, R \begin{bmatrix} \cos\alpha_i(t) \\ \sin\alpha_i(t) \\ 0 \end{bmatrix}. \tag{33}\]
Силы дисбаланса приложены к осям роторов. Плечо \(i\)-й силы относительно центра вращения \(O'\), положение которого задаётся поступательными обобщёнными координатами \(q_t\), определяется выражениями
\[ r_{\text{rot},i} = R\, r_{\text{rot},i}^0, \tag{34}\]
где \(r_{\text{rot},i}^0\) — радиус-вектор начала системы координат \(i\)-го ротора на оси вращения относительно центра \(O'\) в недеформированном состоянии.
Суммарная сила и суммарный момент от двух роторов равны
\[ F_{\text{dis}} = F_{\text{up}} + F_{\text{down}},\qquad M_{\text{dis}} = r_{\text{rot,up}} \times F_{\text{up}} + r_{\text{rot,down}} \times F_{\text{down}}. \tag{35}\]
Вектор обобщённых сил от дисбаланса, соответствующий вектору \(q\), имеет вид
\[ Q_{\text{dis}}(q,p,\alpha,\omega) = \begin{bmatrix} F_{\text{dis}} \\ M_{\text{dis}} \end{bmatrix} = \begin{bmatrix} F_x & F_y & F_z & M_\varphi & M_\theta & M_\psi \end{bmatrix}^{\top}, \tag{36}\]
где \(F_x\), \(F_y\), \(F_z\) — проекции \(F_{\text{dis}}\) на глобальные оси, а \(M_\varphi\), \(M_\theta\), \(M_\psi\) — проекции \(M_{\text{dis}}\) на оси \(X\), \(Y\), \(Z\).
В принятой модели не учитываются: тангенциальная составляющая силы при \(\dot\omega \ne 0\), момент \(\rho_i \times F_i\) из-за смещения дисбалансирующей массы относительно оси ротора и смещение центра масс системы из-за дисбалансирующих масс (они не входят в \(m\), \(r_c\) и \(J\), см. п. 2.1).
3 Моделирование движения системы
Прямая задача состоит в определении движения \(q(t)\) (1) при известных параметрах модели \(k\), \(c\), \(p\) (2), (3), (4) и заданной частоте вращения роторов \(\omega(t)\). Решение используется для проверки модели и методов обработки на данных с известными параметрами.
3.1 Приведение к системе первого порядка
Уравнение (8) является системой шести обыкновенных дифференциальных уравнений второго порядка. Для численного интегрирования оно приводится к системе первого порядка относительно вектора состояния
\[ y(t) = \begin{bmatrix} q(t) \\ \dot q(t) \end{bmatrix} \in \mathbb{R}^{12}:\qquad \dot y = \begin{bmatrix} \dot q \\ \ddot q \end{bmatrix} = f(t,y). \tag{37}\]
Матрица масс \(M(\alpha)\) положительно определена, поэтому ускорение определяется из уравнения (8):
\[ \ddot q = M^{-1}(\alpha)\,F(q,\dot q,t), \tag{38}\]
\[ F = Q_{\text{dis}}(q,p,\alpha,\omega) + Q_{\text{sup}}(q,\dot q) - \big( G(m,e,\omega) + \dot M(\alpha,\omega) \big)\,\dot q. \tag{39}\]
3.2 Вычисление правой части
Правая часть \(f(t,y)\) вычисляется в такой последовательности.
Определяется фаза вращения роторов (5). Для постоянной частоты \(\omega\):
\[ \alpha_i(t) = s_i \pm \omega t \tag{40}\]
Строится матрица масс \(M(\alpha)\) (см. п. 2.1).
Вычисляется вектор \(F\): \(Q_{\text{dis}}\) по п. 2.4, \(Q_{\text{sup}}\) по п. 2.3 (с точной матрицей поворота \(R(\varphi,\theta,\psi)\) для текущего \(q\)), а также слагаемое \((G + \dot M)\dot q\).
Из линейной системы \(M\ddot q = F\) находится \(\ddot q\), после чего
\[ f(t,y) = \begin{bmatrix} \dot q \\ \ddot q \end{bmatrix}. \tag{41}\]
3.3 Начальные условия
Начальное состояние задаётся вектором \(y_0 = [q_0^{\top}, \dot q_0^{\top}]^{\top}\). По умолчанию система находится в покое в положении равновесия:
\[ q_0 = 0,\qquad \dot q_0 = 0. \tag{42}\]
Сила тяжести в уравнение (8) не входит, так как \(q\) отсчитывается от положения статического равновесия под действием веса, поэтому \(q_0 = 0\) соответствует именно ему. Роторы запускаются мгновенно с частотой \(\omega\), из-за чего в решении возникает переходный процесс, затухающий за счёт демпфирования опор. Для получения установившегося режима моделируется промежуток времени, заведомо превышающий время затухания.
3.4 Численное интегрирование
Задача Коши \(\dot y = f(t,y)\), \(y(0) = y_0\) решается на отрезке \([0, T_{\text{sim}}]\) явным методом Рунге–Кутты Дормана–Принса [30] с автоматическим выбором шага (реализация SciPy [31]). Результат записывается в узлах равномерной сетки
\[ t_k = k\,\Delta t,\qquad k = 0,\ldots,N_{\text{sim}},\qquad N_{\text{sim}} = \operatorname{round}\left( T_{\text{sim}} / \Delta t \right), \tag{43}\]
где \(\Delta t = 1/f_s\) — шаг дискретизации, совпадающий с периодом опроса датчиков, чтобы модельные и измеренные сигналы можно было обрабатывать одинаково. Шаг интегрирования ограничивается сверху значением \(\Delta t\).
На каждом шаге локальная погрешность оценивается разностью решений четвёртого и пятого порядков, и шаг принимается, если
\[ \left| e_i \right| \le \text{atol} + \text{rtol}\left| y_i \right|,\qquad i = 1,\ldots,12, \tag{44}\]
где \(\text{rtol}\) и \(\text{atol}\) — относительный и абсолютный допуски. Для устойчивости явного метода необходимо, чтобы произведение шага на максимальную собственную частоту системы было мало (\(\Delta t\,\omega_{\max} \lesssim 3\)), что при частоте опроса порядка нескольких килогерц для рассматриваемых опор выполняется с запасом.
Полученные \(q(t_k)\) и \(\dot q(t_k)\) далее используются для проверки обратной задачи; ускорение \(\ddot q(t_k)\) вычисляется из (38) по найденному решению, а не численным дифференцированием, поэтому оно не содержит шума дифференцирования.
Пример смоделированного сигнала представлен на Рис. 2 и Рис. 3. Параметры дисбаланса модели: \(m_{\text{up}} = m_{\text{down}} \ne 0\), \(e_{\text{up}} = e_{\text{down}} \ne 0\), \(\alpha_{0,\text{up}} = \pi/2\), \(\alpha_{0,\text{down}} = \pi/2\), \(\omega = 4.7\) Гц.
Спектр сигнала на установившемся режиме представлен на Рис. 4. На спектрах видны гармоники частоты вращения роторов (1–3 роторные), а также собственные частоты конструкции (≈ 0.7 Гц).
4 Определение обобщённых координат по измерениям вибродатчиков
Обратная задача (гл. 5) требует знания \(q(t)\), \(\dot q(t)\) и \(\ddot q(t)\) на установившемся режиме. В эксперименте измеряются траектории трёх точек конструкции, поэтому сначала по ним восстанавливается \(q(t)\), затем вычисляются производные. Те же процедуры применяются к результатам моделирования (гл. 3), что позволяет проверить метод на данных с известными параметрами.
4.1 Исходные данные
Используются три трёхосных вибродатчика, установленные в точках \(D_i^{\text{sens}}\), \(i = 1, 2, 3\), заданных в неподвижной системе координат в начальном (недеформированном) положении. Датчики дают смещения точек \(u_i(t_k)\) относительно этого положения на равномерной сетке \(t_k = k\Delta t\), \(\Delta t = 1/f_s\). Положение точек в произвольный момент:
\[ D_i(t) = D_i^{\text{sens}} + u_i(t),\qquad i = 1, 2, 3. \tag{45}\]
Точки не должны лежать на одной прямой:
\[ \left( D_2 - D_1 \right) \times \left( D_3 - D_1 \right) \ne 0. \tag{46}\]
Это условие однозначно определяет положение конструкции как твёрдого тела: три неколлинеарные точки фиксируют и положение, и ориентацию. Для начала отсчёта принимается первый отсчёт записи: \(D_i^0 = D_i(t_0)\). Углы и перемещения, найденные ниже, отсчитываются от этой конфигурации.
4.2 Восстановление координат по измеренным ускорениям
В случае применения трёхосных акселерометров исходными измеряемыми величинами являются ускорения контрольных точек конструкции. Для определения обобщённых координат системы необходимо восстановить траектории этих точек, а затем определить положение и ориентацию конструкции как твёрдого тела. Такой подход обусловлен нелинейной зависимостью координат точек от углов ориентации, поэтому интегрирование выполняется для линейных координат точек, а не непосредственно для угловых обобщённых координат.
Пусть \(a_i(t)\) — измеренный вектор ускорения \(i\)-й точки, \(i = 1, 2, 3\). Восстановление её траектории требует решения задачи двойного интегрирования:
\[ u_i(t) = u_i(t_0) + \dot u_i(t_0)\left( t - t_0 \right) + \int\limits_{t_0}^{t} \int\limits_{t_0}^{\tau} a_i(\xi)\, d\xi\, d\tau, \tag{47}\]
где \(u_i(t)\) — перемещение точки относительно начального положения, \(u_i(t_0)\) и \(\dot u_i(t_0)\) — начальные перемещение и скорость. Для отсчёта перемещений от начальной конфигурации принимается \(u_i(t_0) = 0\) и \(\dot u_i(t_0) = 0\).
Прямое двойное интегрирование измеренных ускорений приводит к накоплению ошибок [32]. В частности, постоянная составляющая смещения нуля акселерометра \(b\) вызывает ошибку восстановленного перемещения, пропорциональную \(b(t - t_0)^2/2\). Поэтому для повышения устойчивости восстановления траекторий необходимо компенсировать смещение нуля, ограничивать низкочастотную составляющую ошибки и использовать методы регуляризации.
Для установившихся периодических колебаний траектории могут восстанавливаться в частотной области [33]. При представлении ускорения рядом Фурье
\[ a_i(t) = \sum_{k=-\infty}^{\infty} A_{i,k}\, e^{ik\omega t}, \tag{48}\]
коэффициенты Фурье перемещений для ненулевых гармоник определяются выражением
\[ U_{i,k} = -\frac{A_{i,k}}{(k\omega)^2},\qquad k \ne 0, \tag{49}\]
где \(\omega\) — основная угловая частота колебаний, \(k\) — номер гармоники, \(A_{i,k}\) и \(U_{i,k}\) — комплексные векторы коэффициентов Фурье ускорения и перемещения соответственно. Восстановление выполняется по гармоникам, соответствующим целому числу периодов записи. Постоянная составляющая перемещения при таком подходе не определяется ускорением и задаётся отдельно выбором начала отсчёта.
После восстановления траекторий трёх контрольных точек их положения определяются как
\[ D_i(t) = D_i^{\text{sens}} + u_i(t),\qquad i = 1, 2, 3. \tag{50}\]
4.3 Определение ориентации (алгоритм Кабша)
Для исключения поступательного движения вводятся геометрический центр точек и центрированные координаты:
\[ g(t) = \frac{1}{3}\sum_{i=1}^{3} D_i(t),\qquad s_i(t) = D_i(t) - g(t),\qquad s_i^0 = D_i^0 - g_0,\qquad g_0 = g(t_0). \tag{51}\]
Конструкция считается абсолютно твёрдой, поэтому \(s_i(t) \approx R(t)\,s_i^0\). Матрица поворота \(R(t) \in SO(3)\) находится из задачи наименьших квадратов (задача Кабша, или задача Вахбы [34, 35])
\[ R(t) = \arg \min_{R \in SO(3)} \sum_{i=1}^{3} \left\| s_i(t) - R\,s_i^0 \right\|^2. \tag{52}\]
Решение строится через корреляционную матрицу и её сингулярное разложение:
\[ H(t) = \sum_{i=1}^{3} s_i^0\, s_i^{\top}(t) = Q_0^{\top} Q(t) = U\Sigma V^{\top},\qquad R(t) = V \operatorname{diag}\left( 1, 1, d \right) U^{\top}, \tag{53}\]
где строки \(Q_0\) и \(Q(t)\) составлены из \(s_i^{0\top}\) и \(s_i^{\top}(t)\), а \(d = \operatorname{sign} \det\left( VU^{\top} \right)\). Множитель \(d\) исключает отражение и обеспечивает \(\det R = 1\) (при \(\det VU^{\top} < 0\) меняется знак последнего сингулярного вектора). Три центрированные точки лежат в одной плоскости, поэтому \(\operatorname{rank} H = 2\), третье сингулярное число равно нулю, и именно поправка \(d\) однозначно определяет третий столбец \(R\).
4.4 Обобщённые координаты
Центр вращения \(O'\) жёстко связан с конструкцией. Его положение в начальный момент \(r_c\) известно из конструкции. Вектор от геометрического центра точек до него постоянен в связанной системе:
\[ r_0 = r_c - g_0. \tag{54}\]
Положение центра вращения в момент \(t\) и поступательные координаты:
\[ P_0(t) = g(t) + R(t)\,r_0,\qquad \begin{bmatrix} x \\ y \\ z \end{bmatrix} = P_0(t) - r_c. \tag{55}\]
Углы определяются из разложения \(R = R_z(\psi)R_y(\theta)R_x(\varphi)\) (см. п. 2.3):
\[ \varphi = \operatorname{atan2}\left( R_{32}, R_{33} \right),\qquad \theta = \operatorname{atan2}\left( -R_{31}, \sqrt{R_{11}^2 + R_{21}^2} \right),\qquad \psi = \operatorname{atan2}\left( R_{21}, R_{11} \right). \tag{56}\]
Углы малы (\(|\theta| \ll \pi/2\)), особенность \(\theta = \pm\pi/2\) не достигается, а перескоков через \(\pm\pi\) нет, поэтому развёртка не требуется. В результате формируется вектор
\[ q(t_k) = \begin{bmatrix} x & y & z & \varphi & \theta & \psi \end{bmatrix}^{\top},\qquad k = 0, 1, \ldots \]
4.5 Вычисление скоростей \(\dot q(t)\) и ускорений \(\ddot q(t)\)
Для идентификации параметров динамической модели необходимо знать обобщённые координаты системы, их первые и вторые производные по времени. Вектор обобщённых координат \(q(t)\) включает три поступательные координаты и три угловые координаты:
\[ q(t) = \begin{bmatrix} x(t) & y(t) & z(t) & \varphi(t) & \theta(t) & \psi(t) \end{bmatrix}^{\top}. \]
Пусть измеренные значения обобщённых координат заданы на равномерной временной сетке
\[ t_k = k\Delta t,\qquad \Delta t = \frac{1}{f_s}, \]
где \(k = 0,\ldots,N-1\) — индекс отсчёта, \(N\) — число отсчётов, \(f_s\) — частота дискретизации, \(\Delta t\) — шаг по времени.
По результатам измерений доступны перемещения конструкции, тогда как векторы обобщённых скоростей \(\dot q(t)\) и ускорений \(\ddot q(t)\) необходимо вычислить численно. Прямое дифференцирование измеренных сигналов чувствительно к высокочастотному шуму, особенно при вычислении второй производной. Поэтому для оценки скоростей и ускорений применяется дифференцирующий фильтр Савицкого–Голея [36], основанный на локальной полиномиальной аппроксимации сигнала методом наименьших квадратов. Для каждой компоненты \(q_k = q(t_k)\) в окрестности отсчёта \(k\) рассматривается окно из \(W = 2m + 1\) точек, где \(m\) — число отсчётов с каждой стороны от центральной точки. На этом окне сигнал аппроксимируется полиномом степени \(n\):
\[ p_k(\tau) = \sum_{l=0}^{n} a_{l,k}\,\tau^l, \tag{57}\]
где \(\tau\) — локальное время относительно центрального отсчёта, \(a_{l,k}\) — коэффициенты полинома, определяемые из условия минимизации суммы квадратов отклонений:
\[ \min_{a_{0,k},\ldots,a_{n,k}} \sum_{j=-m}^{m} \left[ q_{k+j} - \sum_{l=0}^{n} a_{l,k} \left( j\Delta t \right)^l \right]^2. \tag{58}\]
Здесь \(q_{k+j}\) обозначает значение одной выбранной компоненты вектора обобщённых координат в отсчёте \(k + j\). Коэффициенты полинома зависят от выбранной компоненты и центрального отсчёта \(k\).
Оценки первой и второй производных в центральной точке получают дифференцированием аппроксимирующего полинома:
\[ \dot q_k = a_{1,k},\qquad \ddot q_k = 2a_{2,k}. \tag{59}\]
Таким образом, вычисление производных сводится к определению коэффициентов локального полинома. Поскольку аппроксимация линейна относительно измеренных значений, обе производные могут быть представлены в виде линейных комбинаций отсчётов:
\[ \dot q_k = \frac{1}{\Delta t}\sum_{j=-m}^{m} c_j^{(1)}\, q_{k+j},\qquad \ddot q_k = \frac{1}{\Delta t^2}\sum_{j=-m}^{m} c_j^{(2)}\, q_{k+j}, \tag{60}\]
где \(c_j^{(1)}\) и \(c_j^{(2)}\) — безразмерные коэффициенты фильтра для первой и второй производных соответственно. Они определяются размером окна \(W\) и степенью полинома \(n\).
Выбор параметров фильтра представляет собой компромисс между подавлением шума и сохранением динамических составляющих сигнала. Увеличение длины окна улучшает сглаживание, но может приводить к ослаблению амплитуды и искажению фазовых характеристик быстро меняющихся составляющих. Поэтому длительность окна \(W\Delta t\) выбирается существенно меньше характерного периода исследуемых колебаний. Вблизи начала и конца записи симметричное окно недоступно, что приводит к краевым эффектам. Для уменьшения их влияния при идентификации параметров используется внутренний участок записи, исключающий крайние отсчёты, обработанные с односторонней или иной краевой аппроксимацией.
5 Идентификация параметров модели
5.1 Постановка задачи
По измеренным \(q_k\), \(\dot q_k\), \(\ddot q_k\), фазе \(\alpha_k\) и частоте \(\omega_k\) (\(k = 1,\ldots,N\)) требуется найти параметры (2), (3), (4):
\[ \lambda = \begin{bmatrix} p & c & k \end{bmatrix}^{\top} \in \mathbb{R}^{n_\lambda},\qquad n_\lambda = 6 + 2(N_s + 2). \tag{61}\]
Остаток из уравнения движения (8) в момент \(t_k\):
\[ r_k(\lambda) = M(\alpha_k)\ddot q_k + G(m,e,\omega_k)\,\dot q_k + \dot M(\alpha_k,\omega_k)\,\dot q_k - Q_{\text{dis}}(q_k,p,\alpha_k,\omega_k) - Q_{\text{sup}}(q_k,\dot q_k; c, k) \in \mathbb{R}^6. \tag{62}\]
Здесь левая часть вычисляется по измеренным движениям, а правая зависит от параметров. При точной модели и точных данных \(r_k = 0\).
Задача оценки параметров — нелинейный метод наименьших квадратов с ограничениями:
\[ \hat\lambda = \arg \min_{\lambda_{\min} \le \lambda \le \lambda_{\max}} \frac{1}{2}\sum_{k=1}^{N} \left\| W r_k(\lambda) \right\|^2, \tag{63}\]
где \(W = \operatorname{diag}\left( w_1,\ldots,w_6 \right)\) — весовая матрица, приводящая компоненты остатка к сопоставимому масштабу (силы в Н и моменты в Н·м имеют разные порядки). Веса берутся обратно пропорциональными типичным амплитудам компонент \(M\ddot q\) по измерениям. Границы \(\lambda_{\min}\), \(\lambda_{\max}\) задают физически допустимые значения (\(m_i, e_i, c, k \ge 0\)), а также диапазон начальных фаз \(\alpha_{0,i} \in [-\pi, \pi)\).
5.2 Структура задачи и идентифицируемость
Параметры опор входят в \(Q_{\text{sup}}\) линейно: \(Q_{\text{sup}} = -K q - C\dot q\) при малых углах, то есть линейно по \(k\) и \(c\).
Параметры дисбаланса входят в \(Q_{\text{dis}}\) через комбинации
\[ u_i = m_i e_i \cos\alpha_{0,i},\qquad v_i = m_i e_i \sin\alpha_{0,i}, \tag{64}\]
линейно (при заданных \(q_k\) и \(\alpha_k\)). Таким образом, силы определяют только произведение \(m_i e_i\) и фазу \(\alpha_{0,i}\), а \(m_i\) и \(e_i\) по отдельности входят только в \(G\) через \(m_i e_i^2\). Этот вклад мал (\(m_i e_i^2 \ll J_i^z\), п. 2.2) и в пределе вообще не определяет \(m_i\) и \(e_i\) раздельно. Поэтому практически корректной оценкой являются \(m_i e_i\) и \(\alpha_{0,i}\), а разделение \(m_i\) и \(e_i\) возможно только при известном одном из них или при априорных границах.
Нелинейность остатка обусловлена тригонометрической зависимостью от \(\alpha_{0,i}\), вкладом \(m_i e_i^2\) в \(G\) и точной матрицей поворота в \(Q_{\text{sup}}\). Поэтому целевая функция может иметь несколько локальных минимумов, главным образом по фазам \(\alpha_{0,i}\).
Различимость дисбалансов. Верхний и нижний роторы вращаются с одной частотой в противоположных направлениях, поэтому их вклады имеют разные зависимости от \(\alpha\) (\(\cos(\alpha_{0,\text{up}} + \alpha)\) и \(\cos(\alpha_{0,\text{down}} - \alpha)\)), и по составляющим на частоте вращения оба определяются. Для оценки же пары роторов нужны обе плоскости сил/моментов (6 компонент), в противном случае дисбалансы вырождаются в одну эквивалентную силу.
5.3 Алгоритм решения
Задача (63) решается методом доверительной области (Trust Region Reflective, TRF [37]) для задач наименьших квадратов с ограничениями (реализация SciPy [31]). На итерации решается линеаризованная задача
\[ \min_{\delta\lambda} \left\| W r + W J\,\delta\lambda \right\|^2 \quad \text{при} \quad \left\| D\,\delta\lambda \right\| \le \Delta, \tag{65}\]
где \(J = \partial r / \partial \lambda\) — якобиан (вычисляется конечными разностями), \(\Delta\) — радиус доверительной области, \(D\) — диагональная матрица масштабов. Масштабы параметров берутся по столбцам якобиана (x_scale='jac'), так как параметры имеют разные размерности и порядки величин (\(m_i\) в кг, \(e_i\) в м, \(k\) в Н/м).
Для ослабления влияния локальных минимумов и ускорения сходимости расчёт проводится в два этапа.
- Грубый этап. Начальное приближение \(\lambda^{(0)}\) выбирается случайно внутри границ. Минимизация проводится по каждой второй точке выделенного периода (допуски \(10^{-8}\), не более 300 вычислений остатка).
- Уточнение. Найденное решение используется как начальное для расчёта по всем точкам периода (допуски \(10^{-10}\), не более 300 вычислений остатка).
Остановка происходит по относительному изменению функции, параметров или градиента (ftol, xtol, gtol). При необходимости глобальный поиск выполняется методом дифференциальной эволюции [38] в тех же границах, найденное решение затем уточняется методом TRF. Для повышения надёжности начальное приближение можно запускать несколько раз и выбирать решение с наименьшим значением целевой функции.
5.4 Оценка качества результата
Качество оценки определяется тремя способами.
Невязка. Сравниваются левая и правая части уравнения движения по каждой компоненте на выделенном периоде; относительная невязка
\[ \delta_j = \frac{\left\| r_{\cdot j} \right\|}{\left\| \left( M\ddot q \right)_{\cdot j} \right\|}. \tag{66}\]
Доверительные интервалы. В окрестности решения ковариация оценок приближённо
\[ \operatorname{cov}\left( \hat\lambda \right) \approx \hat\sigma^2 \left( J^{\top} W^{\top} W J \right)^{-1},\qquad \hat\sigma^2 = \frac{\sum_k \left\| W r_k \right\|^2}{6N - n_\lambda}. \tag{67}\]
Большие значения дисперсии отдельных параметров (например, \(m_i\), \(e_i\) порознь, п. 5.2) указывают на плохую идентифицируемость.
Проверка на моделировании. Для синтетических данных (гл. 3) с известными \(\lambda_{\text{true}}\) оценка сравнивается с истинным значением, а также сравниваются остатки при \(\hat\lambda\) и \(\lambda_{\text{true}}\) (последний остаток характеризует ошибку самих данных, получаемую при дифференцировании и обработке датчиков).
Дополнительно оценивается вычислительная эффективность алгоритма по времени идентификации \(T_{\text{id}}\) (в секундах), которое зависит от числа точек \(N\), числа идентифицируемых параметров \(n_\lambda\) и вычислительной системы. Условия расчёта и полученное значение приведены в разделе 6.
6 Результаты
Для проверки разработанного алгоритма выполнена идентификация параметров дисбаланса верхнего и нижнего роторов на синтетических данных, сформированных с использованием модели с известными параметрами \(\lambda_{\text{true}}\) [39]. В качестве исходных данных использованы временные реализации обобщённых координат и их производных, а оценка параметров осуществлялась путём минимизации невязки уравнений движения.
Для каждого ротора \(i \in \{\text{up}, \text{down}\}\) рассматривались масса дисбаланса \(m_i\), эксцентриситет \(e_i\) и начальная фаза \(\varphi_{0,i}\). Дополнительно определялся статический момент дисбаланса (31).
На Рис. 5 представлен пример сопоставления левой и правой частей уравнений движения, рассчитанных с идентифицированными параметрами, а также правой части, полученной при использовании истинных параметров модели. Совпадение расчётных зависимостей позволяет визуально оценить качество восстановления динамического отклика системы.
После определения параметров \(\hat\lambda\) выполнялось сравнение полученных оценок с истинными значениями \(\lambda_{\text{true}}\), использованными при формировании синтетических данных. Для количественной оценки точности определения статического момента дисбаланса рассчитывалась относительная погрешность
\[ \varepsilon_U^i = \frac{\left| \hat U_i - U_{i,\text{true}} \right|}{\left| U_{i,\text{true}} \right|} \cdot 100\%,\qquad i \in \{\text{up}, \text{down}\}. \tag{68}\]
Для начальной фазы учитывалась её периодичность. Абсолютная ошибка фазы определялась как
\[ \Delta\varphi_{0,i} = \left| \operatorname{atan2}\left( \sin(\hat\varphi_{0,i} - \varphi_{0,i,true}), \cos(\hat\varphi_{0,i} - \varphi_{0,i,true}) \right) \right|. \tag{69}\]
Исходные данные для идентификации получены с использованием модели вибрационного стенда, разработанной в программном комплексе Exudyn [40] и рассмотренной в работе. На основе сформированных временных реализаций движения системы выполнена идентификация параметров дисбаланса верхнего и нижнего роторов. Результаты идентификации представлены в Табл. 1.
| Параметр | Верхний ротор | Нижний ротор |
|---|---|---|
| Истинная масса \(m_{i,\text{true}}\), у.е. | 1.00 | 1.00 |
| Идентифицированная масса \(\hat m_i\), у.е. | 1.51 | 0.99 |
| Истинный эксцентриситет \(e_{i,\text{true}}\), у.е. | 1.00 | 1.00 |
| Идентифицированный эксцентриситет \(\hat e_i\), у.е. | 0.69 | 1.00 |
| Истинный статический момент \(U_{i,\text{true}}\), у.е. | 1.00 | 1.00 |
| Идентифицированный статический момент \(\hat U_i\), у.е. | 1.05 | 0.99 |
| Истинная начальная фаза \(\varphi_{0,i,\text{true}}\), ° | 0 | 0 |
| Идентифицированная начальная фаза \(\hat\varphi_{0,i}\), ° | −1.62 | 2.93 |
| Погрешность \(U_i\), % | 5 | 1 |
| Ошибка фазы \(\Delta\varphi_{0,i}\), ° | 1.62 | 2.93 |
Из Табл. 1 видно, что погрешность определения статического момента дисбаланса не превышает 5 %, а ошибка определения начальной фазы составляет менее 3°. При этом для верхнего ротора наблюдается заметное отклонение идентифицированных массы и эксцентриситета от истинных значений. Данное расхождение согласуется с особенностями идентифицируемости параметров, рассмотренными в подразделе 5.2. Для практической оценки дисбаланса это произведение и соответствующая начальная фаза являются основными характеристиками.
На Рис. 6 показана нормированная по амплитуде невязка уравнений движения при использовании идентифицированных параметров. Она позволяет оценить величину остаточного рассогласования относительно характерного масштаба каждой компоненты уравнения. Дополнительное сравнение невязок, рассчитанных с идентифицированными параметрами \(\hat\lambda\) и истинными параметрами \(\lambda_{\text{true}}\), приведено на Рис. 7. Оно позволяет оценить, насколько остаточная невязка при идентификации отличается от невязки, обусловленной исходной моделью и используемыми данными.
Значение невязок для некоторых уравнений превышает 5 %, что может быть связано с недостаточной точностью модели, шумом измерений и особенностями численного дифференцирования или нахожнем локального экстремума. Однако в целом результаты идентификации демонстрируют хорошее согласование с истинными параметрами и подтверждают эффективность предложенного подхода к оценке дисбаланса роторов.
Время идентификации параметров составило \(T_{\text{id}} = 243\) с для сигнала с частотой дискретизации \(f_s = 25\) кГц, количество точек \(N = 178\), количество опор \(N_s = 4\), число идентифицируемых параметров \(n_\lambda = 6 + 2 \cdot 4 \cdot 3 = 30\), на вычислительной системе, характеристики которой приведены в Табл. 2.
| Параметр | Значение |
|---|---|
| Процессор | Intel Core i9-12900K, 3.20 ГГц |
| Число ядер / потоков | 16 (8P + 8E) / 24 |
| Оперативная память | 64 ГБ |
| Операционная система | Windows, 64-разрядная (x64) |
| Язык программирования | Python 3.12 |
| Библиотека оптимизации | SciPy (least_squares, метод TRF) |
| Режим вычислений | один поток |
7 Заключение
В работе рассмотрена задача идентификации параметров дисбаланса соосных роторов вертолёта на основе измерений вибрационных перемещений. Разработан подход, основанный на математическом моделировании динамики механической системы и минимизации невязки уравнений движения.
Сформирована математическая модель, учитывающая поступательные и вращательные движения конструкции, упругие и диссипативные свойства опор, а также возмущающие воздействия, обусловленные дисбалансом верхнего и нижнего роторов. Для описания дисбаланса введены параметры массы, эксцентриситета и начальной фазы дисбалансирующих масс.
Рассмотрены вопросы восстановления кинематических характеристик системы по экспериментальным данным. При использовании измеренных перемещений скорости и ускорения могут быть получены численным дифференцированием. При использовании акселерометров требуется предварительное восстановление траекторий измерительных точек по ускорениям с учётом влияния смещения нуля, погрешностей интегрирования и начальных условий. Полученные кинематические характеристики используются для формирования невязки уравнений движения.
Проанализированы особенности идентифицируемости параметров. Показано, что центробежные силы непосредственно определяются статическим моментом дисбаланса и фазой, тогда как раздельное определение массы и эксцентриситета требует дополнительных условий или информативных составляющих модели. Данное обстоятельство необходимо учитывать при интерпретации результатов оптимизации и оценке достоверности восстановленных параметров.
Работоспособность разработанного подхода проверена на синтетических данных с известными параметрами. В рассмотренном примере относительная погрешность определения статического момента дисбаланса составила 5 % для верхнего ротора и 1 % для нижнего, а абсолютная ошибка определения начальной фазы не превысила 3°. При этом для верхнего ротора наблюдалось заметное отклонение оценок массы и эксцентриситета от истинных значений при существенно меньшей погрешности их произведения. Полученные результаты согласуются с рассмотренными особенностями идентифицируемости параметров.
Таким образом, результаты моделирования подтверждают возможность применения предложенного подхода для оценки параметров дисбаланса соосных роторов по данным о движении конструкции. Дальнейшее развитие метода связано с исследованием устойчивости идентификации при наличии шумов измерений, уточнением модели опор, а также ускорением идентификации за счёт использования более эффективных алгоритмов оптимизации и методов регуляризации.
Список литературы
Сведения об авторах
Максим Александрович Бодаков — Центральный институт авиамоторостроения им. П.И. Баранова, Россия, Москва. Тел.: (915) 165-6164, e-mail: mabodakov@ciam.ru
Михаил Сергеевич Куц — Московский государственный технический университет им. Н.Э. Баумана, кафедра «Основы конструирования машин» РК-3, Россия, Москва.