Министерство науки и высшего образования Российской Федерации
Федеральное государственное бюджетное образовательное учреждение
высшего образования
Пермский национальный исследовательский политехнический университет
Факультет
Выпускающая кафедра:
Направление подготовки:
Квалификация:
Прикладной математики и механики
Динамика и прочность машин
15.03.03 «Прикладная механика»
бакалавр
Допускается к защите
Зав.кафедрой ДПМ
__________В.П.Матвеенко
«21» июня 2019г.
ИЗУЧЕНИЕ ЭВОЛЮЦИИ ФОРМЫ ГИБКОГО УПРУГОГО СТЕРЖНЯ
ПРИ ЕГО КВАЗИСТАТИЧЕСКОМ КРУЧЕНИИ
Выпускная квалификационная работа
Выполнила:
студентка группы ДПМ-15-1б
Еленская Наталия Витальевна
(___________)
подпись
Научный руководитель:
д.ф.-м.н., проф. кафедры ДПМ
Келлер Илья Эрнстович
(___________)
подпись
Пермь 2019
Министерство науки и высшего образования Российской Федерации
Федеральное государственное бюджетное образовательное учреждение
высшего образования
Пермский национальный исследовательский политехнический университет
Кафедра «Динамика и прочность машин»
УТВЕРЖДАЮ
Зав. кафедрой ДПМ
Матвеенко В.П.
«05» июня 2019г.
ЗАДАНИЕ
на выполнение выпускной квалификационной работы студента
Фамилия И.О.
Еленская Н.В.
Факультет ПММ
Группа ДПМ-15-1б
Начало выполнения работы
05 июня 2019
Контрольные сроки просмотра работы кафедрой 18 июня 2019
Дата защиты работы на заседании ГЭК 25 июня 2019
1. Наименование работы: Изучение эволюции формы гибкого упругого стержня при его
квазистатическом кручении;
2. Исходные данные к работе
3. Содержание пояснительной записки:
1) Математическая постановка задачи;
2) Исследование статических уравнений нелинейной системы уравнений
равновесия гибкого упругого стержня;
3) Численное решение системы уравнений в среде MATLAB
Дополнительные указания
Основная литература
1.
Bachurikhin V.P., Keller I.E., Merzlyakov A.F., Yurlov M.A. Experimental study of
nonlinear effects under torsion of the uniform cylinder with initially circular cross section. Solid
State Phenomena. Trans Tech Publications, Switzerland. 2016. Vol.243. P.29-34.
2.
Светлицкий В.А. Механика гибких стержней и нитей. —М.:Машиностроение, 1978.
— 222 с, ил. — (Б-ка расчетчика).
Руководитель выпускной квалификационной работы студента
д. ф.-м. н., профессор кафедры ДПМ
(Келлер И. Э.)
(должность, Ф.И.О.)
Задание получил
05 июня 2019
(Еленская Н. В.)
(дата и подпись студента)
КАЛЕНДАРНЫЙ ГРАФИК ВЫПОЛНЕНИЯ
ВЫПУСКНОЙ КВАЛИФИКАЦИОННОЙ РАБОТЫ
Объем
этапа,
в%
№п.
п.
1
2
3
4
Разработка
основных
разделов
выпускной квалификационной работы
Оформление
выпускной
квалификационной работы
Разработка
и
оформление
иллюстративной части материала к
защите диссертации
Представление ВКР на проверку и отзыв
научного руководителя
5
Представление
кафедрой
работы
6
Защита на заседании ГЭК
Сроки выполнения
начало
Примечание
конец
70
20
10
заведующему
Научный руководитель работы:
д.ф.-м.н., профессор кафедры ДПМ
«_____» __________________2019 г.
И. Э. Келлер
ОГЛАВЛЕНИЕ
ВВЕДЕНИЕ .............................................................................................................. 5
ГЛАВА 1. ПОСТАНОВКА ЗАДАЧИ ................................................................... 8
1.1. Наблюдаемые нелинейные эффекты ........................................................ 8
1.2. Основные положения механики гибких стержней ............................... 10
1.3. Вывод нелинейной системы равновесия гибкого упругого стержня для
общего случая статической неопределимости ............................................... 11
ГЛАВА
2
ИССЛЕДОВАНИЕ
СТАТИЧЕСКИХ
УРАВНЕНИЙ
НЕЛИНЕЙНОЙ СИСТЕМЫ УРАВНЕНИЙ РАВНОВЕСИЯ ГИБКОГО
УПРУГОГО СТЕРЖНЯ ........................................................................................ 23
2.1. Общий случай решения статических уравнений системы ................... 23
2.2. Частные решения системы для случая постоянной крутки стержня .. 25
2.3. Частный случай эволюции формы стержня в виде однородной
цилиндрической спирали .................................................................................. 31
ГЛАВА 3 ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМЫ УРАВНЕНИЙ В СРЕДЕ
MATLAB ................................................................................................................ 36
3.1. Использование численных методов для решения нелинейных систем
дифференциальных уравнений......................................................................... 36
3.2. Численное решение нелинейной системы уравнений равновесия
гибкого упругого стержня в случае постановки начальной задачи ............. 36
3.3. Численное решение нелинейной системы уравнений равновесия
гибкого упругого стержня в случае постановки краевой задачи ................. 41
ЗАКЛЮЧЕНИЕ ..................................................................................................... 48
СПИСОК ЛИТЕРАТУРЫ..................................................................................... 49
ПРИЛОЖЕНИЕ А ................................................................................................. 50
ПРИЛОЖЕНИЕ Б .................................................................................................. 51
ПРИЛОЖЕНИЕ В ................................................................................................. 52
ПРИЛОЖЕНИЕ Г .................................................................................................. 53
Введение
Гибкие стержни широко применяют в различных областях техники: их
используют в качестве упругих элементов различных приборов, например, в
качестве
чувствительных
элементов
в
акселерометрах
и
частотных
преобразователях, механических низкочастотных фильтров в электронной
технике, аккумуляторов механической энергии в часовых механизмах и др.
Приборы времени, использующие гибкие стержни, широко распространены
не только как часы, но и как преобразователи стабильных сигналов в
различных устройствах автоматики[1].
Целью
данной
работы
является
формулировка
уравнений
математической модели эволюции формы гибкого стержня при его кручении,
ее численная реализация и верификация. Работа направлена на создание
инструмента для исследования нелинейных явлений, обнаруженных в
эксперименте.
Объект
исследования:
эволюция
формы
гибкого
стержня
с
закрепленными концами при его кручении.
Современное состояние исследования:
В
статье
[2]
обсуждается
кинетическая
аналогия
Кирхгофа,
связывающая основные уравнения для статики упругих стержней и динамики
твердых тел в свете нескольких различных гамильтоновых формулировок,
включая неканоническое описание равновесий стержней. Особое внимание
уделяется восстановлению осевой линии стержня из переменных каркаса,
которые образуют полное пространство конфигурации в интерпретации
твердого тела.
В статье [3] рассматривается геометрическая теория стержней для
случая естественно прямого, линейно упругого, нерастяжимого, круглого
стержня, испытывающего изгиб и кручение, но не сдвигающегося.
Исследуется поведение таких стержней под воздействием конечного момента
5
и натяжения, явно вычисляются геометрические свойства изогнутых
стержней. В едином подходе, использующем динамическую аналогию
Кирхгофа, рассматриваются как классическая спиральная, так и более
недавно исследованная локализованная потеря устойчивости. Особое
внимание
уделяется
последовательной
трактовке
концепций
связи,
скручивания и изгиба.
В статье [4] показано, что равномерный и гиперэластичный, но в
любом случае произвольный нелинейный стержень Коссера, подверженный
соответствующим конечным нагрузкам, имеет равновесия, центральные
линии которых образуют двухпараметрические семейства спиралей.
В статье [5] приведены результаты экспериментов по прогибу,
включающих изгиб и скручивание стержней. Проводится сравнения с
расчетами, основанными на теории стержней Коссера. Стержни зажаты в
совмещенных патронах, и эксперименты проводятся в жестких условиях
нагрузки. Эксперимент продолжается либо скручиванием концов стержня на
определенную величину, а затем корректировкой провисания или фиксацией
провисания
и
изменением
величины
закручивания. Таким
образом,
исследуются часто встречающиеся явления, такие как изгиб с защелкой,
образование петель и выпучивание в и из плоских конфигураций. Эффект
гравитации обсуждается.
В статье [6] представлены результаты испытаний на кручение
сплошных цилиндрических образцов при фиксированном расстоянии между
концами образца. Обнаружено возникновение осевого напряжения, сначала
растягивающего, а затем – сжимающего образец, имеющего пульсационную
составляющую с периодом один оборот. Представленные в работе
экспериментальные данные обращают внимание на нелинейные эффекты,
сопровождающие кручение сплошных металлических цилиндров при
больших деформациях.
6
Задачи:
1.
Вывод нелинейной системы уравнений равновесия гибкого упругого
стержня;
2.
Построение частного решения уравнений в форме однородной
цилиндрической спирали для тестирования программы;
3.
Написание программы в пакете MATLAB для численного решения
краевых задач на кручение гибких стержней;
4.
Численное решение ряда начальных и граничных задач на кручение
гибкого стержня, в том числе воспроизведение решения в виде однородной
цилиндрической спирали.
Научная новизна: заключается в том, что впервые была решена
краевая задача с конкретными заданными условиями.
Практическая значимость: задача имеет фундаментальную ценность.
Объем и структура работы. Работа состоит из введения, трех глав и
заключения. Полный объём работы составляет 49 страниц, включая 22
рисунка и 0 таблиц. Список литературы содержит 10 наименований.
7
Глава 1. Постановка задачи
1.1.
Наблюдаемые нелинейные эффекты
При достаточно больших степенях деформации кручение сплошных
цилиндров с начально круговым поперечным сечением сопровождается
нелинейными эффектами — изменением длины (появлением осевого
напряжения) и развитием спиральной формы рабочей части образца[6].
Первый эффект называется прямым эффектом Пойнтинга и заключается в
следующем: при кручении образца из изотропного материала для сохранения
неизменной длины возникает сжимающая осевая сила.
При достаточно сильном закручивании одного из концов гибкого
упругого стержня (стержня Коссера), форма линии стержня сильно
изменяется и возникают такие эффекты, как закручивание в однородную
цилиндрическую спираль, появление петли[5] и т.д.
Рис. 1 Эволюция формы скрученного стержня[5]
8
В статье [6] был проведен опыт на квазистатическое кручение
цилиндрических образцов с диаметром рабочей части 10 мм и длиной 112 мм
из сплава Al-6%Mg.
Рис. 2 Схема экспериментальной установки[5].
В ходе работы наблюдалось закручивание оси цилиндрической рабочей
части образца в равномерную спираль, число витков которой совпадало с
числом оборотов. При достижении шага спирали значения, равного
диаметру, начинался изгиб рабочей части образца, сопровождающийся
уменьшением диаметра и увеличением длины стержня. Типичные кривые
осевой силы и крутящего момента представлены на рис. 3-4.
Рис. 3 Крутящий момент в зависимости от количества оборотов[6]
9
Рис. 4 Осевая сила в зависимости от количества оборотов[6]
Авторы статьи [1] высказывают
предположение, что сложная
взаимосвязь кручения и изгиба в нелинейной задаче Сен-Венана имеет
отношение к аналогичным эффектам в рамках задачи Кирхгоффа, хорошо
экспериментально и теоретически исследованным для упругих стержней
Таким образом, сформулируем задачу: имеется стержень с длиной
рабочей части l 112 мм, который закручивают с одного конца, вследствие
чего он начинает закручиваться и изгибаться, расстояние между концами
стержня не фиксировано, сам стержень жёстко закреплён. Известны
граничные условия для перемещений и углов поворота: перемещения
стержня на концах нулевые, углы Эйлера, определяющие поворот системы
координат, нулевые в начале, а нулевые два, а отвечающий за поворот
относительно оси стержня равен некоторой константе.
1.2.
Основные положения механики гибких стержней
При выводе уравнений равновесия считается, что поперечные
нормальные сечения стержня, плоские до деформации, остаются плоскими и
после деформации (гипотеза Бернулли), т. е. сдвиги не учитываются.
10
Размеры поперечного сечения остаются малыми по сравнению с длиной
стержня и радиусом кривизны оси стержня (под осью стержня понимается
линия, соединяющая центры тяжести площадей поперечных сечений
стержня). Различные, но статически эквивалентные локальные нагрузки
вызывают в стержне (если не учитывать местные напряжения вблизи точки
приложения нагрузок) одно и то же напряженное состояние (принцип СенВенана). Из гипотезы Сен-Венана следует, что продольные волокна (если
стержень представить состоящим из большого числа плотно прилегающих и
не связанных между собой волокон) не испытывают поперечного сжатия или
растяжения, а также касательных напряжений. Взаимные перемещения
сечений стержня при малых упругих деформациях в общем случае конечны,
т. е. задача является геометрически нелинейной, а физически — линейной
(перемещения точек осевой линии стержня могут быть большими, в то время
как материал стержня работает в пределах закона Гука).[1]
1.3.
Вывод нелинейной системы равновесия гибкого упругого стержня
для общего случая статической неопределимости
Стержень – это длинное тонкое тело, характеризующееся в первую
очередь своей осью – пространственной кривой. Ось определяется
зависимостью радиус-вектора от дуговой координаты. Далее будем
рассматривать стержень с постоянным сечением.
Для постановки задачи и последующего решения вводятся две системы
координат: неподвижная декартовая с единичными векторами
ai
и
подвижная, в терминах репера Френе pi рис. 5. Дуга s осевой линии
отсчитывается
от
некоторой
фиксированной
точки,
выбор
которой
произволен.
11
Рис. 5 Схема ориентации осей локального и
глобального базисов[7]
Записываем уравнения Кирхгоффа – аналог уравнений равновесия в
механике гибких стержней[8,9]:
dQ
0,
ds
d M p Q 0,
1
ds
(1)
где Q Q1 p1 Q2 p2 Q3 p3 – вектор внутренних усилий;
Q1 – осевое усилие;
Q2 , Q3 – перерезывающие усилия;
M M1 p1 M 2 p2 M 3 p3 – вектор внутренних моментов;
M 1 – крутящий момент;
M 2 , M 3 – изгибающие моменты;
p1 – единичный касательный вектор репера Френе.
12
p1 p2
p1 Q 1 0
Q1 Q2
p3
0 Q2 p3 Q3 p2 ,
Q3
Запишем соотношения для репера Френе – ортонормированного базиса,
присоединенного к кривой стержня[10]:
dp1
ds 1 p2 ,
dp2
2 p3 1 p1 ,
ds
dp3
ds 2 p2 ,
(2)
где p1 – единичный касательный вектор;
p2 – единичный вектор главной нормали;
p3 p1 p2 – единичный вектор бинормали к кривой в данной точке;
1 – кривизна стержня;
2 – кручение стержня.
Распишем производную вектора внутренних усилий:
dp
dp dQ
dp dQ
dQ dQ1
p1 Q1 1 2 p2 Q2 2 3 p3 Q3 3 ,
ds
ds
ds
ds
ds
ds
ds
dQ
dQ1
dQ
p1 Q1 1 p2 2 p2 Q2 2 p3 1 p1 3 p3 Q3 2 p2 0.
ds
ds
ds
Сгруппируем слагаемые относительно базисных векторов pi :
dQ3
dQ1
dQ2
Q
p
Q
Q
p
2Q2 p3 0.
1 2 1
1 1
2 3 2
ds
ds
ds
Равенство выполняется, если все выражения в скобках равняются
нулю, тогда:
13
dQ1
ds 1Q2 ,
dQ2
2Q3 1Q1 ,
ds
dQ3
ds 2Q2 .
(3)
Распишем производную вектора внутренних моментов:
dp dM 2
dp
dM dM1
p1 M1 1
p2 M 2 2
ds
ds
ds
ds
ds
dM 3
dp
p3 M 3 3 Q2 p3 Q3 p2 ,
ds
ds
Переносим всё в левую часть:
dM 3
dM1
dM 2
p1 M1 1 p2
p2 M 2 2 p3 1 p1
p3
ds
ds
ds
M 3 2 p2 Q2 p3 Q3 p2 0,
Сгруппируем слагаемые относительно базисных векторов pi :
dM1
dM 2
1M 2 p1
1M1 2 M 3 Q3 p2
ds
ds
dM 3
2 M 2 Q2 p3 0,
ds
Равенство выполняется, если все выражения в скобках равняются
нулю, тогда:
dM 1
ds 1M 2 ,
dM 2
2 M 3 1M 1 Q3 ,
ds
dM 3
ds 2 M 2 Q2 .
(4)
14
Будем считать, что упругие моменты пропорциональны изменениям
кривизны и кручения[8]:
M1 A11 2 o ,
M 2 0,
M A ,
33 1
3
(5)
где A11 , A33 – элементы матрицы жесткости A :
GJ
A 0
0
0
EJ
0
0
0 .
EJ
(6)
Подставляем (5) в (4):
d 2 o
A11
0 1 ,
ds
d0
A33 2 1 A11 2 o 1 Q3 ,
A22
ds
d 1
A33 ds 0 A22 2 Q2 .
(7)
Преобразуем (7):
d 2
do
A11
0,
A11
ds
ds
A33 A11
A11 o
1
1
Q3 0,
2 1
A22
A22
A22
d
1
Q2 .
1
A33
ds
Подставляем элементы матрицы жесткости (6):
15
d 2 d o
ds ds 0,
GJ o
EJ GJ
1
1
Q3 0,
2 1
EJ
EJ
EJ
d
1
Q2 .
1
ds
EJ
Перейдем к безразмерным величинам:
s d s,
d
, Q
EJ
Q.
d2
Модуль сдвига связан с модулем Юнга через коэффициент Пуассона:
GJ
EJ
,
1
где G – модуль сдвига,
E – модуль Юнга,
– коэффициент Пуассона,
J
d4
J
64
– момент инерции сечения стержня,
d4
32
– полярный момент инерции сечения стержня.
Подставим момент инерции и момент сопротивления в формулу связи
модуля сдвига с модулем Юнга
G
E
,
2 1
тогда, подставив полученное соотношение в систему, получим
16
d 2 d o
ds ds 0,
1
2 1
1
1
d 1
Q2 .
ds
o
1 Q3 0,
do
Из первого уравнения системы выразим
:
ds
d
do
2.
ds
ds
Из второго выражения системы выразим Q3 :
1
Q3
1
o
1
1
2 1.
Возьмем производную и подставим в выражение для
dQ3 1
ds 1
dQ3
:
ds
o
d 1 o d 2
d
d
1 1 2 2Q2 .
1
ds
ds
ds
1 ds
d 1 d o
Подставляем в полученное выражение
и
:
ds
ds
d 2
d 2
1
Q1 o
1 Q1 2 2Q2 .
1
ds
1
1 ds
Группируем относительно
d 2
ds 1
1
1
d 2
:
ds
1
o
Q
1
1
Q1 2 2Q2 .
1
1
Приведем к нормальной форме Коши:
17
d 2
ds
1
o
Q1
Q1 2
1
1
2Q2
1
После упрощения система уравнений статики имеет вид:
dQ1
ds 1Q2 ,
1 o
dQ2
Q
1 1
2
ds
1
1
d 1
ds Q2 ,
d 2 1 2Q2 Q1 o Q1 2
,
1 1
ds
do
1 2Q2 Q1 o Q1 2
.
ds
1
1
2
2 ,
(8)
где Q1 – осевое усилие;
Q2 , Q3 – перерезывающие усилия;
1 , 2 – кривизны стержня;
o – крутка стержня.
Движение стержня Коссера относительно его недеформированного
состояния описывается радиус-вектором и тензором поворота. Чтобы
получить
и
построить
в
дальнейшем
форму
стержня,
запишем
дифференциальное уравнение для радиус-вектора:
dr
p1.
ds
(9)
18
В свою очередь p1 связан с ортом базиса в недеформированном
0
состоянии стержня pi векторами соотношениями:
p l p0 ,
i ij j
pi0 l ji p j ,
где lij
(10)
– компоненты соответствующих ортогональных матриц L
1
и L ,
связывающих данные базисы.
0
dr
l1 j p j ,
ds
Матрица перехода L представляется в виде:
L L L L ,
(11)
где матрицы L , L и L записываются через углы Эйлера.
cos sin 0
L sin cos 0 ,
0
0
1
0
0
1
L 0 cos sin ,
0 sin cos
cos
L sin
0
sin
cos
0
(12)
0
0 .
1
Тогда, подставляя (12) в (11) и, пользуясь средствами Wolfram
Mathematica (прил. А), получим матрицу поворота локальной системы
координат:
19
cos cos cos sin sin cos cos sin cos sin sin sin
L cos sin cos cos sin cos cos cos sin sin cos sin ,
sin sin
cos sin
cos
cos cos cos sin sin cos sin cos cos sin sin sin
1
L cos cos sin cos sin cos cos cos sin sin cos sin .
sin sin
cos sin
cos
(13)
Распишем (9) в покомпонентной форме:
dr1
ds l11 ,
dr2
l12 ,
ds
dr3
l13 ,
ds
(14)
где l11 , l12 , l13 – элементы матрицы перехода L .
Тогда l11 cos cos cos sin sin ,
l12 cos cos sin cos sin ,
l31 sin sin .
dr1
ds cos cos cos sin sin ,
dr2
cos cos sin cos sin ,
ds
dr3
sin sin ,
ds
Для
окончательной
постановки
задачи
не
хватает
(15)
уравнений,
связывающих углы Эйлера с кривизнами. Получим эти соотношения из (16):
p ' LLT p,
(16)
где L ' – производная по s от матрицы L
LT – транспонированная матрица L, равная L1.
20
Результатом перемножения этих матриц будет кососимметричная
матрица (17).
0
'cos
'cos
0
'cos sin sin 'sin sin cos
Тогда,
подставляя
(17)
в
(16)
и
'cos sin sin
'sin sin cos (17)
0
расписывая
уравнения
в
покомпонентной форме получим (18):
dp1
ds 'cos p2 'cos sin sin p3 ,
dp2
(18)
'cos p1 'sin sin cos p3 ,
ds
dp3
ds 'cos sin sin p1 'sin sin cos p2 .
Подставим (2) в (20) и сгруппируем относительно pi . Получаем
систему (19)
'cos 1 p2 'cos sin sin p3 0,
(19)
'cos 1 p1 'sin sin cos 2 p3 0,
'cos sin sin p1 'sin sin cos 2 p2 0.
Очевидно, что равенства будут выполняться только в том случае, если
выражения в скобках будут равняться нулю. Отсюда получаем три
независимых уравнения связи углов Эйлера с кривизнами:
'cos 1 ,
'sin sin cos 2 ,
'cos sin sin 0.
(20)
Преобразуем систему так, чтобы в каждом уравнении осталось по
одной производной. В результате получаем недостающие три уравнения:
21
sin
d
,
2
ds
sin
d
2 cos ,
ds
d
ds 2 sin ctg 1.
(21)
Записываем полную систему уравнений:
dQ1
ds 1Q2 ,
1 o
2
dQ2
Q
1 1
2
2 ,
ds
1
1
d 1
ds Q2 ,
d 2 1 2Q2 Q1 o Q1 2
,
ds
1
1
do
1 2Q2 Q1 o Q1 2
,
ds
1 1
dr1
cos cos cos sin sin ,
ds
dr2
cos cos sin cos sin ,
ds
dr3
sin sin ,
ds
d
sin
,
2
ds
sin
d
ds 2 cos ,
d sin ctg .
2
1
ds
(22)
22
Глава 2 Исследование статических уравнений нелинейной
системы уравнений равновесия гибкого упругого стержня
2.1.
Общий случай решения статических уравнений системы
Обратимся к системе статических уравнений (8), представленной на
странице 16. Заменим в ней последнее уравнение уравнением связи 2 с
0
и перепишем систему статических уравнений в новом виде:
dQ1
ds 1Q2 , 23.1
1 o
dQ2
Q
1
1
2
ds
1
1
d 1
ds Q2 , 25.3
d 2 1 2Q2 Q1 o Q1 2
,
ds
1
1
o
d
2
0, 23.5
ds
2
2 ,
23.2
23.4
(23)
Проинтегрируем уравнение (23.5) и получим связь крутки с кручением
0 2 C1 ,
(24)
где C1 – некоторая константа интегрирования.
Подставим
в
уравнение
(23.1)
Q2
из
уравнения
(23.3)
и
проинтегрируем. В результате получим связь осевой силы Q1 с кривизной
стержня 1 :
1
Q1 12 C2 ,
2
(25)
где C2 – некоторая постоянная интегрирования.
23
Продифференцируем уравнение (23.3) и подставим в (23.2) наряду с
уравнениями (24) и (25), получим
d 2 1
1
1
1 12 C2
C
2
1
2
ds
1
1
2
2
2 ,
Раскрываем скобки:
d 2 1 1 3
C
1 C2 1 1 1 2 22 1. (26)
ds
2
1
Рассмотрим уравнение (23.4), подставим в него выражения (24), (25) и
(23.3)
d 2
ds
1 2
d 1 1 2
1
1 C2 2 C1 12 C2
ds 2
2
,
1 1
или же
1
d
1
1 1 2 12 2 C2 2 C1 12 C2
d 2
2
2
. (27)
ds
ds
1 1
Система (23) свелась к системе двух уравнений (уравнения (26) и (27))
с двумя неизвестными функциями и двумя неизвестными константами.
Можно оставить в таком виде и попробовать решить, однако удобнее
привести её к нормальной форме Коши:
dQ2 1 Q1 1 o 2 22 ,
ds
1
1
d 1
Q2 ,
(28)
ds
2 1 2Q2 12 C2 1 2 C1
d
2
,
ds
2 1 1
Система (28) описывает распределение силовых характеристик вдоль
оси стержня для общего случая постановки задачи. В дальнейшем
24
рассмотрим частные случаи решения этой системы дифференциальных
уравнений, полагая, что одна из её искомых функций является константой.
2.2.
Частные решения системы для случая постоянной крутки стержня
Предположим, что
o – константа. Рассмотрим два случая – когда
крутка стержня тождественно равна нулю, и когда она константа, неравная
нулю.
1)
Итак, пусть o 0. Из условия равенства производных получаем, что
2 Const. Подставляем это значение в систему статических уравнений –
одно дифференциальное уравнение превращается в алгебраическое, остается
четыре дифференциальных уравнения
dQ1
ds 1Q2 ,
dQ2
Q
1
1
ds
1
d 1 Q ,
2
ds
d 2 0,
ds
2
2 ,
и дополнительное равенство
1 2Q2 Q1 2
1 1
0.
Из последнего уравнения следует, что 2 0 или Q2
1
Q1.
a) Допустим, что 2 0. Тогда остается система трех уравнений с тремя
неизвестными:
25
dQ1
ds 1Q2 ,
dQ2
1Q1 ,
ds
d 1
ds Q2 .
b) Предположим, что
Q2
1
Q1. Тогда остается система трех
уравнений с тремя неизвестными:
dQ1
Q1 ,
1
ds
1
d 1
Q1 ,
ds
1
d 2
ds 0.
и дополнительное условие на Q1 :
1
dQ1
1
Q1 Const .
ds
2)
Рассмотрим
случай,
когда
константа
o 0. Тогда система
статических уравнений преобразуется к следующему виду:
26
dQ1
ds 1Q2 ,
1 o
dQ2
Q
1 1
2
ds
1
1
d 1
ds Q2 ,
d 2 1 2Q2 Q1 o Q1 2
,
ds
1
1
1 Q Q o Q
2 2
1
1 2
0.
1
1
2
2 ,
Рассмотрим алгебраическое уравнение:
1 2Q2 Q1 o Q1 2
1 1
Равенство
нулю
выполняется
1 2Q2 Q1 o Q12 0.
в
0
двух
случае,
если
Данное выражение содержится как в
алгебраическом уравнении, так и в уравнении на производную на 2 . Из
этого получаем следующее уравнение:
d 2
0,
ds
значит, 2 – константа.
Выразим
Q2
из
условия
равенства
нулю
выражения
1 2Q2 Q1 o Q12 0.
1 2Q2 Q1 o Q12 0.
1 2Q2 Q1 o 2 ,
Q2
o
2
1 2
Q1 ,
27
где дробь при Q1 – некоторая константа. Для упрощения выражений
обозначим её . Тогда получим следующее соотношение:
o
2
1 2
,
Q2 Q1.
Так как 2 оказалась константой, система статических уравнений
уменьшилась:
dQ1
ds 1Q2 ,
dQ2
1
1 Q1
1
ds
d 1
Q2 ,
ds
где
o
2
1
2
2 ,
o , 2 – некоторые константы. Тогда выражение в скобках при 1 в
дифференциальном уравнении относительно Q2 есть сумма Q1 и некоторой
константы, обозначим её . Тогда
dQ1
ds 1Q2 ,
dQ2
1 Q1 ,
ds
d 1
ds Q2 ,
2 o
где
2 .
1
Рассмотрим получившуюся
систему
уравнений
с учетом всех
полученных условий:
28
dQ1
ds 1Q2 ,
dQ2
1 Q1 ,
ds
d 1
ds Q2 ,
где
Q2 Q1 ,
o
2
1 2
,
2 o
2 .
1
Подставим
соотношение для Q2 :
dQ1
ds 1Q1 ,,
dQ1
1 Q1 ,
ds
d 1
ds Q1 ,,
Получили два дифференциальных уравнения относительно Q1.
a) Рассмотрим случай, когда 0. Из этого следует следующее
равенство: o 2 . Тогда Q2 0 Q1 , т.е. Q2 0.
Из условия 0 система преобразуется к следующему виду:
dQ1
ds 0,
1Q1 0,
d
1 0,
ds
Из этой системы получаем, что 1 const , Q1 const , при этом
1Q1 0. Возможны два случая: 1 0, тогда Q1 const , либо Q1 0,
тогда 1 const. Однако, условие 1 0 не допускает наша система
уравнений, значит, остается случай Q1 0, 1 const.
29
b) Рассмотрим случай, когда: , 0.
dQ1
ds 1Q1 ,,
dQ1
1 Q1 ,
ds
d 1
ds Q1 ,,
Получили два уравнения относительно производной от Q1 , то есть,
получаем следующее равенство:
1Q1
1
Q ,
1
1 2
Q1 Q1 0.
Получаем два случая:
1 0,
2
Q1 1 0.
I.
Если 1 0, то получаем систему:
dQ1
0,
ds
Q1 0,
тогда Q1 0. Отсюда Q2 тоже равняется нулю. Получаем следующее:
Q1 0,
Q2 0,
1 0,
C const
1
2
o C2 const
30
Если Q1
II.
2
1 0, то Q1
2 1
Const , но тогда 1
линейно зависит от s :
1
Из
этого
следует,
что
2
1
s.
dQ1
2
sQ1.
ds
1
Приходим
к
противоречию.
2.3.
Частный случай эволюции формы стержня в виде однородной
цилиндрической спирали
Сформулируем вспомогательную задачу о винтовой линии, вид
которой описывается следующим уравнением:
x R(cos
ks
1 R k
2
2
a1 sin
ks
1 R k
2
2
a2 )
s
1 R k
2
2
a3 ,
(20)
где s – длина дуги вдоль спирали, т.е. длина стержня,
k – параметр спирали,
R – радиус спирали.
Запишем геометрические соотношения для спиральной линии:
ks
ks
1
a1 cos
a2
,
sin
2 2
2 2
2 2
1
R
k
1
R
k
1
R
k
Первая производная от p1 имеет вид
dx
Rk
p1
ds
1 R2k 2
dp1
Rk 2
ks
cos
2 2
ds
1 R 2 k 2
1 R k
ks
a
sin
1
2 2
1 R k
a2 ,
Известна связь p1 с p2 :
dp1
1 p2 ,
ds
31
отсюда получаем 1 и p2 :
Rk 2
1
,
1 R2k 2
ks
p2 cos
2 2
1 R k
ks
a
sin
1
2 2
1 R k
a2 .
Первая производная от p2 имеет вид
dp2
k
ds
1 R2k 2
ks
ks
sin
a1 cos
a2 .
2 2
2 2
1
R
k
1
R
k
Известна связь p2 с p3 через p1 :
dp2
2 p3 1 p1 ,
ds
тогда 2 и p3 получим из выражения:
ks
sin
1 R 2 k 2 1 R 2 k 2
k
2 p3
R2k 3
1 R k
3
2
2 p3
2
3
ks
sin
2 2
1 R k
k
1 R2k 2
3
ks
a1 cos
2 2
1 R k
ks
a1 cos
2 2
1 R k
ks
sin
2 2
1 R k
1 R k
a2
ks
a
cos
1
2 2
1 R k
Rk 2
2
a2
2 3
Rk 2
1 R k
2
2 3
a2
a3 ,
тогда получим 2 и p3 :
2
k
,
1 R2k 2
32
a3 .
p3
ks
ks
Rk
sin
a1 cos
a2
a,
2 2
2 2 3
1 R 2 k 2 1 R 2 k 2
1
R
k
1
R
k
1
Получим статические характеристики системы (25), предварительно
задав геометрические параметры спирали: длина дуги l 20, радиус
спирали R 1, параметр спирали k 1. Подставим значения 1 0.5 и
2 0.5 в (30). Сразу же получаем из одного из уравнений значение
перерезывающего усилия Q2 0. Далее, используя средства прикладного
пакета Wolfram Mathematica, определяем из оставшейся системы двух
уравнений
1 2
2
C
C
1 22 ,
2
1 1
2
1 12 C2 1 2 C1 0,
2
константы C1 0.65 и C2 0.125 для заданных условий. Далее, подставляя
полученные постоянные в (24) и (25), находим значения крутки и осевой
силы. В рамках заданных условий значения получили следующие величины
Q1 0, а 0 0.15, что совпадает со значением 0 2 . Таким образом
получаем
решение
системы,
совпадающее
с
частным
случаем
2а
предыдущего параграфа.
Проверим, удовлетворяют ли кинематические уравнения системы (22)
кинематике спирали. Для этого подставим 1 и 2 в (21) и, используя
средства MATLAB, численно решим полученную систему уравнений
совместно с уравнениями (15).
Для решения задачи зададим следующие входные условия: Условия на
левом конце стержня:
r1 0 r2 0 r3 0 0, 0
3
, 0
2
, 0
4
.
33
Численное решение получим по методу Рунге-Кутты 4-го порядка, для
этого используем процедуру ode45, листинг программы представлен в
приложении 2. Построим полученные функции изменения радиус-вектора и
углов поворота вдоль оси стержня, графические зависимости представлены
на рис. 6-7.
Рис. 6 Изменение радиус-вектора вдоль оси стержня
Рис. 7 Изменения углов поворота вдоль оси стержня.
34
Форма осевой линии стержня представлена на рис. 8
Рис. 8 Эволюция формы стержня при заданных условиях
Таким образом, получено решение поставленной задачи и убедились,
что выведенные кинематические уравнения (15) и (21) совместны с
кинематикой спирали.
35
Глава 3 Численное решение системы уравнений в среде
MATLAB
3.1.
Использование численных методов для решения нелинейных
систем дифференциальных уравнений
В предыдущей главе были рассмотрены несколько частных случаев
решения нелинейной системы уравнений равновесия гибкого упругого
стержня. Непосредственное решение аналитическими методами затруднено
ввиду сложности системы (24), необходимо использовать численные методы
решения систем дифференциальных уравнений. Вычисления буду вестись в
пакете прикладных программ MATLAB, который представляет собой
высокоуровневый язык и интерактивную среду для программирования,
численных расчетов и визуализации результатов.
Далее будем использовать два численных метода: метод пристрелки,
позволяющий свести краевую задачу к задаче Коши, и метод Рунге-Кутты 4го порядка, использующийся для решения начальных задач. В большинстве
своем задачи, связанные со стержнями, – краевые. Чтобы решить такую
задачу, необходимо свести её к начальной, то есть получить из граничных
условий, заданных для части функций, входящих в систему уравнений,
начальные условия на все искомые функции. В нашем случае – перейти от
условий на правой границе, заданных для кинематической составляющей
системы уравнений равновесия, к начальным условиям на уравнения статики.
Для этого напишем программу в среде MATLAB, способную сводить
краевую задачу к задаче Коши и решающую её.
3.2.
Численное решение нелинейной системы уравнений равновесия
гибкого упругого стержня в случае постановки начальной задачи
36
Проверим корректность работы программы на начальной задаче,
решение которой мы точно знаем. Такой задачей является частный случай
эволюции формы стержня в виде однородной цилиндрической спирали.
Воспользуемся статическими начальными относительно осевой линии
стержня условиями, полученными в параграфе 2.3 предыдущей главы и
выбранными там же кинематическими и геометрическими.
Сформулируем задачу следующим образом: имеется нелинейная
система дифференциальных уравнений (22), физический параметр –
коэффициент Пуассона
0.3, геометрический параметр – длина осевой
линии l 20, и одиннадцать начальных условий, представленных ниже:
Q1 0, Q2 0, 1 0.5, 2 0.5, 0 0.15,
r1 0 r2 0 r3 0 0, 0
3
, 0
2
, 0
4
.
Используя метод Рунге-Кутты 4-го порядка, записанный в процедуру
ode45, решим поставленную начальную задачу. Отобразим полученное
решение на рис. 9, построив на одном графике распределение силовых
факторов вдоль осевой линии стержня, а на рис. 10 – изменение
кинематических характеристик стержня.
Рис. 9 Распределение силовых характеристик по длине стержня
37
Рис. 10 Изменение координат радиус вектора и углов
поворота вдоль оси стержня
Форма стержня представлена на рис. 11
Рис. 11 Форма осевой лини стержня при заданных условиях
38
Решение совпадает с решением из предыдущей главы, полученная
форма
равновесия
стержня
–
однородная
цилиндрическая
спираль.
Программа для решения начальных задач работает корректно.
Решим новую задачу, полагая, что 1 и 2 не являются константами.
Начальные условия возьмем из предыдущей задачи, заменив только значения
1 и Q1 на левом конце стержня; Q1 0.375,
1 1. Решаем
поставленную задачу, результаты вычислений представлены на рис. 12-13.
Рис. 12 Распределение силовых характеристик по длине стержня
при начальных условиях Q1 0.375, 1 1.
39
Рис. 13 Изменение кинематических характеристик по длине
стержня при начальных условиях Q1 0.375, 1 1.
Кривая осевой линии стержня представлена на рис.14.
Рис. 14 Форма равновесия стержня при начальных условиях
Q1 0.375, 1 1.
Анализируя
графики,
представленные
на
рис. 12,
видим,
что
полученные силовые характеристики есть функции, изменяющиеся вдоль
40
осевой стержня, получена форма равновесия стержня для поставленных
начальных условий. Программа для решения начальных задач работает
корректно, переходим к краевым задачам.
3.3.
Численное решение нелинейной системы уравнений равновесия
гибкого упругого стержня в случае постановки краевой задачи
Для того, чтобы перейти от краевой задачи к начальной, будем
использовать метод пристрелки. Будем считать, что нам известны
перемещения на правом конце стержня, а на левом – осевая и
перерезывающая силы, кривизна стержня, перемещения и углы поворота
локального базиса. Тогда в качестве неизвестных выступают значения
кручения и крутки на левом конце стержня.
Для решения поставленной задачи была написана программа в среде
MATLAB, реализующая с помощью метода пристрелки переход от краевой
задачи к задаче Коши
Проверим корректность работы программы на краевой задаче, решение
которой мы точно знаем. Такой задачей является частный случай эволюции
формы стержня в виде однородной цилиндрической спирали. Частично
используем статическими начальными относительно осевой линии стержня
условиями, полученными в параграфе 2.3 предыдущей главы. Два начальных
условия оставим неизвестными: на 2 и на
0 – их подберем с помощью
метода пристрелки. Кинематические краевые условия выберем таким
образом, чтобы спираль была сориентирована вдоль одной из неподвижных
осей.
Сформулируем задачу следующим образом: имеется нелинейная
система дифференциальных уравнений (22), физический параметр –
коэффициент Пуассона
0.3,
геометрический параметр – длина осевой
линии l 17, и одиннадцать смешанных условий, представленных ниже:
Q1 0, Q2 0, 1 0.5, r2 0 0,
41
r1 0 r1 17 0, r3 0 r3 17 0, 0
4
, 0
Методом пристрелки подбираем значения 2 и
2
, 0 0.
0 на левом конце
стержня. На рис. 15 представлена функция невязки
Рис. 15 Сходимость краевой задачи к задаче Коши при подборе
начальных условий a и b
Задача сводится к начальной при следующих значениях a и b :
a 0.5 2 0 0.5
b 0.15 0 0 0.15
Решаем полученную начальную задачу. Распределение силовых и
кинематических параметров представлено на рис. 16-17.
42
Рис. 16 Распределение силовых характеристик по длине стержня
0
при начальных условиях 2 0.5, 0.15.
Рис. 17 Изменение кинематических характеристик по длине
0
стержня при начальных условиях 2 0.5, 0.15.
Получили решение в конфигурации однородной цилиндрической
спирали, геометрия осевой Коссера представлена на рис.18:
43
Рис. 18 Форма равновесия стержня тестовой краевой задачи
Таким образом на тестовой задаче убедились, что алгоритм для
перехода от краевой задачи к начальной работает корректно.
Решим задачу для новых смешанных условий: пусть
1 0.4,
остальные условия возьмем из предыдущей задачи, в роли подбираемых
неизвестных условий оставим 2 и 0 . По методу пристрелки перейдем к
начальной задаче, функция невязки представлена на рис. 19.
Рис. 19 Сходимость краевой задачи к задаче Коши при подборе
начальных условий a и b
44
Задача сводится к начальной при следующих значениях a и b :
a 0.27 2 0 0.27,
b 0.54 0 0 0.54.
Решаем полученную начальную задачу. Распределение силовых и
кинематических параметров представлено на рис. 20-21.
Рис. 20 Распределение силовых характеристик по длине стержня
0
при начальных условиях 2 0.27, 0.54.
45
Рис. 21 Изменение кинематических характеристик по длине
стержня при начальных условиях 2 0.27, 0 0.54.
Построим конфигурацию стержня поставленной задачи, вид кривой
представлен на рис.22.
Рис. 22 Форма равновесия стержня поставленной краевой задачи
46
Таким образом, получено решение краевой задачи путем сведения её к
начальной методом пристрелки. В дальнейшем, написанный алгоритм можно
использовать для решения начальных и краевых задач, а также для задач со
смешанным типом граничных условий.
47
Заключение
В работе исследуется эволюция формы гибкого упругого стержня при
его квазистатическом кручении. Выводится нелинейная система уравнений
равновесия гибкого упругого стержня для общего случая статической
неопределимости, которая необходима для объяснения нелинейных эффектов,
обнаруженных экспериментально и описанных в статье [6]. Аналитические
преобразования проводятся с использованием прикладного пакета программ
Wolfram Mathematica. Рассматриваются возможные специализации этой
системы, соответствующие определенным гипотезам о виде распределения
крутящего момента по длине стержня. Реализуется численное решение
краевой задачи о кручении стержня, жестко закрепленного с обоих концов, в
программе MATLAB.
Выводы по работе:
1.
Выведена нелинейная система уравнений равновесия гибкого стержня
при его кручении.
2.
Получено
точное
частное
решение
системы,
соответствующее
конфигурации стержня в виде однородной цилиндрической спирали.
3.
В пакете MATLAB написана программа численного решения системы
дифференциальных уравнений методом Рунге-Кутты четвертого порядка.
При этом решение краевой задачи сводилось к решению задачи Коши
методом
пристрелки.
Корректная
работа
программы
подтверждена
воспроизведением решения в виде однородной цилиндрической спирали.
4.
Получены и исследованы решения ряда начальных и краевых задач.
Численная процедура проведена с путем подбора недостающих начальных
условий.
5.
В результате создан инструмент для исследования нелинейных
эффектов при кручении гибкого стержня.
48
Список литературы
1.
Светлицкий
В.А.
Механика
гибких
стержней
и
нитей.
—
М.:Машиностроение, 1978. — 222 с.
2.
Kehrbaum S., Maddocks J. H., Elastic Rods, Rigid Bodies, Quaternions and
the Last Quadrature. // Philosophical Transactions: Mathematical, Physical and
Engineering Sciences – 1997 –Vol. 355 – No. 1732 – P. 2117–2136.
3.
Van der Heijden G.H.M., Thompson J.M.T. Helical and Localised Buckling
in Twisted Rods: A Unified Analysis of the Symmetric Case. // Nonlinear
Dynamics. – 2000. – Vol.21. – P. 71–99.
4.
Choualeb N., Maddocks J.H. Kirchhoff’s Problem of Helical Equilibria of
Uniform Rods. // Journal of Elasticity. – 2004. – Vol.77. – P. 221–247.
5.
Goss V.G.A., van der Heijden G.H.M., Thompson J.M.T., Neukirch S.
Experiments on Snap Buckling, Hysteresis and Loop Formation in Twisted Rods.
// Experimental Mechanics. – 2005. – Vol.45, – No.2. – P.101-111.
6.
Bachurikhin V.P., Keller I.E., Merzlyakov A.F., Yurlov M.A. Experimental
study of nonlinear effects under torsion of the uniform cylinder with initially
circular cross section. // Solid State Phenomena. – 2016 – Vol.243. – P.29-34.
7.
Елисеев В. В. Механика упругих тел. СПб.: Изд-во СПбГТУ, 1999. -
341 с.
8.
Светлицкий В. А. Статистическая механика и теория надежности:
учебник для втузов / В. А. Светлицкий. – М.: Изд-во МГТУ им. Н. Э.
Баумана, 2002. – 504 с.
9.
Светлицкий В. А. Строительная механика машин. Механика стержней:
учебник для вузов : в 2 т. / В. А. Светлицкий. - Москва: Физматлит, 2009.
10.
Келлер И.Э. Тензорное исчисление: Учебное пособие. – СПб.:
Издательство «Лань». 2012. – 176 с.
49
Приложение А
F= {{Cos[θ[s]], -Sin[θ[s]],0},{Sin[θ[s]], Cos[θ[s]],0}, {0, 0, 1}}
A ={{Cos[ϕ[s]], -Sin[ϕ[s]], 0},{Sin[ϕ[s]], Cos[ϕ[s]], 0},{0, 0, 1}}
B = {{1, 0, 0},{0, Cos[ψ[s]], -Sin[ψ[s]]},{0, Sin[ψ[s]], Cos[ψ[s]]}}
A//MatrixForm
B//MatrixForm
F//MatrixForm
L =F.B.A
L//MatrixForm
Lt = Transpose [L]
Lt//MatrixForm
Ld = D[L,s]//FullSimplify
Ld//MatrixForm
LL = Ld.Lt//FullSimplify
LL//MatrixForm
p= {{p1},{p2},{p3}}
p//MatrixForm
P=LL.p
P//MatrixForm
G = {{K1*p2},{K2*p3-K1*p1},{-K2*p2}}
PP=P-G//FullSimplify
PP//MatrixForm
SYSR = {{-Cos[θ[s]] *Sin[ψ[s]],Sin[θ[s]] ,0 ,0},{Cos[ψ[s]],0,1,K1},
{Cos[ψ[s]],0,1,K1},{-Sin[θ[s]] *Sin[ψ[s]],-Cos[θ[s]],0,-K2},
{Cos[θ[s]]*Sin[ψ[s]],-Sin[θ[s]],0,0},{Sin[θ[s]]* Sin[ψ[s]],Cos[θ[s]],0,K2}}
SYSR //MatrixForm
MatrixRank[SYSR]
SYS = {{-Cos[θ[s]] *Sin[ψ[s]],Sin[θ[s]] ,0 },{Cos[ψ[s]],0,1},{Cos[ψ[s]],0,1},{Sin[θ[s]] *Sin[ψ[s]],-Cos[θ[s]],0},{Cos[θ[s]]*Sin[ψ[s]],-Sin[θ[s]],0},{Sin[θ[s]]*
Sin[ψ[s]],Cos[θ[s]],0}}
SYS //MatrixForm
MatrixRank[SYS]
50
Приложение Б
function dy = su4(t,y)
R = 1;
k = 1;
k1 = R*k^2/(1+R^2*k^2);
k2 = k/(1+R^2*k^2);
dy=zeros(6,1);
dy(1) = cos(y(6)).*cos(y(4)) - cos(y(5)).*sin(y(6)).*sin(y(4));
dy(2) = - cos(y(4)).*cos(y(5)).*sin(y(6)) - cos(y(6)).*sin(y(4));
dy(3) = sin(y(6)).*sin(y(5));
dy(4) = -k2.*sin(y(6))./sin(y(5));
dy(5) = -k2.*cos(y(6));
dy(6) = k2.*sin(y(6)).*cot(y(5)) - k1;
end
clc
clear all
close all
C=0;
I=0;
[T,Y] = ode45(@su4, [0 20], [0 0 0 pi/3 pi/2 pi/4]);
[C,I] = max(T);
figure
plot(T,Y(:,1),T,Y(:,2),T,Y(:,3),'LineWidth', 1.5)
legend({'u_{1}', 'u_{2}', 'u_{3}'},'FontSize',14,'Location','bestoutside')
xlabel('\its','FontSize',14)
ylabel('\itu_i','FontSize',14)
grid on
figure
plot(T,Y(:,4),T,Y(:,5),T,Y(:,6),'LineWidth', 2)
legend({'\phi', '\psi', '\theta'},'FontSize',14,'Location','bestoutside')
xlabel('\its','FontSize',14)
ylabel('\it\phi_i','FontSize',14)
grid on
figure
plot3(Y(:,1),Y(:,2),Y(:,3),'LineWidth', 2.5)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('z','FontSize',14)
grid on
daspect([1,1,1])
% view([0,1,0])
51
Приложение В
function dy = su3(t,y)
nu = 0.3;
C1 = nu/(1+nu);
C2 = 1/(1+nu);
C3 = 1+nu;
dy=zeros(11,1);
dy(1) = y(3).*y(2);
dy(2) = -y(3).*(y(1) - C2.*y(5).*y(4) + C1.*y(4).*y(4));
dy(3) = -y(2);
dy(4) = (C3.*y(4).*y(2) - y(1).*y(5) + nu.*y(1).*y(4))./(C3.*y(3));
dy(5) = -(C3.*y(4).*y(2) - y(1).*y(5) + nu.*y(1).*y(4))./(C3.*y(3));
dy(6) = cos(y(11)).*cos(y(9)) - cos(y(10)).*sin(y(11)).*sin(y(9));
dy(7) = - cos(y(9)).*cos(y(10)).*sin(y(11)) - cos(y(11)).*sin(y(9));
dy(8) = sin(y(11)).*sin(y(10));
dy(9) = -y(4).*sin(y(11))./sin(y(10));
dy(10) = -y(4).*cos(y(11));
dy(11) = y(4).*sin(y(11)).*cot(y(10)) - y(3);
clc
clear all
close all
R = 1;
k = 1;
nu = 0.3;
k1 = R*k^2/(1+R^2*k^2);
k2 = k/(1+R^2*k^2);
k0 = - k2 + 0.65;
q1 = -1/2*k1^2 + 0.125;
q2 = 0;
C=0;
I=0;
[T,Y] = ode45(@su3, [0 17.7572], [q1 q2 k1 k2 k0 0 0 0 pi/4 pi/2 0]);
[C,I] = max(T);
disp([Y(end,7),Y(end,8)])
figure
plot(T,Y(:,1),'k', T,Y(:,2),'g--', T,Y(:,3),'y', T,Y(:,4),'r--',T,Y(:,5),
'LineWidth',1.5)
legend({'Q_{1}', 'Q_{2}', '\kappa_{1}', '\kappa_{2}', '\kappa^{o}'},...
'FontSize',14,'Location','bestoutside')
xlabel('\its')
ylabel('\itQ_i, \it\kappa_{i}')
grid on
figure
plot(T,Y(:,6),'b', T,Y(:,7), T,Y(:,8),'c', T,Y(:,9),
T,Y(:,10),T,Y(:,11),'y','LineWidth', 1.5)
legend({'u_{1}', 'u_{2}', 'u_{3}','\phi', '\psi',
'\theta'},'FontSize',14,'Location','bestoutside')
xlabel('\its')
ylabel('\itu_i, \it\phi_i')
grid on
figure
plot3(Y(:,6),Y(:,7),Y(:,8),'LineWidth', 2.5)
daspect([1,1,1])
xlabel('x')
ylabel('y')
zlabel('z')
grid on
52
Приложение Г
function dy = su3(t,y)
nu = 0.3;
C1 = nu/(1+nu);
C2 = 1/(1+nu);
C3 = 1+nu;
dy=zeros(11,1);
dy(1) = y(3).*y(2);
dy(2) = -y(3).*(y(1) - C2.*y(5).*y(4) + C1.*y(4).*y(4));
dy(3) = -y(2);
dy(4) = (C3.*y(4).*y(2) - y(1).*y(5) + nu.*y(1).*y(4))./(C3.*y(3));
dy(5) = -(C3.*y(4).*y(2) - y(1).*y(5) + nu.*y(1).*y(4))./(C3.*y(3));
dy(6) = cos(y(11)).*cos(y(9)) - cos(y(10)).*sin(y(11)).*sin(y(9));
dy(7) = - cos(y(9)).*cos(y(10)).*sin(y(11)) - cos(y(11)).*sin(y(9));
dy(8) = sin(y(11)).*sin(y(10));
dy(9) = -y(4).*sin(y(11))./sin(y(10));
dy(10) = -y(4).*cos(y(11));
dy(11) = y(4).*sin(y(11)).*cot(y(10)) - y(3);
end
clc
clear all
close all
R = 1;
k = 1;
nu = 0.3;
% k1 = R*k^2/(1+R^2*k^2);
k1 = 0.4;
k2 = k/(1+R^2*k^2);
k0 = - k2 + 0.65;
% q1 = -1/2*k1^2 + 0.125;
q1 = 0;
q2 = 0;
a = -0.5;
z = 1;
Z = zeros(1,5);
for i = 1:100
b = -0.15;
for j = 1:100
C=0;
I=0;
[T,Y] = ode45(@su3, [0 10], [q1 q2 k1 a b 0 0 0 pi/4 pi/2 0]);
[C,I] = max(T);
f=abs(Y(I,7))+abs(Y(I,8));
A(i)=a;
B(j)=b;
F(i,j)=f;
if f < z
z = f;
Z = [a,b];
disp(z)
disp(Z)
end
b = b + 0.01;
end
53
a = a + 0.01;
end
figure
surf(B,A,F)
xlabel('b')
ylabel('a')
zlabel('f')
C = 0;
I = 0;
[T,Y] = ode45(@su3, [0 10], [q1 q2 k1 Z(1,1) Z(1,2) 0 0 0 pi/4 pi/2 pi]);
[C,I] = max(T);
figure
plot(T,Y(:,1),'k', T,Y(:,2),'g--', T,Y(:,3),'y', T,Y(:,4),'r--',T,Y(:,5),
'LineWidth',1.5)
legend({'Q_{1}', 'Q_{2}', '\kappa_{1}', '\kappa_{2}', '\kappa^{o}'},...
'FontSize',14,'Location','bestoutside')
xlabel('\its')
ylabel('\itQ_i, \it\kappa_{i}')
grid on
figure
plot(T,Y(:,6),'b', T,Y(:,7), T,Y(:,8),'c', T,Y(:,9),
T,Y(:,10),T,Y(:,11),'y','LineWidth', 1.5)
legend({'u_{1}', 'u_{2}', 'u_{3}','\phi', '\psi',
'\theta'},'FontSize',14,'Location','bestoutside')
xlabel('\its')
ylabel('\itu_i, \it\phi_i')
grid on
figure
plot3(Y(:,6),Y(:,7),Y(:,8),'LineWidth', 2.5)
daspect([1,1,1])
xlabel('x')
ylabel('y')
zlabel('z')
grid on
54
Отзывы:
Авторизуйтесь, чтобы оставить отзыв