Петлевой маятник - численное моделирование. ТЮФ Looping pendulum.

Петлевой маятник - численное моделирование. ТЮФ Looping pendulum.

В этот раз на примере петлевого маятника - интересной и красивой механической задачи из предстоящего в 2018 году Турнира Юных физиков я покажу, как можно решать сложные, не имеющие аналитического решения задачи с помощью компьютера. У меня уже есть один пост, посвященный задачке Турнира Юных Физиков. Сейчас всё будет интереснее, ведь мы не просто получим решение, но ещё и смоделируем его на компьютере!
Думаю, что и этот пост будет не последним по этой теме и разборы некоторых механических задач будут традиционно появляться здесь.

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

Connect two loads, one heavy and one light, with a string
over a horizontal rod and lift up the heavy load by pulling
down the light one. Release the light load and it will sweep
around the rod, keeping the heavy load from falling to the
ground. Investigate this phenomenon.


Описание движения системы


Рассмотрим динамику данной системы:





На маленький груз, так-же как и на большой, действует сила натяжения и сила тяжести. При этом, естественно, силы натяжения на грузы различны, так как между верёвкой и блоком есть трение. Найти связь между \(T\) и \(T_1\) достаточно просто. Рассмотрим небольшой участок стержня:




При скольжении сила трения на небольшой кусочек нити будет \(F_{тр} = \mu N\). Сила реакции опоры

$$
N = Tsin(d\theta/2) + (T+dT)sin(d\theta/2)
$$
При малых углах \(sin(\alpha) \approx \alpha\). То есть

$$
N = Td\theta + dTd\theta/2
$$
Вторым членом \(dTd\theta/2\) мы можем пренебречь, так как он второго порядка малости.
Считаем нить невесомой, так что сумма сил, действующих на рассматриваемый нами кусочек, равна нулю:

$$
T + \mu Td\theta = T+dT
$$
Отсюда приходим к выражению

$$
\frac{dT}{T} = \mu d\theta
$$
при интегрировании которого получаем формулу Эйлера:



$$
T_1 = T e^{\mu\theta} \qquad (1)
$$


Отлично, ведь одна из подзадач решена! С помощью результата (1) можно качественно объяснить, почему маленький груз так просто останавливает большой. Так же нужно объяснить, почему он будет вращаться по спирали. Это достаточно просто и станет понятно немного позже.
Теперь рассмотрим динамику маленького и большого грузов более подробно:




Для большого груза запишем второй закон Ньютона:

