CC BY
4Нелинейная динамика. 2016. Т. 12. № 1. С. 99-120. Полнотекстовая версия в свободном доступе http://nd.ics.org.ru Б01: 10.20537/па1601007
ОРИГИНАЛЬНЫЕ СТАТЬИ
УДК: 532.3
М8С 2010: 70Е15, 65Рхх
Хаотическая динамика в задаче о падении тела винтовой формы в жидкости
В данной работе изучается процесс свободного падения трехлопастного винта в жидкости. Исследование проводится в постановках идеальной и вязкой жидкостей. Для случая идеальной жидкости исследована устойчивость равноускоренных вращений (решений Стеклова). При исследовании движения в постановке вязкой жидкости была построена феноменологическая модель вязких сил и моментов. Построена карта показателей Ляпунова и бифуркационные деревья. Показано, что в зависимости от параметров системы возможны два типа движения: квазипериодическое и хаотическое, переход к хаосу происходит через каскад бифуркаций удвоения периода.
Ключевые слова: идеальная жидкость, вязкая жидкость, движение твердого тела, динамическая система, устойчивость движения, бифуркации, карта показателей Ляпунова
Получено 10 февраля 2016 года После доработки 04 марта 2016 года
Работа Е. В. Ветчанина выполнена в рамках базовой части государственного задания вузам, а также поддержана грантами РФФИ № 15-08-09093-а и №15-38-20879 мол_а_вед. Работа В. А. Тененева выполнена при поддержке гранта РФФИ № 14-01-00395-а.
Тененев Валентин Алексеевич tenenev@istu.ru Илалетдинов Ленар Фаритович ano-lenar@yandex.ru
Ижевский государственный технический университет 426069, Россия, г. Ижевск, ул. Студенческая, д. 7
Ветчанин Евгений Владимирович eugene186@mail.ru
Удмуртский государственный университет 426034, Россия, г. Ижевск, ул. Университетская, д. 1 Ижевский государственный технический университет 426069, Россия, г. Ижевск, ул. Студенческая, д. 7
В. А. Тененев, Е. В. Ветчанин, Л. Ф. Илалетдинов
1. Введение
Задача о движении твердого тела в идеальной жидкости относится к классическим задачам гидродинамики и изучалась в различных постановках. В рамках модели идеальной жидкости движение тела по инерции описывается уравнениями Кирхгофа [1] и изучалось, например, Стекловым [2] и Ламбом [3]. Систематические исследования движения твердых тел в поле тяжести были проведены Стекловым [4] и Чаплыгиным [5].
В современной работе [6] были исследованы асимптотическое поведение и различные динамические эффекты, возникающие при падении различных тел в идеальной жидкости. Продемонстрировано, что в плоскопараллельном случае падение тела широкой стороной устойчиво, при этом амплитуда малых колебаний относительно устойчивого режима убывает как , а их частота растет пропорционально времени Ь. Для тел с винтовой симметрией исследована линейная устойчивость решений, соответствующих равноускоренным падениям с вращением, когда одна из осей тела остается вертикальной. Кроме того, показано, что при наличии начального толчка возможны режимы движения, типичные для диссипативных систем.
Следует иметь в виду, что равноускоренные движения, получаемые в рамках модели идеальной жидкости, практически не наблюдаются в действительности. Безусловно, чтобы смоделировать различные динамические эффекты, обусловленные как влиянием вязкости, так и влиянием формы тела, необходимо совместное численное решение уравнений Навье-Стокса и уравнения движения тела [7-11].
Однако расчет сопротивления жидкости непосредственно из уравнений Навье - Стокса весьма трудоемок, а жесткие требования к вычислительным ресурсам ограничивают возможности параметрических исследований. В связи с этим наиболее разумным представляется применение феноменологических моделей вязкого сопротивления, при этом движение тела описывается небольшим количеством обыкновенных дифференциальных уравнений. Такой подход открывает возможность использования для анализа динамики инструментария современной теории динамических систем.
Существенные успехи в исследовании динамических систем с сопротивлением, задаваемым феноменологически, были достигнуты при изучении задачи Максвелла о падающей пластинке [12]. Для исследования динамики падающей пластинки было построено множество моделей (например, [13-16]), однако все они приводят к различным результатам. Подробный анализ этих моделей с построением карт режимов, фазовых портретов аттракторов и бифуркационных диаграмм был проведен в недавней работе [17]. В настоящей работе мы также применяем методику, использованную в [17].
В данной работе мы рассматриваем падение винтового тела в жидкости и обсуждаем результаты, полученные в рамках моделей идеальной и вязкой (с феноменологическим законом сопротивления) жидкостей. Трехмерный характер обтекания винта создает существенные трудности при выборе вида закона сопротивления. В работах [18-21] для описания вращений тела с неподвижной точкой использовался линейный по угловым скоростям закон сопротивления, в [22] аналогичная задача рассматривается также при квадратичном законе сопротивления. Для определения закона сопротивления в широком диапазоне скоростей обычно используют квадратичное сопротивление [23]
где р — плотность жидкости, V — скорость движения тела относительно жидкости, 5 — площадь миделева сечения, Си — коэффициент сопротивления, в общем случае зависящий от числа Рейнольдса Кв.
Для определения вязкого сопротивления мы на основе численного моделирования стационарного обтекания винта строим феноменологическую модель. Результаты расчетов показали, что гидродинамическое сопротивление рассматриваемого тела хорошо аппроксимируется квадратичным по поступательным и угловым скоростям законом сопротивления с постоянными коэффициентами.
Отметим, что тела с винтовой симметрией могут иметь применение в мобильной робототехнике, поскольку такая форма дает преимущества при самопродвижении тела [24].
В разделе 2 приведены геометрические, массовые и инерционные характеристики исследуемого тела, а также уравнения движения в идеальной и вязкой жидкостях. В разделе 3 приведены результаты исследования в рамках модели идеальной жидкости. Проведен анализ устойчивости частных решений уравнений движения при движении из состояния покоя. С помощью численных экспериментов показано, что найденные частные решения соответствуют простым аттракторам, а какие-либо иные режимы движения не реализуются. В разделе 4 исследуется движение винта в вязкой жидкости. Указаны частные решения, соответствующие движениям с постоянными поступательной и угловой скоростями. С помощью численных расчетов показано, что частные решения неустойчивы. Построена карта показателей Ляпунова, продемонстрировано, что в зависимости от параметров движение системы соответствует либо простому, либо хаотическому аттрактору. Показано, что переход к хаосу происходит через каскад бифуркаций удвоения периода.
2. Математическая модель
2.1. Геометрия тела
Рассмотрим падение в безграничном объеме жидкости однородного тяжелого винтового тела, состоящего из центрального шара и трех лопастей в форме сплюснутых эллипсоидов вращения (см. рис. 1). Центр масс винта совпадает с центром шара, ось винтовой симметрии совпадает с диаметром шара.
Винт обладает следующими характеристиками:
— радиус центрального шара К = 0.1 м,
— радиус лопасти а = 0.2 м,
— толщина лопасти й = 0.04 м,
— расстояние от центра сферы до центра лопасти Н = 0.15 м,
— объем V = 1.238575 • 10"3 м3.
Ориентацию лопастей винта относительно центрального шара будем задавать углом Ф, причем Ф = 0° соответствует максимальному перекрытию плоскости 0\е\е2 (см. рис. 1), Ф = 90° — минимальному.
Центральный тензор инерции 3, тензор Л1 присоединенных масс, тензор Л2 присоединенных моментов инерции, тензор В, связанный с винтовой симметрией, зависят от величины угла Ф. В таблицах 1, 2 приведены значения этих тензоров для Ф = 20°, 30°, 45°, 60°.
Рис. 1. Геометрия тела.
Таблица 1. Центральный тензор инерции
ф J•Ю4/рь Ф J•104/рь
20° diag(2.12001, 2.12001, 4.16675) 30° diag(2.17355, 2.17355, 3.784)
45° diag(2.253855, 2.253855, 3.56706) 60° diag(2.32331, 2.32331, 3.33956)
Здесь рь — плотность материала винта.
Таблица 2. Присоединенные массы
Ф Ai•103/Pf Л2•103/Pf В•103/pf
20° diag(5.490, 5.490, 54.820) diag(0.841, 0.841, 0.191) diag(1.311, 1.311, -2.771)
30° diag(9.401, 9.40145.606) diag(0.734, 0.734, 0.399) diag( 1.825, 1.825, -3.702)
45° diag(16.917, 16.917, 30.559) diag(0.542, 0.542, 0.781) diag(2.168, 2.168, -4.196)
60° diag(24.727, 24.727, 16.395) diag(0.348, 0.348, 1.148) diag( 1.912, 1.912, -3.553)
Здесь pf — плотность жидкости.
2.2. Выбор систем координат. Кинематические соотношения
Для описания движения введем три системы координат: неподвижную Оахух, подвижную Охух, полученную из системы Оахух параллельным переносом, и подвижную Ов\в2ез, жестко связанную с телом (см. рис. 1). Будем считать, что
1. начала О подвижных систем координат совпадают с центром масс винта,
2. ось Оез подвижной системы координат совпадает с осью винтовой симметрии тела.
Положение начала подвижной системы координат относительно неподвижной будем задавать вектором г = (х,у,х)Т. Для описания ориентации подвижной системы координат введем единичные векторы а, в и 7, направленные вдоль осей Ох, Оу, Ох. Их
компоненты в подвижной системе координат образуют ортогональную матрицу перехода
( \
а1 71
Q
а2 @2 Y2 \аз вз Y3 У
е2
Таким образом, будут справедливы следующие кинематические соотношения
а = а х П, в = в х П, 7 = 7 х П, г = V, (2.1)
Рис. 2.
где V — поступательная скорость винта, П — угловая скорость винта.
Начальную ориентацию тела будем задавать следующими векторами:
а(0) = (cos е, 0, sin e)T, 0(0) = (0,1,0)T, 7(0) = (- sin e, 0, cos e)T,
где e — угол тангажа, соответствующий повороту винта вокруг оси Oy (см. рис. 2). Мы не рассматриваем повороты вокруг осей Ox и Oz в начальный момент времени, так они не влияют на динамику.
Для анализа пространственного движения введем углы Эйлера стандартным образом [25, 26]. Компоненты матрицы Р будут связаны с углами Эйлера следующими соотношениями:
^ cos ф cos ф — sin ф cos в sin ф sin ф cos ф + cos^ cos в sin ф sin в sin ф
Q
(2.2)
— cos ф sin ф — sin ф cos в cos ф — sin ф sin ф + cos ф cos в cos ф sin в cos ф sin ф sin в — cos ф sin в cos в
где ф — угол прецессии, в — угол нутации, ф — угол собственного вращения.
2.3. Уравнения движения тела в идеальной жидкости
Движение произвольного тела в идеальной жидкости описывается уравнениями Кирхгофа [1, 3]
p = p X П — UY,
(2.3)
M = M X П + p X V,
где p = C V + ВП — импульс, M = Вт V + АП — кинетический момент, C = mE + Ai, А = J + Л2, m — масса тела, / = (рь — pf)Vg — вес тела в жидкости, g = —9.807 -стандартное ускорение свободного падения, причем g Ц Oaz. Коэффициенты матриц А, В, C могут быть аппроксимированы следующими соотношениями:
cii = С22 = —4.6764 • 10"3pf + pbV + 4.8542 • 10-4pfФ, cas = 74.279 • 10"3pf + рьV — 9.6606 • 10"4pf Ф, aii = a22 = (1.09 74 • 10-3pf + 2.01995 • 10-4рь) +
+ (—1.2418 • 10-5pf + 5.1027 • 10-7рь)Ф, (2.4)
ass = (—3.0673 • 10-4pf + 4.4726 • 10-4рь) + (2.4167 • 10-5pf — 1.9567 • 10-6рь)Ф, &ii = b22 = (—0.51062 + 0.11603Ф — 1.2601 • 10-3Ф2) • 10-3pf, bas = (0.60381 — 0.21822Ф + 2.4818 • 10-3Ф2) • 10-3pf.
Аппроксимации (2.4) и исходные данные для pf = 1000 кг/м3, / =1 Н приведены на рисунке 3.
Система уравнений (2.1), (2.3) полностью описывает движение рассматриваемого тела в идеальной жидкости. Нетрудно заметить, что уравнения относительно p, M, y отделяются.
Согласно работе [6], система (2.1), (2.3) допускает следующие первые интегралы движения:
(p, a) = Pi, (p, в)= P2, (p, Y)+ /t = P3, (2.5)
(M, y) = a,
a2 = в2 = Y2 = 1, (a, в) = (a, Y) = (в, Y) = 0,
где Pi, P2, P3 — проекции начального импульса (начального толчка) на неподвижные оси, a — проекция начального кинетического момента на ось Oaz неподвижной системы координат.
Ф
Ф
Ф
28 24
ч 20
U 16 12
8
60 50 S3 40
о
30
20 -
0
20
1.1 1
0.9 0.8 0.7 0.6
1.6 1.4 1.2 1 0.8
°.6о
Рис. 3. Аппроксимации коэффициентов матриц A, B, C.
40 Ф
60 80
В силу динамической симметрии ац = а22, Ьц = 622, оц = С22 уравнения (2.1), (2.3) допускают интеграл Лагранжа
Мз = C = const. (2.6)
Мы будем рассматривать движение из состояния покоя, в этом случае F\ = P2 = P3 = = 0, а = C = 0.
2.4. Уравнения движения тела в вязкой жидкости
Для полного описания движения тела в вязкой жидкости требуется определение переменных во времени полей давления и скорости жидкости на основании уравнений Навье-Стокса. Такая задача в трехмерной постановке является трудоемкой, а подробное исследование динамики не представляется возможным за разумное время. Расчеты показали, что гидродинамическое сопротивление рассматриваемого винтового тела определяется и поступательной, и угловой скоростями, поэтому для описания вязкого трения используется следующая феноменологическая модель:
Fj = fj \vj\vj + fjК, Gj = gj\vj\vj + gj,
(2.7)
где — ]-я компонента силы сопротивления, Gj — ]-я компонента момента сопротивления, ¡], íjj , 9], — коэффициенты сопротивления.
С учетом силы и момента (2.7) уравнения движения примут вид
p = p х П - HY - Fs, M = M х П + p х V - Gs
(2.8)
где Fs = (^1 ,^2,^3) — стационарная сила вязкого сопротивления, Gs = ^1, G2, Gз) -стационарный момент вязкого сопротивления.
Система уравнений (2.1), (2.8) полностью описывает движение рассматриваемого винта в вязкой жидкости. Как и в случае идеальной жидкости, нетрудно заметить, что уравнения
относительно р, М, 7 отделяются. Кроме того, они обладают следующей симметрией:
Р1 ^ -VI, Р2 ^ -Р2, Рз ^ Рз
Mi - —Mi, Yi - - -7i,
M2 - >-M2, Y2 - - -Y2,
M3 M3, Y3 - » Y3,
t - t.
(2.9)
Для определения коэффициентов /'!■, , д, было смоделировано движение винта, представленного на рисунке 1, при различных поступательных и угловых скоростях на основе уравнений Навье-Стокса и без учета обратного влияния жидкости на движения тела. В результате обработки результатов численных экспериментов были получены следующие зависимости:
fV = fV = (312.9 - 3.54Ф)
- 4.6764 + 0.48542Ф
> 0,
fi = fg = (0.085 - 0.0737Ф) fV = 312.9 - 3.54Ф > 0,
74.279 - 0.96606Ф - 510.62 + 116.03Ф - 1.2601Ф2
603.81 - 218.22Ф + 2.4818Ф2 fg = 0.085 - 0.0737Ф < 0,
> 0,
gl
V_ o o r 510.62 + 116.03Ф - 1.2601Ф2 „
(Jo — —Zo.O--— > U,
603.81 - 218.22Ф + 2.4818Ф2
9i = 92 = (0.542 - 0.0268Ф + 0.0008Ф2)
2 109.74 - 1.2418Ф
- 30.673 + 2.4167Ф
> 0,
gV = -23.5 < 0, gg = 0.542 - 0.0268Ф + 0.0008Ф2 > 0.
На рисунке 4 приведены графики функций (2.7) для Ф = 45 и ( = 0.01.
Сц
ЕС
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
iiiiii - • Численное моделирование — Аппроксимирующая зависимость
0
0.14 0.12 0.1 ^ 0.08 " 0.06 0.04 0.02 0
0.02
0.04 0.06
Vi, V2
0.08 0.1
• Численное моделирование --Аппроксимирующая зависимость
>
0
0.02
0.04 0.06
üJi, U>2
0.08 0.1
1.6 1.4 1.2 1
гЭТ 0.8 ^ 0.6 0.4 0.2 0 0.2
0
0.05 0
-0.05 Ü -0.1 0.15 -0.2 -0.25
i i i i i i - • Численное моделирование — Аппроксимирующая зависимость
/ <
0.02
0.04 0.06
V3
0
0.02
0.08
0.04 0.06
Рис. 4. Значения сил и моментов в зависимости от поступательной скорости при Ф = 45 и (
0.08 0.1
•
« Ч ислс зннс >е
моделирование _ — Аппроксимирующая зависимость ...... i
0.1
= 0.01.
3. Движение в идеальной жидкости
Рассмотрим движение винта в идеальной жидкости. В этом случае движение описывается системой уравнений (2.1), (2.3). Из работ [6, 27] известно, что система уравнений (2.1), (2.3) обладает частными решениями (решения Стеклова) вида
Yk = ±1, Yi = Yj = 0, i = j = k = i,
vk = —№kkt, Uk = —Ц■bkkt, akk
Vi = Vj = 0,
ui = Uj = 0,
Ckk
akk0kk - b2
bkk = -
bkk
(3.1)
kk
akkckk - Ь2
kk
Указанные решения соответствуют равноускоренным движениям с сохранением ориентации одной из осей подвижной системы координат относительно неподвижной системы координат.
Анализ устойчивости решений Стеклова для случая
«11 = «22 = азз = ап, Ъц = 622 = Ьзз = Ьп, еп = С22 = С33 = еп (3.2) проведен в работах [6, 27]. Для случая винта, обладающего динамической симметрией,
aii = a22 = азз, Ьп = 622 = Ьзз, оц = 022 = 033,
подобный анализ ранее не проводился.
В силу симметрии (3.3) можно указать два семейства частных решений:
Yi = Y2 = 0, Yз = ±1, Y2 + Y2 = 1, Yз = 0, Y2 = const, Yз = const.
(3.3)
(3.4)
(3.5)
Решения (3.4) линейно устойчивы при условии сц > Сзз, решения (3.5) — при сц < Сзз. Подробный анализ устойчивости решений (3.4), (3.5) приведен в приложении к статье. Граница устойчивости решений Стеклова показана на рисунке 5.
100
=t 10
20 25 30 35 40 45 50 55 60 Ф
Рис. 5. Граница устойчивости решений (3.4), (3.5). Область 1 соответствует неравенству сц > С33, область 2 соответствует неравенству сц < С33.
1
Исследуем вопрос воспроизведения устойчивости решений (3.4), (3.5) в численных расчетах. Для численного интегрирования уравнений движения был выбран метод Вернера, обеспечивающий шестой порядок точности [28]. Величина шага интегрирования по времени определялась из условия
At = Ato min 1,
1
2п|П|
(3.6)
где Ato — максимальный шаг интегрирования. Вычислительный алгоритм был реализован на языке С++.
Согласно рисунку 5, при j = 1, Ф = 30° решения (3.4) устойчивы. При задании начального угла е = ±10"10 численно воспроизводится устойчивое равноускоренное падение. Если же задать е = ±п/2 ^ 10"10, то в расчетах воспроизводится потеря устойчивости. Поведение величин Yi, Y2 + Y2, Y3 в различные моменты времени приведено на рисунке 6.
1
£ 0 -1
10
15
20
25
0.1 £ 0 -0.1
30
995.9 995.92 995.94 995.96 995.98 996
10
15
20
25
30
995.9 995.92 995.94 995.96 995.98 996
0.8 ^ 0.4 0
10
15
20
25
_ г-
1
0.996
0.992
30
995.9 995.92 995.94 995.96 995.
996
Рис. 6. Поведение величин yi, Yi + Y2, Y3 при j = 1, Ф = 30° и начальном угле е = п/2 - 10
-10
Из рисунка 6 видно, что малое отклонение от положения равновесия приводит к быстрой потере устойчивости, а винт стремится занять положение 73 = 1, соответствующее устойчивому решению (3.4).
Аналогичная ситуация наблюдается при ц = 1, Ф = 60°, однако в этом случае устойчивыми являются решения семейства (3.5). При задании начального угла е = 10_10 численно воспроизводится неустойчивость движения. Поведение величин 71, 7^ + 7!, 73 в различные моменты времени приведено на рисунке 7.
Из рисунка 7 видно, что винт во время движения стремится занять устойчивое положение 73 = 0, соответствующее семейству решений (3.5).
Мы провели поиск устойчивых решений, отличных от (3.4) и (3.5). Для этого при фиксированном значении ц плоскость (Ф,е) покрывалась равномерной сеткой, для каждого узла сетки Ф^, е ^ осуществлялся расчет траектории системы до установления. Каких-либо устойчивых режимов кроме (3.4) и (3.5) найти не удалось.
1
£ 0 -1
10
15
20
25
0.1 0
г-
30
995.9 995.92 995.94 995.96 995.98 996
0
10
1
? 0 -1
15
20
25
СЯ1М
+ 0.99
(МгЧ
30
0.98
995.9 995.92 995.94 995.96 995.98 996
0.1 £ 0 -0.1
10 15 20 25 30 995.9 995.92 995.94 995.96 995.98 996
Рис. 7. Поведение величин 71, 72 + 7!, 73 при / =1, Ф = 60° и начальном угле е = 10~10.
4. Движение в вязкой жидкости
4.1. Частные решения
Рассмотрим движение винта в вязкой жидкости, в этом случае движение описывается системой уравнений (2.1), (2.8). Система уравнений (2.1), (2.8) допускает частные решения
71 = 1, 72 = 7з = 0,
У1 =
/93
ГО пш ¡191 - ¡1 91
, У2 = Уз = 0, Ш1 =
т
ГО „Ш ¡191 - ¡191
, Ш = Ш = 0,
71 = 72 = 0, 7з = 1,
У1 = У2 = 0, Уз = -4
/93
¡з 9з — ¡39] '
Ш = Ш2 = 0, Шз = —
/9 ]
(4.1)
(4.2)
/з 93 — ¡3 9]
соответствующие движениям с постоянными поступательной и угловой скоростями при сохранении ориентации одной из осей подвижной системы координат относительно неподвижной системы координат.
Анализ устойчивости решений (4.1), (4.2) аналитическими методами требует поиск корней характеристического многочлена шестой степени, что весьма затруднительно. Однако численные расчеты показывают, что оба решения (4.1), (4.2) неустойчивы.
В качестве примера рассмотрим процесс потери устойчивости решения (4.1) при / = 1, Ф = 30° (см. рис. 8)
Из рисунка 8 видно, что движение быстро выходит на стационарный режим (при Ь ~ 5 с). Поступательная и угловая скорости стационарного движения воспроизводятся с точностью порядка 10_15, что соизмеримо с величиной машинной погрешности. По истечении достаточно короткого промежутка времени (при Ь ~ 72 с) после выхода на стационарный режим движения происходит потеря устойчивости. После некоторого переходного процесса скорость У1 изменяется периодически.
1
0.20 0.10 » 0.00 -0.10 -0.20
0
20 40 60
Рис. 8. Зависимость скорости от времени £ при л =1, Ф = 30° и начальном угле е = п/2.
80 100 120 140 160 180 200 t
Рассмотренный пример мотивирует исследование различных предельных режимов, возникающих в системе.
4.2. Карта показателей Ляпунова. Бифуркационные диаграммы
В системах с диссипацией [14, 17] могут существовать различные притягивающие режимы (как периодические, так и хаотические). Для исследования возможных режимов движения была построена карта динамических показателей Ляпунова [31, 32] на плоскости параметров л, Ф (см. рис. 9).
Рис. 9. Карта режимов 400 х 400 узлов.
При построении карты режимов плоскость (^ л, Ф) покрывалась равномерной сеткой 400 х 400 узлов. Для каждого узла сетки производился расчет траектории и оценка показателей Ляпунова. Все наблюдаемые в системе режимы являются притягивающими (сумма показателей отрицательна). Если старший показатель равен нулю, то наблюдается предельный цикл, а соответствующий узел отмечается оттенком серого". Интенсивность рассчитывалась по формуле
А2 - А™п
gray
1 -
\ max \ mm ' Л2 — Л2
(4.3)
аДля читателя печатной версии: здесь и далее полноцветные версии рисунков см. в эл. версии статьи — http://nd.ics.org.ru/nd1601007/
где Л2 — второй показатель Ляпунова, А™1П и А™ах — наименьшее и наибольшее, соответственно, значения второго показателя в области регулярного движения, оцененные численно.
Если старший показатель положителен, то наблюдается хаотический аттрактор, а соответствующий узел сетки отмечается оттенком красного. Интенсивность рассчитывалась по формуле
!гв4 1
А1 - А™ш
Атах
Ашт '
(4.4)
где А1 — старший показатель Ляпунова, А^"1" и Атах — наименьшее и наибольшее, соответственно, значения старшего показателя в области регулярного движения, оцененные численно.
Анализ карты режимов показывает, что рассматриваемое винтовое тело при законе вязкого сопротивления (2.7) совершает либо регулярное движение, либо хаотическое. Все реализующиеся режимы являются асимптотически устойчивыми. Далее рассмотрим более подробно возможные режимы движения.
4.2.1. Регулярные режимы
При ц = 3 и Ф = 40° реализуется периодическое движение. Тело движется по спирали достаточно сложной формы. Проекция траектории центра масс тела на плоскость Оху показана рисунке 10; видно, что движение в плоскости Оху ограничено. Скорость падения тела (V, 7) непостоянна (см. рис. 11).
-0.18 -0.20 -0.22 -0.24
со
^ -0.26 -0.28 -0.30
-0.32 500
д д д д д д д д
510
520
530
540
550
Рис. 10. Проекция траектории центра масс на плоскость Оху.
Рис. 11. Скорость падения ( V, 7).
х
Вращение тела в пространстве является квазипериодическим движением с круговыми частотами ш ^ 0.56549 с 1 и ш2 ~ 0.23038 с 1. На рисунке 12 показан график изменения аппликаты апекса оси Оез = (аз,вз,1з) во времени. Видно, что величина 73 регулярно изменяет знак, это свидетельствует о периодических переворотах тела.
На рисунке 13 показана зависимость ф от времени. Из рисунка 13 видно, что тело совершает регулярную прецессию.
4.2.2. Бифуркации и переход к хаосу
Для исследования бифуркаций, приводящих к возникновению хаоса, были построены однопараметрические бифуркационные диаграммы неподвижных точек сечения Пуанкаре.
1.00 0.80 0.60 0.40 0.20 ¡2 0.00 -0.20 -0.40 -0.60 -0.80 -1.00
500 510
520 530
t
Рис. 12. Аппликата апекса оси Овз.
540 550
4.00 3.50 3.00 2.50 2.00 1.50 1.00 0.50
500 510
520 530
t
540
550
Рис. 13. Скорость падения ( V, y).
0.10 0.00 -0.10 -0.20 -0.30 -0.40 -0.50 -0.60 -0.70
Ai Л2 1 4........ - Аз----- , , , , , , , , , , , ,
7"7?- У ;' Wf? /ЛА ^ -- ___
^ '' ' \ ,л / \ / \
|;; 'i:/ ■ } j ■ V v' 1 i J Л .
' 1 ! 1J 1 I t ' 1 г 1 ' \ ' и —
\ i' 1'' 1
. . . . . . . . 1 ' . . . . . . . , , , ,
30
32
¡F 0.00
30
32
34
34
36 Ф
36 Ф
38
38
40
40
42
42
Рис. 14. Зависимость показателей Ляпунова от параметра Ф и бифуркационная диаграмма для переменной 73. Маркеры соответствуют точкам бифуркаций.
Сечение Пуанкаре рассчитывалось по методу Эно [30, 31]. В качестве секущей в фазовом пространстве была выбрана поверхность рз = 0.
При расчете бифуркационных диаграмм применяют схему с наследованием начальных условий для ускорения сходимости [33], однако численные расчеты показали, что при одних и тех же значениях параметров сосуществуют различные аттракторы. В связи с этим при построении бифуркационной диаграммы осуществлялось сканирование с различными начальными условиями е € [0,^/2].
На рисунке 14 показана бифуркационная диаграмма для переменной 73 и зависимость первых трех показателей Ляпунова от параметра Ф. На диаграмме 14 маркерами обозначены точки, в которых система терпит бифуркацию «суперкритическая вилка».
На рисунке 15 показана бифуркационная диаграмма для значений Ф € [45.4, 45.68]. На диаграмме 15 метками «1» и «2» обозначены ветки, соответствующие различным режи-
0.05 0.00 -0.05 -0.10 -0.15 -0.20 -0.25 -0.30 -0.35
Ах-
45.4
0.12 0.10 0.08 0.06 0.04 0.02
-Л2----Лз
45.45
45.5
45.55 Ф
45.4
45.45
45.5
45.55 Ф
45.6
45.6
45.65
1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 : 2
: тзйЯВтИР—--
•лящщЩ
бифуркация «седло—узел» х 'Я
бифуркация удвоения периода • , , , , 1 , .^Гт1®, ,
45.65
Рис. 15. Бифуркационная диаграмма для переменной 73. Маркеры соответствуют точкам бифуркаций.
мам движения. Режим «2» рождается при Ф ~ 45.51 вследствие седло-узловой бифуркации. При движении по параметру Ф слева направо режим «2» претерпевает каскад бифуркаций удвоения периода, в результате которого происходит переход к хаосу. Режим «2» перестает существовать при Ф ~ 45.658.
Проекции аттракторов, соответствующие различным режимам движения при Ф = 45.64, показаны на рисунках 16 и 17.
Траектория аттрактора, изображенного на рисунке 16, обладает следующими показателями Ляпунова:
Л1 = 2.39 • 10"4 ± 6.87 • 10"4, Л2 = -1.56 • 10"3 ± 4.74 • 10_3, Л3 = -0.109 ± 1.29 • 10"3, Л4 = -0.664 ± 1.62 • 10"3, Л5 = -0.664 ± 9.28 • 10"4, (4.5) Л6 = -1.15 ± 4.74 • 10"3, Л7 = -1.18 ± 1.58 • 10"4, Л8 = -1.35 ± 1.28 • 10"4.
Оба старших показателя являются нулевыми, что свидетельствует о квазипериодичности движения.
Траектория аттрактора, изображенного на рисунке 17, обладает следующими показателями Ляпунова:
Л1 = 2.73 • 10"2 ± 2.45 • 10"3, Л2 = -1.48 • 10"3 ± 4.49 • 10"3, Л3 = -0.287 ± 6.94 • 10"3, Л4 = -0.493 ± 6.24 • 10"3, Л5 = -0.829 ± 3.49 • 10"3, (4.6)
4
Л6 = -1.02 ± 3.81 • 10"3, Л7 = -1.19 ± 5.41 • 10"4, Л8 = -1.35 ± 3.23 • 10
Старший показатель является положительным, что свидетельствует о хаотичности движения.
Рис. 16. Проекции простого аттрактора при ¡л = 3, Ф = 45.64°.
5. Заключение
Проведенное исследование показало, что модели идеальной и вязкой жидкостей дают качественно различные результаты в расчетах свободного падения винтового тела. При описании движения в рамках теории идеальной жидкости в системе наблюдаются два простых аттрактора, соответствующие равноускоренным падениям. Применение модели вязкой жидкости приводит к возникновению в системе хаоса при некоторых значениях параметров. Переход к хаосу происходит через каскад бифуркаций удвоения периода. Кроме того, при одних и тех же значениях параметров могут сосуществовать несколько устойчивых режимов движения.
Приведем ряд открытых вопросов, касающихся динамики тел винтовой формы:
1. Экспериментальные исследования свободного падения винта и верификация теоретических результатов, полученных в рамках моделей идеальной и вязкой жидкостей.
2. Исследование движений уравновешенного винта при наличии начального толчка.
3. Исследование стабилизации движений тела за счет вращения внутренних роторов.
4. Исследование движений винта, управляемых вращением внутренних роторов, перемещением внутренних точечных масс и изменением угла разворота лопастей.
Авторы благодарят А. В. Борисова, И. С. Мамаева, А. А. Килина, С. П. Кузнецова за ценные обсуждения результатов.
Приложение. Устойчивость решений Стеклова
Уравнения (2.1), (2.3) обладают частными решениями [6, 27]
1к = ±1, И = Ц = 0, г = 2 = к = г, (П.1)
называемыми решениями Стеклова и соответствующими равноускоренным падениям с равноускоренным вращением. Для случая различных коэффициентов
«11 = «22 = азз = ап, Ъц = Ь22 = Ьзз = Ъц, еп = С22 = С33 = еп (П.2)
вопрос устойчивости исследован в работах [6, 27]. Однако для винта, обладающего динамической симметрией вида
ац = а22 = азз, Ъц = Ъ22 = Ьзз, еп = С22 = сзз, (П.3)
вопрос устойчивости решений (П.1) не изучался.
Представим уравнения (2.1), (2.3) в гамильтоновой форме
м х дм+рх др' р-рх дм ад 7- = (- - = § > - = (П4)
Я = ± (лм, м) + (вм, р) + | (ср, р) + /лхз, А = (А - БтС-1Б)-1, С = (С - БА-1БТ)-1, Б = -СВА-1.
На нулевом уровне интеграла импульса (2.5)
(р, а) = (р, 0) = (р, 7) + ^ = 0 уравнения движения (П.4) могут быть редуцированы к виду
—ж. —ж. (п-5)
Я = ± ^АМ, М) - fit [ВМ, 7J + ^2i2 (С7, 7
Из системы (П.4) видно, что уравнения относительно Ж, 7 отделяются. В последующем анализе мы будем изучать изменение только величин М, 7. Запишем уравнения (П.5) в явном виде
Mi = (азз - аи)М2Мз + ¡t(bn - Ьзз)М73 + M3Y2) + /J.2t2(c33 - си)y273, M2 = (ап - азз)МзМ1 + ¡t(&33 - М(Мз71 + M173) + ¡2t2(cn - C33)7зYi,
M3 = 0,
- 3 , г ^ (П.6)
71 = аззM3Y2 - aiiM2J3 + ¡t(bn - Ьзз)Y2Y3, Y2 = anM1Y3 - аззM3Y1 + ¡t(b33 - bn )Y3Y1,
Y3 = an(M2Y1 - M1Y2).
Здесь первое уравнение выражает интеграл Лагранжа, являющийся следствием динамической симметрии [25]. Так как мы рассматриваем движение из состояния покоя, то M3 = 0, а уравнения (П.6) примут вид
Mi = цфи - Ьзз)M2Y3 + f2t2 (С33 - СЦ)ъъ, M2 = Kc33 - Cii)Mi73 + H2t2(cii - С33)737i,
Yi = -CiiM273 + /J.t(bii - b33)73 72, (п.7)
Y2 = CiiMi73 + Kc33 - Cii)Y3Yi, 73 = Сц (M2Yi - Mi72).
Решения Стеклова для системы (П.7) имеют вид
Yi = ±1, 72 = 73 = 0, Mi = M2 = 0, (П.8)
72 = ±1, 7i = 73 = 0, Mi = M2 = 0, (П.9)
73 = ±1, 7i = 72 = 0, Mi = M2 = 0. (П.10)
Замечание 1. Вследствие динамической симметрии (П.3) решения (П.8), (П.9) идентичны и принадлежат семейству решений вида
7i + 72 = 1, 73 = 0, 71 = const, 72 = const, M2 = M3 = 0. (П.11)
В дальнейшем для анализа устойчивости будем пользоваться теоремой [29] об асимптотическом поведении систем вида
Z = ( A + V(t) + R(t)) Z. (П.12)
Теорема 1. Пусть A — постоянная матрица с различными характеристическими корнями fj, j = 1,...,n. Пусть матрица V дифференцируема и удовлетворяет условию
J\V'(t)\dt< œ (П.13)
0
и пусть V(t) ^ 0 при t ^œ. Пусть матрица R интегрируема и
У |R(t)|dt < œ. (П.14)
0
Обозначим корни уравнения det(A + V(t) - AE) = 0 через Xj(t), j = 1,...,n. Очевидно, что
можно, если это необходимо, переставить fj так, чтобы lim Aj (t) = fj. Для данного k
t—
положим,
Dkj(t) = Re(Afc(t) - Aj (t)). (П.15)
Допустим, что все j, 1 ^ j ^ n, попадают в один из двух классов Ii и I2, где
t
j G Ii, если J Dkj (r )dr ^ œ при t ^œ (П.16)
¿2
У Бу(т)йт > -К (¿2 ^ ¿1 ^ 0), (П.17)
¿1
¿2
3 € 12, если J (т)йт <К (¿2 ^ ¿1 ^ 0); (П.18)
¿1
.здесь к фиксировано и К — постоянная. Пусть рк — характеристический вектор А, соответствующий Лк, так что
Арк = Лк Рк. (П.19)
Тогда существует решение фк системы (П.12) и число ¿о, 0 ^ ¿о ^ ж, такие, что
рк. (П.20)
Нш фк (¿) ехр
Ь—>оо
-/^к (т )йт
. ¿0
П.1. Устойчивость решения 73 = 1
Линеаризуем систему уравнений (П.7) в окрестности решения (П.10), тогда система относительно возмущений примет вид
¿71 = 1лг(Ъп - Ъзз)5^2 - Т1116М2, ¿72 = КЬзз - Ьц)^71 + ап5М1,
(П.21)
¿ММ 1 = Л2^2(сзз - сц)5ъ + лФп - Ъзз)ЬМ2, ¿ММ2 = ц212(Ъи - сзз)¿7l + ц,г(Ъзз - Ьц^Мь Вблизи решения 7з = 1 в силу интеграла 72 = 1 можно записать следующее выражение:
7з = 1 - ¿73 + о (¿73), ¿7з = т^(¿7? + ¿71) С помощью замены переменных и масштабирования времени
£ = (¿¿71, ¿¿72, 5М\, ¿А/, )т , т = \{2 (П.22)
система уравнений (П.21) может быть представлена в виде, идентичном виду (П.12),
Щ = + (п.23)
йт 2т
М = (р1 °1\ Н = Ц, Е1 = С2 = л ( 0 Ъ11 - Ъзз
\Р2 с^ у 0 2 0^ \Ъзз -Ъ11 0
С1 = | 0 -С11), ^ = Л2 ( V Сза-С11 ас11 0 с11 - сзз 0
где Е2 — единичная матрица размерности 2, 02 — нулевая матрица размерности 2.
и
Согласно теореме 1, поведение системы (П.23) при Ь будет определяться собствен-
ными числами матрицы М. Решение системы (П.23) будет устойчивым, если все собственные числа матрицы М будут иметь неположительную действительную часть. Характеристический многочлен матрицы М имеет форму биквадратного уравнения
Л4 + аЛ2 + Ь = 0,
(П.24)
а = 2^ ПЬзз -Ьи) + а\1 (сц - сзз)),
Ь = (аи(сзз - сп) + (Ьзз - Ьп)2
Следовательно, для того чтобы решение системы (П.23) было устойчивым, все корни уравнения (П.24) должны лежать на мнимой оси. Для этого необходимо и достаточно выполнение трех неравенств а > 0, Ь > 0, а2 — 4Ь > 0. Нетрудно убедиться, что все три неравенства выполняются при условии
ац > Сзз. (П.25)
Таким образом, неравенство (П.25) является условием устойчивости решения (П. 10).
П.2. Устойчивость решения = 1
Линеаризуем систему уравнений (П.7) в окрестности решения (П.8), тогда система относительно возмущений примет вид
¿7з = ац Ш2, ¿72 = цЬ{Ьзз - С11 М7з, 5Ы1 = 0, Ш2 = 1л2Ь2(с11 - Сзз)5^з.
(П.26)
Вблизи решения 71 = 1 в силу интеграла 72 = 1 можно записать следующее выражение:
71 = 1 - ¿71 + о(5Ъ), ¿71 = + ¿71) 0.
С помощью замены переменных и масштабирования времени
Я ¿.\/,)' . т =
12
(П.27)
система уравнений (П.26) может быть представлена в виде, идентичном виду (П.12),
йт 2т
(
М
\
0 0 С1
» (Сзз -Сп) 0 0 (Сц - Сзз) 0 0 у
(
Н
1 0 0 0 1 0 0 0 0
\
(П.28) (П.29)
Согласно теореме 1, поведение системы (П.28) при Ь ^ж будет определяться собственными числами матрицы М. Решение системы (П.28) будет устойчивым, если все собственные числа матрицы М будут иметь неположительную действительную часть. Характеристический многочлен матрицы М имеет форму кубического уравнения
Следовательно, для того чтобы решение системы (П.28) было устойчивым, все корни уравнения (П.30) должны лежать на мнимой оси. Для этого необходимо и достаточно выполнение неравенства
Таким образом, неравенство (П.31) является условием устойчивости решения (П.8).
Список литературы
[1] Кирхгоф Г. Механика. Лекции по математической физике. Москва: АН СССР, 1962. 404 с.
[2] Стеклов В. А. О некоторых возможных движениях твердого тела в жидкости // Тр. отд. физ. наук общ-ва любителей естествознания, антропологии и этнографии, 1895, т. 7, №2, с. 10-21.
[3] Ламб Г. Гидродинамика. Москва-Ленинград: ОГИЗ, 1947. 929с.
[4] Стеклов В. А. О движении тяжелого твердого тела в жидкости: 1 // Сообщ. Харьковск. матем. общ-ва. Сер. 2, 1891, т. 2, №5-6, с. 209-235.
[5] Чаплыгин С. А. О движении тяжелых тел в несжимаемой жидкости // Полн. собр. соч. Москва-Ленинград: АН СССР, 1933. Т. 1, с. 133-150.
[6] Borisov A. V., Kozlov V. V., Mamaev I. S. Asymptotic stability and associated problems of dynamics of falling rigid body // Regul. Chaotic Dyn., 2007, vol. 12, no. 5, pp. 531-565.
[7] Mougin G., Magnaudet J. The generalized Kirchhoff equations and their application to the interaction between a rigid body and an arbitrary time-dependent viscous flow // Int. J. Multiphas. Flow, 2002, vol.28, no. 11, pp. 1837-1851.
[8] Childress S., Spagnolie S.E., Tokieda T. A bug on a raft: Recoil locomotion in a viscous fluid // J. Fluid Mech., 2011, vol.669, pp. 527-556.
[9] Рамоданов С. М., Тененев В. А. Движение тела с переменной геометрией масс в безграничной вязкой жидкости // Нелинейная динамика, 2011, т. 7, №3, с. 635-647.
[10] Ветчанин Е.В., Мамаев И. С., Тененев В. А. Движение тела с переменной геометрией масс в вязкой жидкости // Нелинейная динамика, 2012, т. 8, №4, с. 815-836.
[11] Vetchanin E.V., Mamaev I.S., Tenenev V.A. The self-propulsion of a body with moving internal masses in a viscous fluid // Regul. Chaotic Dyn., 2013, vol. 18, nos. 1-2, pp. 100-117.
[12] Maxwell J. K. On a particular case of descent of a heavy body in a resisting medium // Cambridge and Dublin Math. Journ., 1854, vol. 9, pp. 145-148.
[13] Козлов В. В. К задаче о падении тяжелого твердого тела в сопротивляющейся среде // Вестн. Моск. ун-та. Сер. Матем. Механ., 1990, №1, с. 79-86.
[14] Tanabe Y., Kaneko K. Behavior of a falling paper // Phys. Rev. Lett., 1994, vol.73, no. 10, pp.1372-1375.
[15] Belmonte A., Eisenberg H., Moses E. From flutter to tumble: Inertial drag and froude similarity in falling paper // Phys. Rev. Lett., 1998, vol.81, no. 2, pp. 345-348.
[16] Andersen A., Pesavento U., Wang Z.J. Analysis of transitions between fluttering, tumbling and steady descent of falling cards // J. Fluid Mech., 2005, vol.541, pp. 91-104.
A3 - Ла1Л2(си - сзз) = 0.
(П.30)
C11 < C33.
(П.31)
[17] Кузнецов С. П. Движение падающей пластины в жидкости: Конечномерные модели и феномены сложной нелинейной динамики // Нелинейная динамика, 2015, т. 11, № 1, с. 3-49.
[18] Leipnik R. B., Newton T.A. Double strange attractors in rigid body motion with linear feedback control // Phys. Lett. A, 1981, vol.86, no. 2, pp. 63-67.
[19] Акуленко Л. Д., Лещенко Д. Д., Черноусько Ф. Л. Быстрое движение вокруг неподвижной точки тяжелого твердого тела в сопротивляющейся среде // МТТ, 1982, №3, с. 5-13.
[20] Черноусько Ф. Л., Акуленко Л. Д., Лещенко Д. Д. Эволюция движений твердого тела относительно центра масс. Москва-Ижевск: Институт компьютерных исследований, 2015. 308 с.
[21] Inarrea M., Lanchares V. Chaotic pitch motion of an asymmetric non-rigid spacecraft with viscous drag in circular orbit // Int. J. Nonlinear Mech., 2006, vol.41, no. 1, pp. 86-100.
[22] Кошляков В. Н. Задачи динамики твердого тела и прикладной теории гироскопов: Аналитические методы. Москва: Наука, 1985. 288 с.
[23] Шлихтинг Г. Теория пограничного слоя. Москва: Наука, 1974. 712 с.
[24] Воинов О. В. Инерциальное движение тела в идеальной жидкости из состояния покоя // Прикл. мех. и техн. физ., 2008, т. 49, №4, с. 214-219.
[25] Борисов А. В., Мамаев И. С. Динамика твердого тела: гамильтоновы методы, интегрируемость, хаос. Москва-Ижевск: Институт компьютерных исследований, 2005. 576 с.
[26] Маркеев А. П. Теоретическая механика. Москва-Ижевск: НИЦ «Регулярная и хаотическая динамика», 2007. 592 с.
[27] Дерябин М. В. Об устойчивости равноускоренных вращений тяжелого твердого тела в идеальной жидкости // МТТ, 2002, №5, с. 30-34.
[28] Hairer E., N0rsett S. P., Wanner G. Solving ordinary differential equations: Vol. 1. Nonstiff problems. 2nd ed. (Springer Ser. Comput. Math., vol.8.) Berlin: Springer, 1993. 528pp.
[29] Коддингтон Э. А., Левинсон Н. Теория обыкновенных дифференциальных уравнений. Москва: ИЛ, 1958. 475 с.
[30] Henon M. On the numerical computation of Poincare maps // Phys. D, 1982, vol. 5, nos. 2-3, pp.412-414.
[31] Кузнецов С. П. Динамический хаос. 2-е изд. Москва: Физматлит, 2006. 356 с.
[32] Borisov A.V., Kazakov A. O., Sataev I.R. The reversal and chaotic attractor in the nonholonomic model of Chaplygin's top // Regul. Chaotic Dyn., 2014, vol. 19, no. 6, pp. 718-733.
[33] Борисов А. В., Казаков А. О., Кузнецов С. П. Нелинейная динамика кельтского камня: Него-лономная модель // УФН, 2014, т. 184, №5, с. 493-500.
Chaotic dynamics in the problem of the fall of a screw-shaped body in a fluid
Valentin A.Tenenev1, Evgeny V. Vetchanin2, Lenar F. Ilaletdinov3
i,2,3izhevsk State Technical University Studencheskaya 7, Izhevsk, 426069 Russia 2 Udmurt State University Universitetskaya 1, Izhevsk, 426034 Russia
1tenenev@istu.ru, 2eugene186@mail.ru, 3ano-lenar@yandex.ru
This paper is concerned with the process of the free fall of a three-bladed screw in a fluid. The investigation is performed within the framework of theories of an ideal fluid and a viscous fluid. For the case of an ideal fluid the stability of uniformly accelerated rotations (the Steklov solutions) is studied. A phenomenological model of viscous forces and torques is derived for investigation of the motion in a viscous fluid. A chart of Lyapunov exponents and bifucation diagrams are
120
В. A. Тененеe, Е. В. BemnaHUH, JI. H^a^emduHoe
computed. It is shown that, depending on the system parameters, quasiperiodic and chaotic regimes of motion are possible. Transition to chaos occurs through cascade of period-doubling bifurcations.
MSC 2010: 70E15, 65Pxx
Keywords: ideal fluid, viscous fluid, motion of a rigid body, dynamical system, stability of motion, bifurcations, chart of Lyapunov exponents
Received February 10, 2016, accepted March 04, 2016
Citation: Rus. J. Nonlin. Dyn., 2016, vol. 12, no. 1, pp. 99-120 (Russian)
