Численное решение задачи динамики с обходом сингулярностей Текст научной статьи по специальности «Математика»

УДК 004.77

И. Е. Каспирович, Р. Г. Мухарлямов

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ДИНАМИКИ С ОБХОДОМ СИНГУЛЯРНОСТЕЙ

Ключевые слова: параметры Родрига-Гамильтона, численное интегрирование, стабилизация связей, множители Лагранжа.

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

Keywords: Rodrigues-Hamilton parameters, numerical integration, constraint stabilization, Lagrange multipliers.

During numerical integration of rigid body dynamics equations singularities caused by a structure of these equations can appear. Some simplest numerical schemes of singularity skips are considered in the case of sphere rolling on rotating plane. The method of Lagrange multipliers is applied for accounting of nonholonomic constraints caused by the condition of rotating without slipping. Solution of Euler equations for rigid bodies can be obtained analytically. However, singularities appear during numerical integration of angular variables. When the angle of intrinsic rotation coincides with precession the system of motion equations is getting degenerated that causes these singularities. Some local modifications of algorithm that realizes the process of numerical integration are considered. Singularity points are avoided in these modifications. Deviations caused by such skips can be minimized using the method of constraint stabilization.

Введение

Рассматривается задача численного решения уравнений относительно угловых переменных в задаче о качении шара на вращающейся плоскости [1]. Само качение, с физической точки зрения, представляет собой один из наиболее простых примеров задач о движении твердого тела. Изменение угловых переменных шара от времени описывается нелинейной автономной системой дифференциальных уравнений. Если выбрать в качестве угловых переменных компоненты нормированного кватерниона [2], то система уравнений движения шара имеет сингулярность при совпадении углов прецессии и собственного вращения. Предлагаемый метод устранения данной проблемы обусловлен искусственными пропусками точек, в которых имеются сингулярности. Стоит отметить, что данный метод предполагает значительные отклонения решения, полученного численным интегрированием, от аналитического. Чем больше было совершено пропусков при численном интегрировании, тем внушительнее отклонение. Минимизировать такие отклонения возможно методом стабилизации связей [3-5].

Параметры Родрига-Гамильтона

Лагранжиан шара единичной массы, катящегося без проскальзывания на вращающейся плоскости вне поля тяжести, имеет следующий вид:

1 = 1-(х2+у2) + 1-(ш2+ш2 + ш I), (1) где х, у - координаты центра масс (плотность шара принимается за константу), J - момент инерции относительно центральной оси, х = а ш = (ш1, ш2, ш3)Т - вектор угловой скорости шара.

Условие качения без проскальзывания накладывает на движение шара две неголономные связи:

/1=х-аш2+ау=0, /2=у+аы1-Пх=0, ( )

где а- радиус шара, а П -угловая скорость вращения плоскости.

Угловые скорости шара выражаются как квазискорости посредством углов Эйлера <р, ф, в равенствами:

= СОБфв + БШ вБ\пфф, Ш2 = Бтфв — БШбсО Бфф, Ш3 = ф + СОБвф.

Лагранжиан (1) в углах Эйлера принимает вид: I = 1 (х2 + у2) + ^ (в2 + ф2 + ф2 + 2ффсОБв).

Лагранжиан содержит угловые переменные, что позволяет получить уравнения движения для заданных угловых переменных. Однако, в случае Эйлера (все компоненты диагонального тензора инерции равны) возникает вырождение системы уравнений движения.

Система уравнений движения в углах Эйлера записывается в виде:

в + ффБтв = 0, ф + фсоБв — фвБтв = 0, ф + фсоБв — фвБтв = 0, где при в = п второе и третье уравнения становятся эквивалентными. Следовательно, разрешимость вторых производных становится невозможной. Для устранения этой неразрешимости можно ввести параметры Родрига - Гамильтона. При этом составляющие угловой скорости выражаются через компоненты нормированного кватерниона Л = {А0Д} следующим образом:

^1—2Я0Я1—2Я1ЯО+2ЯЗЯ2—2Я2ЯЗ, ^2=2ЯОЯ2—2Я2ЯО+2Я1ЯЗ—2ЯЗЯ1, ^З=2Я0ЯЗ—2ЯЗЯО+2Я2Я1—2Я1Я2.

Условие нормировки кватерниона Я0 2 + л2 = 1 позволяет перейти от параметров Родрига - Гамильтона к угловым переменным а, р, у

Я0 = СО5аСО5^СО5у, Я1=со5асо5^51пу, А2=со5а5тр, Л3=5т а.

Угловые скорости ^ выражаются через переменимте а, р, у как линейные формы:

= кУ + тФ + п№> I = 1,2,3,

(3)

где 11,т.1,щ - функции угловых переменных а,р,у. Лагранжиан (1) в переменных а, р, у запишется

Ь = ±(х2 + у2) + 2](а2 + соБ2а$2 + ^а^/Зу2). (4)

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

Для решения системы уравнений движения, следующих из выражения(4) с неголономными связями (2) применяется метод произвольных множителей Лагранжа и1,и2

X = и1, У = и2,

а = —^Бт2а($2 + соб2РУ2) + ~~(и2п1 — ЩЩ),

у = 2tgaay + 2Ъ%13у13 +

(5)

4/С052ЯС052Д

(и2т1 — и1т2), (и211—и112).

При применении методов численного интегрирования к данной системе возникают сингулярности при а = пп^2'Р = ПП^2 ,п в правой части выраже-нийдляД и у. Однако, данный тип сингулярности не приводит к вырождению системы уравнений.

Устранение сингулярностей

Для устранения сингулярностей при численном интегрировании методом Эйлера системы уравнений (5) рассматриваются два численных метода. Первый заключается в пропуске точки сингулярности при приближении значений координат а или р к ^у^ Второй метод предполагает разбиение численного решения на участки. При этом задача Коши для нового участка состоит в установлении значений всех координат в начальный момент времени нового участка в виде тех же значений для конечного момента времени предыдущего участка за исключением а или р, которые по аналогии пропускают значение п12 Та-

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

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

(6)