$$
\left\{
\begin{array}{ccc}
MA_x = 0\\
MA_y = Te^{\mu\theta} - Mg\\
\end{array}
\right. \qquad (2)
$$
Найти силу натяжения \(T\) достаточно просто: нужно только заметить, что маленький груз вращается по окружности относительно точки касания нити со стержнем. Это значит, что в проекции на саму нить сумма сил должна обеспечивать центростремительное ускорение \(a_ц = v_{t}^2/s\):

$$
T + mg cos\theta = m\frac{v_{t}^2}{s}
$$
Отсюда
$$
T = m\frac{v_{t}^2}{s} - mg cos\theta \qquad (3)
$$
Где \(v_{t} = \frac{d\theta}{dt}s\)
Закон изменения угла \(\theta\) можно найти, записав 2 закон Ньютона для вращательного движения:

$$
I\frac{d\left(\frac{d\theta}{dt}\right)}{dt} = \sum{[r\times F]}
$$
Где \(I\) - момент инерции и в нашем случае \(I = ms^2\), а сумма из правой части - это общий момент сил относительно выделенной точки.
В итоге, получаем следующую формулу (вращение относительно точки контакта нити со стержнем):

$$
\frac{d\left(\frac{d\theta}{dt}\right)}{dt} = \frac{gsin\theta}{s} \qquad (4)
$$
Как изменяются координаты большого груза и угол \(\theta\) мы поняли. Осталось только найти закон изменения \(s\), который получается из условия нерастяжимости нити. Допустим, что длинна нити равна \(L\), тогда \(L = s + R\theta + H\), где \(H\) - это расстояние от большого груза до точки касания. В дифференциальном виде:

$$
ds = -dH - Rd\theta \qquad (5)
$$
Мы получили, что изменение угла и \(H\) приводит к изменению расстояния от маленького грузика до точки контакта нитки со стержнем. В свою очередь это приводит к увеличению угловой скорости из-за закона сохранения момента импульса, который в нашем случае равен \(I\omega\). То есть при уменьшении \(s\) угловая скорость малого груза будет сильно расти (\(\omega \sim 1/s^2\)). Это и приводит к закручиванию груза.
В принципе это всё, что нужно для численного моделирования.

\(dH\) мы знаем из (2), а \(d\theta\) находим численным интегрированием (4). И об этом я не поленюсь рассказать немного поподробнее.

Математическое отступление


Допустим, что угловая скорость \(\frac{d\theta}{dt} = \omega\), а угловое ускорение \(\frac{d\omega}{dt} = \varepsilon\).
Для угла \(\theta\) мы имеем уравнение второго порядка (4), а значит и начальных условий должно быть два. Но почему? Первое, что нам необходимо знать - это начальная угловая скорость \(\omega\), так как

$$
\omega = \int_{0}^{t}{\varepsilon dt} + C \qquad (6)
$$
Нужно найти константу \(C\). При \(t = 0\) интеграл в (6) равен нулю, а угловая скорость равна начальной. То есть \(C = \omega_0\)
Теперь зная \(\omega\) можно найти и угол, но тогда появляется второе необходимое начальное условие - это сам угол \(\theta\):

$$
\theta = \int_{0}^{t}{\omega dt} + C \qquad (7)
$$
При \(t = 0\) угол \(\theta = C = \theta_0\)
Всё бы ничего, но решить (6) и уж тем более (7) не всегда получится. И наш случай как раз такой. Это значит, что мы не сможем таким способом получить точное решение. Но можно решать (6) и (7) численно. Как это делается и зачем для этого нужен компьютер?

Численное интегрирование


Теперь понятно, что для определённости системы, не зависимо от способа решения (6) и (7), обязательно нужно два условия - \(\theta\) и \(\omega\). Допустим, что угол \(\theta = 2\pi/3\). Начальную угловую скорость разумно задать нулевой (в этом случае мы просто отпускаем грузик, не толкая его). Теперь мы можем вычислить угловое ускорение \(\varepsilon\) с помощью (4). Угловое ускорение будет постоянно меняться, но если взять маленький промежуток времени \(\Delta t\) (\(\Delta t\) много меньше времени эксперимента), то в этом промежутке можно считать \(\varepsilon\) постоянной величиной, а значит изменение угловой скорости можно будет легко найти, ведь \(\varepsilon\) выносится из интеграла (6) как константа. После вычисления новой угловой скорости можно вычислить новый угол, опять же считая, что \(\omega\) с хорошей точностью не меняется в промежутке времени \(\Delta t\). После пересчёта всех необходимых параметров алгоритм повторяется.

Так, step-by-step, компьютер производит расчёт. Программа останавливается, например, когда \(s\) становится меньше миллиметра. Я нарисовал алгоритм работы программы для лучшего понимания:





Ну и наконец - результат работы алгоритма:




А вот что будет, если уменьшать радиус стержня. Видно, что количество "завитушек" увеличивается. При этом расстояние между ними будет одинаково и равно \(2\pi R\) (Когда нить уже не проскальзывает)






Друзья! Я очень благодарен вам за то, что вы интересуетесь моими работами, ведь каждый пост на сайте даётся очень непросто. Я буду рад любому отклику и поддержке с вашей стороны.







Если у вас остались вопросы или пожелания, то вы можете оставить комментарий (регистрироваться не нужно)

Гость из Сингапура:

Почему у вас \(dH = -A dt\)? Не должно ли быть \(dH = \frac{A dt^2}{2}\)?
Единицы \(А\) - \(m/s^2\) , а единицы \(t\) - \(s\), но единицы \(H\) - \(m\). Единицы не согласованы.
-------------------------------------------
В статье опечатка, уже исправил, большое спасибо!
Вместо \(A\) должно быть \(V\). Где \(V\) - это скорость большого груза.
Есть теорема (о связи производной и дифференциала), которая утверждает, что \(dH(t_0)(h) = -V(t_0)h + o(h)\) где \(o(h) \to 0\) при \(h \to 0\)
Минус перед производной стоит потому, что при уменьшении координаты \(y\) груза увеличивается \(H\).
При моделировании все было написано правильно.


Дата: 23-10-2018 в 20:47

Гость из Беларуси:

Здравствуйте. Подскажите, как писалось моделирование, то есть как выражались координаты x и y?
-------------------------------------------
Если вы знаете угол \(\theta\), расстояние \(s\) и \(H\), то координаты легко определяются из правил тригонометрии.




Зная параметры системы на предыдущем шаге, программа рассчитывает, что будет на следующем.


Дата: 26-10-2018 в 05:03

Гость из Сингапура:

Почему прекращается массовая масса? Извините, я не понял объяснений.
-------------------------------------------
Write it in english, please :)


