Устойчивость стационарных движений

Исследуем устойчивость стационарных движений диска относительно переменных θ, ψ', θ', φ' по уравнениям линейного приближения уравнений динамики (10).

Приведём уравнения (10) к безразмерному виду. Для этого поделим обе части второго уравнения (10) на g/ρ и положим единицу измерения времени, равной

t* = (ρ / g)^1/2.

Сохраним за безразмерными переменными прежние обозначения. Введём новые переменные ωψ, ωθ, ωφ согласно формулам

ωψ = ψ', ωθ = θ', ωφ = φ'

и преобразуем полученные из (10) уравнения в безразмерной форме к системе четырёх дифференциальных уравнений первого порядка относительно переменных θ, θ' = ωθ,

12 cos θ ωφ' + (7 + 5 cos 2θ ) ωψ' - 4 sin θ ωθ (ωφ + 5 cos θ ωψ ) = 0,

5 ωθ' + 5 cos θ sin θ ωψ² + 6 sin θ ωφ ωψ + 4 cos θ = 0, (13)

3 (ωφ' + cos θ ωψ') - 5 sin θ ωθ ωψ = 0.

Представим (13) в матричном виде

F(X_v,X) = 0, (14)

F - вектор-столбец 4×1 левых частей (13), X = (θ, ωψ, ωθ, ωφ)^T, X_v = (θ', ωψ', ωθ', ωφ')^T.

Примем за невозмущённое стационарное движение диска. Обозначим через X₀, X_v0 векторы X, X_v, отвечающие невозмущённому движению. Тогда, в соответствии с определением (11), будем иметь

X₀ = (θ₀, ωψ₀, ωθ₀, ωφ₀)^T, X_v0 = (0, 0, 0, 0)^T.

Линеаризовав (14) в окрестности невозмущённых значений переменных X₀, X_v0, получим

A_v X_v + A_q (X-X₀) = 0, (15)

A_v = dF/dX_v, A_q = dF/dX - матрицы 4×4, вычисленные при X_v = X_v0, X = X₀.

Характеристическое уравнение системы (15) имеет вид

||λ A_v + A_q || = 0. (16)

Далее представлена программа Mathematica 7 построения полинома в левой части (16). Предварительно следует исполнить программу составления уравнений движения диска.

Программа составления полинома в левой части (16)

