В этот раз на примере петлевого маятника - интересной и красивой механической задачи из предстоящего в 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.
_

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

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



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



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

N=Tsin(dθ/2)+(T+dT)sin(dθ/2)N = Tsin(d\theta/2) + (T+dT)sin(d\theta/2)

При малых углах sin(α)≈αsin(\alpha) \approx \alpha. То есть

N=Tdθ+dTdθ/2N = Td\theta + dTd\theta/2

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

T+μTdθ=T+dTT + \mu Td\theta = T+dT

Отсюда приходим к выражению

dTT=μdθ\frac{dT}{T} = \mu d\theta

при интегрировании которого получаем формулу Эйлера:

T1=Teμθ(1)T_1 = T e^{\mu\theta} \qquad (1)

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



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

{MAx=0MAy=Teμθ−Mg(2)\left\{ \begin{array}{ccc} MA_x = 0\\ MA_y = Te^{\mu\theta} - Mg\\ \end{array} \right. \qquad (2)

Найти силу натяжения TT достаточно просто: нужно только заметить, что маленький груз вращается по окружности относительно точки касания нити со стержнем. Это значит, что в проекции на саму нить сумма сил должна обеспечивать центростремительное ускорение aц=vt2/sa_ц = v_{t}^2/s:

T+mgcosθ=mvt2sT + mg cos\theta = m\frac{v_{t}^2}{s}

Отсюда

T=mvt2s−mgcosθ(3)T = m\frac{v_{t}^2}{s} - mg cos\theta \qquad (3)

Где vt=dθdtsv_{t} = \frac{d\theta}{dt}s
Закон изменения угла θ\theta можно найти, записав 2 закон Ньютона для вращательного движения:

Id(dθdt)dt=∑[r×F]I\frac{d\left(\frac{d\theta}{dt}\right)}{dt} = \sum{[r\times F]}

Где II - момент инерции и в нашем случае I=ms2I = ms^2, а сумма из правой части - это общий момент сил относительно выделенной точки.
В итоге, получаем следующую формулу (вращение относительно точки контакта нити со стержнем):

d(dθdt)dt=gsinθs(4)\frac{d\left(\frac{d\theta}{dt}\right)}{dt} = \frac{gsin\theta}{s} \qquad (4)

Как изменяются координаты большого груза и угол θ\theta мы поняли. Осталось только найти закон изменения ss, который получается из условия нерастяжимости нити. Допустим, что длинна нити равна LL, тогда L=s+Rθ+HL = s + R\theta + H, где HH - это расстояние от большого груза до точки касания. В дифференциальном виде:

ds=−dH−Rdθ(5)ds = -dH - Rd\theta \qquad (5)

Мы получили, что изменение угла и HH приводит к изменению расстояния от маленького грузика до точки контакта нитки со стержнем. В свою очередь это приводит к увеличению угловой скорости из-за закона сохранения момента импульса, который в нашем случае равен IωI\omega. То есть при уменьшении ss угловая скорость малого груза будет сильно расти (ω∼1/s2\omega \sim 1/s^2). Это и приводит к закручиванию груза.
В принципе это всё, что нужно для численного моделирования.

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

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

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

ω=∫0tεdt+C(6)\omega = \int_{0}^{t}{\varepsilon dt} + C \qquad (6)

Нужно найти константу CC. При t=0t = 0 интеграл в (6) равен нулю, а угловая скорость равна начальной. То есть C=ω0C = \omega_0
Теперь зная ω\omega можно найти и угол, но тогда появляется второе необходимое начальное условие - это сам угол θ\theta:

θ=∫0tωdt+C(7)\theta = \int_{0}^{t}{\omega dt} + C \qquad (7)

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

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

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

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



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



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