Дата: 26-10-2018 в 21:27

Гость из Кореи:

Hi, I am a high school student in Korea. I am solving this problem using energy conservation.
I divided this problem to two situations.
1) Yellow box is stopped after falling.
2) The second is yellow box is falling.
I found O.D.E about First situations, but I don't know how to found about Second's O.D.E.....
Do you have solution about it?
-------------------------------------------
I don't think that energy conservation is applicable for this case at least because of the friction in second case.
By the way, the second case is more important than the first one :)


Дата: 01-11-2018 в 15:48

Гость из Кореи:

While the yellow box drops, do we have to use momentum conservation? Or use your method in the top?
-------------------------------------------
I took into account the conservation of angular momentum. It is important!


Дата: 07-11-2018 в 10:11

Гость из Кореи:

Can you express \(d\theta\) over time while yellow box drop for me...?
-------------------------------------------
$$
\frac{d\omega}{dt} = \frac{gsin\theta}{s}
$$
$$
d\theta = \omega dt
$$


Дата: 07-11-2018 в 17:06

Гость из Кореи:

Do you have equation about S, except (S+Rθ+H=L)?
-------------------------------------------
Nope


Дата: 08-11-2018 в 17:14

Гость из Бразилии:

Какую программу вы использовали для этой графики?
-------------------------------------------
AS3 by Adobe flash pro
You can also use javascript.


Дата: 11-11-2018 в 18:29

Тайвань:

Hello, I'm trying to write the simulation in Visual Basic according to your step. But I'm stuck in some problem.
1. You said that is important of the conservation angular momentum, but where should I write it in to my program.
2. My H value doesn't increase. I'm not sure whether I should write
\(v_2 = v_2 + a_2 dt\)
\(h = -v_2 dt\)
or
\(h = -(1/2 * a_2 dt^2)\)

Thank you!!
-------------------------------------------
1) You should take into account the conservation of angular momentum when you change the length \(l\): The angular velocity should be increased when length becomes smaller. \(\textbf{L}=[\textbf{r}\times\textbf{p}] = const\)
2) You should write \(h := h - v_2dt\)


Дата: 05-12-2018 в 13:03

Москва:

Здравствуйте. Вы записывает уравнение движения относительно точки касания верёвки с блоком, т. е. неинерциальной системы. При этом расстояние до этой точки уменьшается из-за намотки и проскальщывания нити - возникает сила кориолиса. Вы уверены, что здесь её не надо учитывать?
-------------------------------------------
Не уверен.
Учтите и проверьте :)


Дата: 20-12-2018 в 15:44

Москва:

Для каких начальных параметров Вы делали моделирование?
-------------------------------------------
А кто тоже любит старый добрый flash player?
По обозначениям все так же, как в статье.




Дата: 20-12-2018 в 18:25

Москва:

Как может сохраняться момент импульса, если он непрерывно изменяется из-за влияния силы тяжести?
-------------------------------------------
Сохранение момента я учитывал только во время небольшого изменения длины нити. В этот малый промежуток влияние силы тяжести стремится к нулю.


