Журнал вычислительной математики и математической физики, 2022, T. 62, № 3, стр. 521-528

1. ВВЕДЕНИЕ

В последние годы прошлого столетия в публикациях появилось обобщение классического метода Ньютона–Эйлера [1] на случай упругодеформируемых звеньев кинематической цепи многозвенных динамических систем. Эту новую модификацию метода Ньютона–Эйлера стало принято называть обобщенным методом Ньютона–Эйлера (см. [2], [3]). Вслед за этим вышла в свет статья [4], в которой на основе сформулированного обобщенного метода Ньютона–Эйлера и нововведенных элементов символьного исчисления был наведен алгоритм численного решения первой задачи динамики (Inverse Dynamic Problem) упругих манипуляторов.

Дальнейшее развитие обобщенного метода Ньютона–Эйлера привело к построению алгоритмов численного решения второй и смешанной задач динамики (Direct and Mixt Dynamic Problems) упругих манипуляторов, которые нашли отражение в публикациях [5], [7]. Кроме того, обобщенный метод Ньютона–Эйлера в контексте решения смешанной задачи динамики был распространен в публикации [6] на класс упругих древовидных динамических систем (Arborescent Dynamic Systems). Впоследствии было показано [8], что применение обобщенного метода Ньютона–Эйлера позволяет существенно сократить машинное время численного решения задач управления упругих манипуляторов по сравнению с методами, использующими процедуру обращения их матрицы масс.

В публикации [9] на примере исследования особенностей символьного построения первой и второй задач динамики манипуляторов с абсолютно жесткими звеньями были раскрыты методологическая и формальная стороны метода Ньютона–Эйлера, откуда стало неопровержимо явствовать существование новой элементарной расчетной единицы – суперэлемента кинематической цепи. Таким образом, в свете указанного формально-методологического осмысления проблематика оптимизации методов динамического анализа может вполне однозначно сводиться к выбору тех методов, которые используют в своей конституции наипростейшие виды суперэлементов. В то же самое время необходимо отметить, что метод Ньютона–Эйлера в задачах динамики упругих манипуляторов не утратил своей актуальности по настоящее время [10]–[12].

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

2. ОБОБЩЕННАЯ ДИНАМИЧЕСКАЯ МОДЕЛЬ УПРУГИХ МАНИПУЛЯТОРОВ

Рассматривается манипулятор произвольной кинематической цепи (см. фиг. 1), состоящий из n линейно деформируемых звеньев всевозможной геометрической конфигурации ${{C}_{1}},{{C}_{2}},...,{{C}_{n}}$, связанных между собой n вращательными или поступательными соединительными узлами (кинематическими парами), идентифицируемыми логической переменной [2], [3]:

Из теории деформируемого тела известно, что малые линейные и угловые деформации упругого звена ${{C}_{j}}$ определяются перемещениями некоторой точки ${{E}_{{0j}}} \in {{C}_{j}}$ из недеформированной ${{S}_{{0j}}}$ в деформированную конфигурацию ${{S}_{j}}$ (фиг. 1) и исчисляются в жестко связанной со звеном ${{C}_{j}}$ подвижной системе отсчета ${{Q}_{j}}$ на основе следующих разложений [2], [3]:

Координатные функции Рэлея–Ритца в (2.2) ${}^{j}{{\Phi }_{{dj,i}}}$ и ${}^{j}{{\Phi }_{{rj,i}}},\;i = 1,...,{{N}_{j}},$ $j \in \{ 1,...,n\} $, где ${{N}_{j}}$ – число интерполяционных вектор-полиномов упругого звена ${{C}_{j}}$, подчинены векторному равенству [2], [3]

и образуют в точке ${{A}_{{0j + 1}}},\;j \in \left\{ {1,...,n} \right\}$, обобщенную матрицу форм [2], [3]:

Динамическая модель упругих манипуляторов выражается системой нелинейных алгебродифференциальных уравнений, включающих в себя: систему нелинейных дифференциальных уравнений динамического равновесия упругого звена ${{C}_{j}},\;j \in \left\{ {1,...,n} \right\}$ (см. [2], [3]):

и систему векторных алгебраических уравнений, выражающих кинематические зависимости между скоростями (линейными и угловыми) и ускорениями (линейными и угловыми) в характерных точках ${{A}_{j}},\;j \in \left\{ {1,...,n} \right\},$ упругих звеньев ${{C}_{j}}$:

3. ФОРМУЛИРОВКА И ЧИСЛЕННОЕ РЕШЕНИЕ СМЕШАННОЙ ЗАДАЧИ ДИНАМИКИ УПРУГИХ МАНИПУЛЯТОРОВ ОБОБЩЕННЫМ МЕТОДОМ НЬЮТОНА–ЭЙЛЕРА

Разобьем множество звеньев кинематической цепи (фиг. 1) $\{ M\} = \{ {{C}_{1}},{{C}_{2}},...,{{C}_{n}}\} $ на два подмножества $\{ {{M}_{1}}\} = \{ C_{1}^{'},C_{2}^{'},...,C_{n}^{'}\} $ и $\{ {{M}_{2}}\} = \{ C_{1}^{{''}},C_{2}^{{''}},...,C_{n}^{{''}}\} $, где $l + m = n$ по следующему признаку:

Тогда смешанная задача динамики упругих манипуляторов сформулируется следующим образом. Предполагается, что заданы законы изменения приводных усилий ${{\Gamma }_{i}},\;i = 1,...,l,$ и обобщенных ускорений $\ddot {q}_{k}^{r}{\text{,}}\;k = 1,...,m,$ принадлежащих соответственно подмножествам $\left\{ {{{M}_{1}}} \right\}$ и $\left\{ {{{M}_{2}}} \right\}$. Требуется определить законы изменения приводных усилий ${{\Gamma }_{k}},$ $k = 1,...,m,$ и обобщенных ускорений $\ddot {q}_{i}^{r}{\text{,}}\;i = 1,...,l,$ принадлежащих соответственно подмножествам $\left\{ {{{M}_{2}}} \right\}$ и $\left\{ {{{M}_{1}}} \right\}$.

Предлагаемый алгоритм решения смешанной задачи динамики упругих манипуляторов без обращения их матрицы масс составлен с использованием элементов итерационного символьного исчисления, нашедшего успешное применение в анализе упругих манипуляторов [4]. Он включает в себя три последовательно производимых для каждого момента времени $t \in \left[ {0,\;T} \right]$ рекурсивных цикла (см. [5]).

1-й рекурсивный цикл при j = 1, …, n требует вычисления следующих матриц: ${}^{j}{{R}_{{j - 1}}},\;{}^{j}{{T}_{{j - 1}}},\;{}^{j}{{J}_{{rrj}}},$ ${}^{j}{{J}_{{rej}}},\;{}^{j}{{J}_{{erj}}},\;{}^{j}{{J}_{{eej}}},$ векторов ${}^{j}C_{j}^{r},\;{}^{j}C_{j}^{e},\;{}^{j}{{h}_{j}},\;{}^{j}{{a}_{j}}$ и скаляра ${{\tau }_{j}}\;\;{\text{при}}\;\;{{\chi }_{j}} = 0$, что позволяет вслед за этим осуществить кинематический расчет цепи:

подчиняя (3.2) инициализациям: ${}^{0}{{V}_{0}} = 0$ и $\dot {q}_{0}^{e} = 0$.

2-й рекурсивный цикл при j = n, …, 1 после проведения инициализации:

предусматривает определение вспомогательных величин ${}^{j}J_{{rrj}}^{*},\;{}^{j}J_{{rej}}^{*},\;{}^{j}J_{{erj}}^{*},\;{}^{j}J_{{eej}}^{*},$ ${}^{j}C_{j}^{{r*}}$ и ${}^{j}C_{j}^{{e*}},$ $j = \left( {n - 1} \right),...,1$, по следующему алгоритмическому циклу.

Шаг 1. ${{K}_{j}} = {}^{j}J_{{rrj}}^{*} - {}^{j}J_{{rej}}^{*} \cdot {}^{j}J_{{eej}}^{{* - 1}} \cdot {}^{j}J_{{erj}}^{*};$ ${}^{j}{{m}_{j}} = - {}^{j}J_{{rej}}^{*} \cdot {}^{j}J_{{eej}}^{{* - 1}} \cdot {}^{j}C_{j}^{{e*}} + {}^{j}C_{j}^{{r*}}.$

Шаг 2. Если ${{\chi }_{j}} = 0$ или $j \ne 1$, то

$\{ {{G}_{j}} = {{K}_{j}} \cdot {}^{j}T_{{j - 1}}^{{}};\;\;{{H}_{j}} = {{K}_{j}} \cdot {}^{j}{{R}_{{j - 1}}} \cdot {}^{{j - 1}}{{\Phi }_{{j - 1}}};\;\;{}^{j}{{p}_{j}} = {{K}_{j}} \cdot {}^{j}{{h}_{j}} + {}^{j}{{m}_{j}}\} .$

Шаг 3. Если ${{\chi }_{j}} = 0,$ то ${{Z}_{j}} = \frac{1}{{{}^{j}a_{j}^{T} \cdot {{K}_{j}} \cdot {}^{j}{{a}_{j}}}}$.

Шаг 4. Если $j \ne 1$ и ${{\chi }_{j}} = 0,$ то $\{ {{U}_{j}} = {{Z}_{j}} \cdot {{K}_{j}} \cdot {}^{j}{{a}_{j}} \cdot {}^{j}a_{j}^{T} \cdot {{K}_{j}};\;\;{{\tilde {G}}_{j}} = {{G}_{j}} - {{U}_{j}} \cdot {}^{j}{{T}_{{j - 1}}};$

${{\tilde {H}}_{j}} = {{H}_{j}} - {{U}_{j}} \cdot {}^{j}{{R}_{{j - 1}}} \cdot {}^{{j - 1}}{{\Phi }_{{j - 1}}};\;\;{}^{j}{{\tilde {p}}_{j}} = {}^{j}{{p}_{j}} + {{K}_{j}} \cdot {}^{j}{{a}_{j}} \cdot ({{\tau }_{j}} - {}^{j}a_{j}^{T} \cdot {}^{j}{{p}_{j}}) \cdot {{Z}_{j}}\} .$

Шаг 5. Если $j \ne 1$ и ${{\chi }_{j}} = 1,$ то $\{ {{\tilde {G}}_{j}} = {{G}_{j}};\;\;{{\tilde {H}}_{j}} = {{H}_{j}};\;\;{}^{j}{{\tilde {p}}_{j}} = {}^{j}{{p}_{j}} + {{K}_{j}} \cdot {}^{j}{{a}_{j}} \cdot \ddot {q}_{j}^{r}\} .$

Шаг 6. Вычисление:

3-й рекурсивный цикл при j = 1,…,n использует результаты вычислений, произведенных в процессе предшествующих расчетных циклов, в целях определения для подмножества $\left\{ {{{M}_{1}}} \right\}$ обобщенных ускорений:

затем для множества $\left\{ M \right\}$ – ускорений звеньев, ускорений варьируемых параметров и, если требуется, внутренних усилий в парах: с учетом в (3.3) и (3.4) инициализаций ${}^{0}{{V}_{0}} = {{\left[ {0,\;0,\;g,\;0,\;0,\;0} \right]}^{T}}$ и $\ddot {q}_{0}^{e} = 0$ и, в заключение, для подмножества $\left\{ {{{M}_{2}}} \right\}$ – приводных усилий: где ${{F}_{{sj}}}$ и ${{F}_{{vj}}}$ – коэффициенты сухого и вязкого трений в парах.

Вследствие изложенного расчетного алгоритма численного решения смешанной задачи динамики упругих манипуляторов, определяемого формулами (3.1)–(3.7), представляется возможным дать схематическую иллюстрацию формирования (см. фиг. 2а) и функционирования (фиг. 2б) оптимального суперэлемента кинематической цепи, соответствующего методу Ньютона–Эйлера [9], [13].

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

Метод численного интегрирования нелинейных уравнений движения базируется на общей схеме численного интегрирования Ньюмарка:

где ${{x}_{j}} = {{[q_{j}^{r},{{(q_{j}^{e})}^{T}}]}^{T}}$; $\beta $ и $\gamma $ – постоянные ($0 \leqslant \beta \leqslant 0.5\;\;{\text{и}}\;\;0 \leqslant \gamma \leqslant 1$). В случае когда $\beta = \gamma = 0$, система уравнений (4.1) вырождается в явную схему интегрирования Ньюмарка, использование которой зачастую сопряжено с необходимостью выбора чрезвычайно малого шага интегрирования $\Delta t$ (см. [5], [11]). В случае же $\beta \ne 0$ или $\gamma \ne 0$ речь идет о неявной схеме численного интегрирования Ньюмарка (14) системы нелинейных дифференциальных уравнений.

Решение системы линеаризованных уравнений движения упругих манипуляторов осуществляется посредством метода Ньютона [5]. Точное вычисление якобиана (матрицы итераций) в расчетной формуле метода Ньютона, как известно, ускоряет процесс сходимости, но предполагает громоздкий расчет тангенциальных матриц жесткости и демпфирования, особенно при большом общем числе степеней свободы рассматриваемых систем. В то же время ассимиляция якобиана с матрицей масс манипулятора, будучи оправданной для жестких систем, вызывает серьезные сомнения при наличии в кинематической цепи упругих звеньев [13]. В этой связи представляет особый практический интерес компромиссный подход к преодолению указанных трудностей вычислительного характера, заключающийся в замещении точного выражения матрицы итерации ее редуцированным аналогом [5] (фиг. 3).

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

Рассмотрим численное решение смешанной задачи динамики упругих манипуляторов обобщенным методом Ньютона–Эйлера на примере изображенного на фиг. 4 пространственного пятизвенного манипулятора. Звенья ${{C}_{1}}$ и ${{C}_{4}}$ манипулятора полагаются жесткими, а звенья ${{C}_{2}}$ и ${{C}_{3}}$ – линейно-упругими, а опорная стойка ${{C}_{0}}$ – неподвижной и недеформируемой. В точке, совпадающей с центром тяжести звена ${{C}_{4}}$, устанавливается груз массой $m{\kern 1pt} * = 200\;кг$, который подлежит транспортировке в процессе движения манипулятора [14].

Длины звеньев заданного манипулятора (фиг. 4) полагаются равными: ${{L}_{1}} = 0.3\;{\text{м}},$ ${{L}_{2}} = 4.0\;{\text{м}},$ ${{L}_{3}} = 5.0\;{\text{м}},$ ${{L}_{4}} = 0.5\;{\text{м}}$ при плотности $\rho = {{\rho }_{i}} = 2699\;{\text{кг/}}{{{\text{м}}}^{3}},$ $i = 1,...,4$. Жесткие звенья ${{C}_{1}}$ и ${{C}_{4}}$ принимаются сплошными, круглого поперечного сечения радиусов ${{r}_{1}} = {{r}_{4}} = 0.05\;{\text{м}}$. Упругие звенья ${{C}_{2}}$ и ${{C}_{3}}$ предполагаются полыми, кольцевого поперечного сечения с радиусами внешних и внутренних окружностей, соответственно равными $r_{{}}^{e} = r_{2}^{e} = r_{3}^{e} = 0.05\;{\text{м}}$ и $r_{{}}^{i} = r_{2}^{i} = r_{3}^{i} = 0.04\;{\text{м}}$, а величины их модулей Юнга будут полагаться равными $E = {{E}_{2}} = {{E}_{3}} = 6.895 \times {{10}^{{10}}}\;{\text{Н/}}{{{\text{м}}}^{2}}$.

Обратная задача динамики изображенного на фиг. 4 пространственного упругого манипулятора формулируется при следующих законах движения звеньев:

где пределы конечных промежутков ${{T}_{s}} = 15\;{\text{с}}$ и $T = 25\;{\text{с}}$.

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

Задача предусматривает определение образующегося вследствие упругих деформаций отклонения точки Р исполнительного органа манипулятора (см. фиг. 4) от фиктивного положения, занимаемого той же точкой в предположении абсолютной жесткости всех звеньев. Это отклонение в проекциях на оси неподвижной системы координат XYZ выражается через упругие перемещения звеньев ${{C}_{2}}$ и ${{C}_{3}}$ следующими геометрическими соотношениями:

Численная реализация динамического анализа исследуемого манипулятора производится посредством метода “предиктора-корректора” при поддержке неявной схемы численного интегрирования Ньюмарка с аппроксимацией в форме редуцированной матрицы итераций $J = [M]$ (см. [5]). Результаты вычисления отклонения точки Р в проекциях на оси координат X, Y и Z при шаге интегрирования $\Delta t = 0.1\;{\text{с}}$ представлены на фиг. 5. Сопоставление итогов численного решения задачи с соответствующими показателями в публикации [14], в которой используется процедура обращения матрицы масс, убеждаемся в высокой точности моделирования, осуществленного с минимальным расчетным временем ${{T}_{c}} = 6\;с$, тогда как зарегистрированное в [14] хронометрическое показание (на машине того же поколения) ${{T}_{c}} = 127\;с$.

4. ЗАКЛЮЧЕНИЕ

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