1/1 = КЪ + к2[2, 1/2 = ЙзЛ +

где £ I - произвольные константы, называемые параметрами возмущения.

Для асимптотической устойчивости решения системы (6) вещественные части собственных значений 1 = 1,2 матрицы коэффициентов возмущения К должны быть строго отрицательны Re(Яг) < 0, что приводит к следующему неравенству:

Re(trK ± ЛхЧГ—4НёК) < 0,

(7)

Если заданы значения параметров к2 и к3 одного знака, то неравенство (7) значительно упрощается, и для к1и к4 можно получить следующую систему неравенств:

>

кг + к4 < —2^к2к3.

(8)

Тогда множители Лагранжа из уравнений (3), (5) выражаются с помощью уравнений возмущений связей (2), (6) как функции обобщенных координат, скоростей и параметров возмущения к^.

Метод пропусков

Рассматривается пример обхода одной точки сингулярности для системы (5). Правые части этой системы обращаются в бесконечность при р = ^у^ Систему уравнений (5) можно переписать в следующем виде:

"а = 9а.

^ = д^ "у = дг + ^Ру^

а = va' У = ^у

(9)

где функции да, др и ду + 2tg^yДпредставляют собой правые части уравнений (5) с определенными множителями Лагранжа и параметрами возмущений, удовлетворяющими неравенствам (8).

Численный алгоритм, реализующий метод Эйлера, для уравнений (9) запишется следующим образом:

(¡+1) (О I а) Vа ="а + Т9а(1),

1>(!+1) = и® + тп (0

(¡ + 1) (0 , г

$ + т(ду

(1) +

а)

а(1+1) = а(о +

рц+ъ = р(я + Ц0,

у(1+1) = у(С) + т„(0,

где т - шаг интегрирования.

Пусть на д + 1-ом шаге угол принадлежит

интервалу (^/2 - £,п/2 + £), тогда для того, чтобы избежать сингулярности, предлагается заменить на Д(м+1) = + пропустив тем самым данный интервал. Параметры^ и рзадаются вручную (рис.1).

ключением бугорка, соответствующего интервалу,

V2"

который содержит значение Д = п/~

Рис. 1 - Схема обхода

Численный алгоритм, реализующий данную схему, выглядит следующим образом:

if е(п/2-£,п/2 + s) then =

= +

Очевидно, что S > 28, чтобы на д + 2-ом шаге угловой параметр снова не оказался внутри интервала (^/2 — £, V2 + Тогда для уравнения, имеющего слагаемое, которое вызывает сингулярность, разностная схема изменит свой первоначальный вид на следующий:

„(¿+2) _ Гм+1) , Т(2 tgft(^ + tg£ ы+1) ы+1) + VY _VY +T{2l — tgp(v)tg8V(S VV +дУ

Результат применения изложенного метода к интегрированию системы (5) представлен на рисунке 2.

Задавались следующие значения параметров возмущения k2 _ — 1, к3 _ — 1, fc4 _ —1, кг менялся от -12.5 до -50. Параметры пропуска выбирались соответственно 8 _ 0.87, р _ 0.05. Изменение величины р в момент пропуска t = 2.95с представляет собой ступеньку. При увеличении модуля параметра кг график ведет себя более плавно. Однако, существует граница для допустимого значения модуля этого параметра, связанная с оценками отклонения от уравнений связей.

Рис. 2 - Зависимость р от времени

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

Рис. 3 - Траектория центра масс

Метод деления интервала

При данном подходе при вхождении угла в

тот же интервал, цикл, реализующий алгоритм, останавливается. Далее запускается новый цикл, в котором начальные значения всех координат принимаются за значения на (д + 1)-ом шаге, за исключением д(м+1) =п^2 + 8. Цикл для временной переменной при этом проходит непрерывно. Далее решение при всех циклах склеивается, и получается решение на всем интервале интегрирования.

Результат применения изложенного метода к интегрированию системы (5) представлен на рисунке 4. Видно, что характерная ступень для графика 2 на рисунке 4 почти отсутствует. Однако, происходит укру-чение линии после склеивания циклов. Это же картина характерна и для траектории центра масс.

Рис. 4 - Зависимость^ от времени

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

Заключение

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

Работа выполнена при поддержке РФФИ (проект 16-08-00558)

Литература

1. Ю. И. Неймарк, Н. А. Фуфаев, Динамика неголономных систем, Наука, Москва,1967г, 498с.

2. В. Н. Кошляков, УМЖ, 40, 2, 154-164 (1988)

3. Р. Г. Мухарлямов, Вестник Казан. технол. ун-та, 16, 12,278--285 (2013)

4. R.G. Mukharlyamov, C. T. Deressa, Вестник КНИТУ, 17, 11, 236-243 (2014)

5. Н.В. Абрамов, А.А. Ахметов, Р. Г. Мухарлямов, Вестник Казан. технол. ун-та, 17,20, 319-325 (2014)

© И. Е. Каспирович - студент-магистр, кафедра теоретической физики и механики, Российский университет дружбы народов, kaspirovich.ivan@mail.ru; Р. Г. Мухарлямов - профессор, доктор физико-математических и естественных наук, кафедра теоретической физики и механики, Российский университет дружбы народов, robgar@mail.ru.

© 1 E. Kaspirovich - B.Sc., Department of Theoretical Physics and Mechanics. Peoples' Friendship University of Russia, kaspiro-vich.ivan@mail.ru; R. G. Mukharlyamov - Doctor of Sc., Full Professor, Department of Theoretical Physics and Mechanics. Peoples' Friendship University of Russia, robgar@mail.ru.