Дата: 20-12-2018 в 18:29

Китай:

Why do you think that angular momentum is conservation?
-------------------------------------------
I took into account the conservation of angular momentum during the change of the length. The effect of gravity tends to zero.


Дата: 21-12-2018 в 05:01

Беларусь:

dθ - что это за угол? между чем и чем он?


Дата: 21-12-2018 в 19:50

Guest from Canada:

I think you can go even more indepth if you consider the motion of the rope along the rod itself because we know that the rope cannot continue wrapping over itself, so it must wrap in a spiral pattern. I think this is actually an important part of this IYPT problem.


Дата: 27-12-2018 в 07:18

Guest from Canada:

I believe that it is important to consider the rope as wrapping in a spiral shape rather than over itself. This causes the rope to move along the rod rather than just tangentially. Any thoughts?


Дата: 27-12-2018 в 07:20

Анонимно:

Hello.
We tried to solve your differential equation using python3 and numpy, but even with the same parameters, we are not getting the same result as yours. Could you send us your source code so that we can examine the difference in the algorithm?
Our code is on https://www.codepile.net/pile/OVZ9mkEa


Дата: 14-01-2019 в 05:14

Анонимно:

help me please. I am student i will pass away soon if you help me i will not died.


Дата: 17-01-2019 в 16:15

Анонимно:

How was the conservation of angular momentum applied at the end?
-------------------------------------------
Nothing special: just like in the beginning


Дата: 21-01-2019 в 13:49

Анонимно:

Здравствуйте. Используя ваш алгоритм и ваши же начальные данные, результат получается совсем другим. Могли бы вы опубликовать код программы?


Дата: 09-03-2019 в 20:07

China:

hello,i am a post-graduate student in China. Can I recite your simulation flash in my report ? Certainly I will mark this source if i can get your admission.
-----------------------------------------------
Yes, of course!


Дата: 19-03-2019 в 17:51

Анонимно:

I seem to be really struggling with the program here. I am using VPython and I think I am implementing the algorithm correctly, except for the conservation of angular momentum part. I don't quite understand why angular momentum is conserved because isn't gravity applying a net torque on the system? (Also, doesn't I*omega clearly go from zero, to a non-zero value?). Or is angular momentum only conserved over the time interval dt, and we recalculate angular momentum each time (this is what I am attempting here, I don't know if I am doing it correctly however). I have included the functional loop of my code below, and the trajectory of the small mass is always a straight line back to the origin. Would you be able to take a look at it and let me know what I am missing?

while(s > 0.002):
rate(100)
T = m1*(omega**2)*s - m1*g*cos(theta)
alpha = (g*sin(theta))/s
if m2*g >= T*(e**(mu*theta)):
A_y = (T*(e**(mu*theta))/m2) - g
else:
A_y = 0
omega += alpha*dt
theta += omega*dt
v_y += v_y + A_y*dt
Y += abs(v_y)*dt
ds = -r*(omega*dt) - abs(v_y)*dt
s = s + ds
L = (m1*s**2)*omega
omega = L/(m1*s**2)

m1 is the lighter mass, m2 is the heavier mass. Also, I am using g = +9.81 (and including a minus sign wherever necessary).

Thank you for your help in advance!


Дата: 01-04-2019 в 19:48


Для будущих авторов | 12.10.18: Если вы хотите стать автором статей на сайте и получить подтвержденный аккаунт, то обращайтесь на почту! support@ilinblog.ru

Обновления | 21.08.18: Добавлена возможность комментировать статьи. Сайт адаптирован под мобильные устройства.

Обновления | 19.01.18: Добавлена возможность добавления математических формул в статьи посредством языка latex. Пример использования тут. Также добавлена возможность редактирования статей.

Информация о пользователях | 28.10.17: Расширена функциональность страницы пользователей, теперь можно добавить статус и личную информацию.

Обновления | 30.08.17: С недавнего времени появилась возможность создавать свои собственные публикации, которые будут отображаться на странице пользователя. Авторы будут иметь возможность создавать посты на главную страницу сайта.



500
800
65
88