(* Приведение уравнений Чаплыгина к безразмерной форме, определение замены переменных ωψ = ψ', ωθ = θ', ωφ = φ' и составление системы уравнений ( 13 ) *)

sub8 = { ψ' [ t ] -> ωψ [ t ] , θ'[ t ] -> ωθ [ t ] , φ'[ t ] -> ωφ [ t ] } ;

sub9 = Join [ sub8 , D [ sub8 , t ] ] ;

Eqch = Join [ { { θ' [ t ] - ωθ [ t ] } } , 2 /(3 m ρ ^2) EqChaplygin /. sub9 /. { g -> ρ ν ^2 } ] /. ν-> 1 // Simplify ;

(* Получение уравнения стационарных движений диска ( 12 ) *)

sub10 = { ωψ -> ( ωψ0 E^ ( 0 #1 ) & ) , θ ->( θ0 E ^ ( 0 #1 ) & ) , ωθ ->(0 #1 & ) , ωφ -> ( ωφ0 E ^ ( 0 #1 ) & ) } ;

Eqch0 = Eqch /. sub10 // Expand

(* Определение векторов X = (θ, ωψ, ωθ, ωφ) T , Xv = (θ', ψω', ωθ', ωφ' ) T , нахождение матрицы Av , Aq *)

X = { { θ [ t ] }, { ωψ [ t ] }, { ωθ [ t ] }, { ωφ [ t ] } };

Xv = D[ X,t] ;

Av = DV [ Eqch , Xv ] /. sub10 // Simplify ;

Aq = DV [ Eqch , X ] /. sub10 // Simplify ;

Row [ { "Av= " , Av // TraditionalForm } ]

Row [ { "Aq= " , Aq // TraditionalForm } ]

(* Составление полинома в левой части ( 16 ) *)

Poly = Det [ Av λ + Aq ] // Simplify

Результат выполнения программы

Список, содержащий левую часть уравнения стационарных движений

Матрица A_v

Матрица A_q

Характеристический полином

В обычной математической нотации уравнение (16) имеет вид

λ² (5 λ² + a) sin ² θ₀ = 0, (17)

где

а = 12 ωφ₀²+ 5 ωψ₀² + 14 ωφ₀ ωψ₀ cos θ₀ - 4 sin θ₀. (18)

Согласно теореме Ляпунова о неустойчивости неравенство

a ≥ 0 (19)

есть необходимое условие устойчивости стационарных движений диска. Достаточность условия (19) устойчивости (со строгим знаком неравенства) доказана в [8] путём построения функции Ляпунова.

Рассмотрим различные стационарные движения диска. Уравнение стационарных движений (12), записанное в безразмерной форме, имеет вид

cos θ₀ sin θ₀ ωψ₀²+ 6/5 sin θ₀ ωψ₀ ωφ₀ + 4/5 cos θ₀ = 0. (20)

Предположим сначала, что ωψ₀ ≠ 0. Выразив ωφ₀ из уравнения (20), получим

ωφ₀ = - сtg θ₀ (4 + 5 ωψ₀² sin θ₀) / (6 ωψ₀). (21)

Подставив (21) в (18), представим условие (19) в виде неравенства

(5 ωψ₀⁴ (2 - cos 2θ₀) + 12 ωψ₀² cos 2 θ₀ / sin θ₀ + 16 ctg² θ₀ ) / (3 ωψ₀ ²) ≥ 0. (22)

Условие (22) (с исключённым знаком равенства) выполняется при любых значениях ωψ₀ в случае

d = 8 (-21 - 20 cos 2θ₀ + 19 сos 4θ₀) сosec² θ₀ < 0, (23)

d - дискриминант квадратного относительно ωψ₀² трёхчлена в числителе левой части (22). Решение неравенства (23) имеет вид

0 < θ₀ < 1,246, 1,9 < θ₀ < π, (24)

В случае

1,246 ≤ θ₀ ≤ 1,9 (25)

дискриминант d ≥ 0 и неравенство (22) выполняется при

ωψ₀² ≤ (-6 cos 2θ₀ / sin θ₀ - 0,5 d ^1/2) / (10 -5 cos 2θ₀ ), ωψ₀² ≥ (-6 cos 2θ₀ / sin θ₀ +

+ 0,5 d ^1/2) / (10 - 5 cos 2θ₀ ). (26)

В размерных переменных условия (26) имеют вид

ωψ₀² ≤ (g / ρ) (-6 cos 2θ₀ / sin θ₀ - 0,5 d ^1/2) / (10 -5 cos 2θ₀ ), ωψ₀² ≥ (g / ρ) (-6 cos 2θ₀ / sin θ₀ + (27)

+ 0,5 d ^1/2) / (10 - 5 cos 2θ₀ ).

В случае верчения диска при θ₀ = π/2, ωφ₀ = 0 из (26) найдём

ωψ₀² ≥ 4/ 5

или в размерных переменных

ωψ₀² ≥ 4/ 5 g / ρ .

Рис. 6. Карта сечений поверхности (20) плоскостями ωφ₀ = const

Полученным результатам можно дать геометрическую интерпретацию . Уравнение (20) определяет в пространстве переменных ωφ₀, ωψ₀, θ₀, поверхность, которую характеризует на плоскости параметров θ_0,ωψ₀ карта [12] её сечений плоскостями

ωφ₀ = const,

показанная на рис. 6. Причём прямая

θ0 = π/2,

уравнение которой получается из (20) при ωφ₀ = 0, отвечает верчению диска, а точки других сечений - качению по окружностям (прецессиям).

На рис. 7 представлен, полученный в [13], график кривой

a( θ_0,ωψ₀) = 0, (28)

отмеченный красным цветом, разделяющий плоскость параметров θ₀,ωψ₀ на две области: внутри и вне "восьмёрки". При значениях параметров из области, закрашенной розовым цветом, стационарные движения диска неустойчивы, из незакрашенной области - устойчивы.

Рис. 7. Область неустойчивости (закрашена розовым цветом) стационарных движений

на плоскости параметров θ_0,ωψ₀

На рис. 8 карта сечений поверхности (20) и кривая (28) изображены вместе. Согласно условию (19) точкам внутри "восьмёрки" соответствуют неустойчивые стационарные движения, а точкам вне "восьмёрки" - устойчивые. Из рис. 8 видно, что в случае (24) стационарные движения диска устойчивы при любых значениях ωψ₀ , а в случае (25) - при выполнении условий (26) с исключёнными знаками равенств, что совпадает с известными результатами [14].

Рис. 8. Карта сечений поверхности (20) плоскостями ωφ₀ = const и область неустойчивости стационарных движений (закрашена розовым цветом)

Пусть теперь ωψ₀ = 0. Решению уравнения (20) при этом отвечает качение диска по прямой, в котором θ₀ = π/2. Из (18), (19) найдём необходимое (с добавлением знака равенства) и достаточное условие устойчивости этого стационарного движения в виде

ωφ₀² > 1/ 3

или в размерных переменных

ωφ₀² > g /(3 ρ).

Уравнение стационарных движений (20) имеет также решение

ωφ₀ = ωψ₀ = 0, θ₀ = π/2, (29)

соответствующее покою диска, при котором его плоскость вертикальна. В этом случае первое и третье уравнения движения диска (13) описывают плоское движение перевернутого маятника, точка подвеса которого совпадает с точкой контакта диска и неподвижной плоскости,

θ' = ωθ,

5 ωθ' + 4 cos θ = 0,

Вертикальное положение этого маятника, разумеется, неустойчиво.

Далее представлена программа, реализующая проведённое исследование устойчивости.

Программа исследования устойчивости стационарных движений

(* Нахождение коэфффициента a при λ ^ 2 в характеристическом полиноме *)

a = Coefficient [ Poly , λ , 2 ]

(* Выражение ωφ0 из уравнения стационарных движений и подстановка ωφ0 в a, af = a (θ0 , ωψ0 ) *)

sub11 = Solve [ Eqch0〚 3 , 1 〛== 0 , ωφ0 ] 〚 1 〛 // Simplify ;

af = a /. sub11 // Simplify

(* Нахождение дискриминанта d *)

d = Discriminant [ af 〚 4 〛 /. { ωψ0 ^ 4-> p ^ 2 , ωψ0 ^ 2-> p } , p ] // Simplify

(* Решение неравенства d > 0 *)

Reduce [ d > 0 && 0 < θ0 < Pi , { θ0 } , Reals ] // N

(* Решение уравнения a(θ0, ωψ0 ) = 0 относительно ωψ0 ^ 2 , p = ωψ0 ^ 2 *)

s = Solve [ af 〚 4 〛 == 0 /. { ωψ0 ^ 4 -> p ^ 2 , ωψ0 ^ 2-> p } , { p } , Reals ] // Simplify

(* Решение уравнения a(θ0, ωψ0 ) = 0 относительно ωψ0 ^ 2 в случае верчений, p = ωψ0 ^ 2 *)

s /. { θ0 -> Pi / 2 , p -> ωψ0 ^ 2 } // Simplify

(* Решение неравенства a > 0 относительно ωφ0 ^ 2 в случае качений, p = ωψ0 ^ 2 *)

Reduce [ a > 0 /. { θ0 -> Pi / 2 , ωψ0 -> 0 , ωφ0 ^ 2-> p } , p ] // Simplify

(*Решение неравенства a > 0 относительно ωφ0 ^ 2 в случае покоя диска в вертикальной плоскости *)

Reduce [ a > 0 /. { θ0-> Pi / 2 , ωφ0 -> 0 , ωψ0 -> 0 } , ωφ0 ] // Simplify

(* Построение графика кривой a (θ0, ωψ0 ) = 0 *)

p1 = Show [ RegionPlot [ { af < 0 } , { θ0 , 0 , Pi } , { ωψ0 , - 1.5 , 1.5 } , BoundaryStyle -> Directive [ Red ] , PlotStyle-> Lighter [ Pink , .6 ] , FrameLabel -> { θ0 , ωψ0 } ] ]

(* Построение карты сечений поверхности ( 20 ) плоскостями ωφ0 = const *)

p2 = Show [ Table [ ContourPlot [ Eqch0 〚 3 , 1 〛 == 0 , { θ0 , 0 , Pi } , { ωψ0 , - 1.5 , 1.5 } , ContourLabels -> ( Text [ Framed [ Style [ StringForm [ "ωφ0 = `` " , ωφ0 ] , FontSize -> 7 ] , FrameMargins -> 1 , RoundingRadius-> 8 , Background-> Lighter [ Blue , .9 ] ] , { #1 , #2 } ] & ) ] , { ωφ0 , { - 2 , - 1.2 , - .5 , 0 , .5 , 1.2 , 2 } } ] ]

(* Построение карты сечений поверхности ( 20 ) и кривой ( 27 ) на одном рисунке *)

Show [ p1 , p2 ]