Петлевой маятник - численное моделирование. ТЮФ 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
$$
при интегрировании которого получаем формулу Эйлера:

Отлично, ведь одна из подзадач решена! С помощью результата (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\) (Когда нить уже не проскальзывает)

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