К основному контенту

Движение тела в поле тяжести с учётом сопротивления воздуха

Это творческое задание для мастер-класса по информатике для школьников при ДВФУ.
Цель задания - выяснить, как изменится траектория тела, если учитывать сопротивление воздуха. Также необходимо ответить на вопрос, будет ли дальность полёта по-прежнему достигать максимального значения при угле бросания в 45°, если учитывать сопротивление воздуха.
В разделе "Аналитическое исследование" изложена теория. Этот раздел можно пропустить, но он должен быть, в основном, понятным для вас, потому что большую часть из этого вы проходили в школе.
В разделе "Численное исследование" содержится описание алгоритма, который необходимо реализовать на компьютере. Алгоритм простой и краткий, поэтому все должны справиться.

Аналитическое исследование

Введём прямоугольную систему координат так, как показано на рисунке. В начальный момент времени тело массой m находится в начале координат. Вектор ускорения свободного падения направлен вертикально вниз и имеет координаты (0, -g).
 - вектор начальной скорости. Разложим этот вектор по базису: . Здесь , где - модуль вектора скорости, - угол бросания.

Запишем второй закон Ньютона: .
Ускорение в каждый момент времени есть (мгновенная) скорость изменения скорости, то есть производная от скорости по времени: .
Следовательно, 2-й закон Ньютона можно переписать в следующем виде:
, где - это равнодействующая всех сил, действующая на тело.
Так как на тело действуют сила тяжести и сила сопротивления воздуха, то
.

Мы будем рассматривать три случая:
1) Сила сопротивления воздуха равна 0: .
2) Сила сопротивления воздуха противоположно направлена с вектором скорости, и её величина пропорциональна скорости: .
3) Сила сопротивления воздуха противоположно направлена с вектором скорости, и её величина пропорциональна квадрату скорости: .

Вначале рассмотрим 1-й случай.
В этом случае , или .
Запишем это равенство в скалярном виде:

Из следует, что (равноускоренное движение).
Так как (r - радиус-вектор), то .
Отсюда .
Эта формула есть не что иное, как знакомая вам формула закона движения тела при равноускоренном движении.
Так как , то .
Учитывая, что и , получаем из последнего векторного равенства скалярные равенства:
Проанализируем полученные формулы.
Найдём время полёта тела. Приравняв y к нулю, получим
Дальность полёта равна значению координаты x в момент времени t0:
Из этой формулы следует, что максимальная дальность полёта достигается при .
Теперь найдём уравнение трактории тела. Для этого выразим t через x
   
и подставим полученное выражение для t в равенство для y.
Полученная функция y(x) -- квадратичная функция, её графиком является парабола, ветви которой направлены вниз.
Про движение тела, брошенного под углом к горизонту (без учёта сопротивления воздуха), рассказывается в этом видеоролике.

Теперь рассмотрим второй случай: .
Второй закон приобретает вид ,
отсюда .
Запишем это равенство в скалярном виде:

Мы получили два линейных дифференциальных уравнения.
Первое уравнение имеет решение
в чём можно убедиться, подставив данную функцию в уравнение для vx и в начальное условие .
Здесь e = 2,718281828459... -- число Эйлера.
Второе уравнение имеет решение
Так как , , то при наличии сопротивления воздуха движение тела стремится к равномерному, в отличие от случая 1, когда скорость неограниченно увеличивается.
В следующем видеоролике говорится, что парашютист сначала движется ускоренно, а потом начинает двигаться равномерно (даже до раскрытия парашюта).

Найдём выражения для x и y.
Так как x(0) = 0, y(0) = 0, то

Отсюда

Нам осталось рассмотреть случай 3, когда .
Второй закон Ньютона имеет вид
, или .
В скалярном виде это уравнение имеет вид:
Это система нелинейных дифференциальных уравнений. Данную систему не удаётся решить в явном виде, поэтому необходимо применять численное моделирование.


Численное исследование

В предыдущем разделе мы увидели, что в первых двух случаях закон движения тела можно получить в явном виде. Однако в третьем случае необходимо решать задачу численно. При помощи численных методов мы получим лишь приближённое решение, но нас вполне устроит и небольшая точность. (Число π или квадратный корень из 2, кстати, нельзя записать абсолютно точно, поэтому при расчётах берут какое-то конечное число цифр, и этого вполне хватает.)

Будем рассматривать второй случай, когда сила сопротивления воздуха определяется формулой. Отметим, что при k = 0 получаем первый случай.
Скорость тела подчиняется следующим уравнениям:

В левых частях этих уравнений записаны компоненты ускорения .
Напомним, что ускорение есть (мгновенная) скорость изменения скорости, то есть производная от скорости по времени.
В правых частях уравнений записаны компоненты скорости. Таким образом, данные уравнения показывают, как скорость изменения скорости связана со скоростью.

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

Заменим в уравнениях ускорение (мгновенную скорость изменения скорости) на среднюю скорость изменения скорости, рассматривая движение тела на промежутке времени :
Теперь подставим полученные аппроксимации в наши уравнения.
Отсюда
Полученные формулы позволяют нам вычислить значения функций в следующем узле сетки, если известны значения этих функций в предыдущем узле сетки.
При помощи описанного метода мы можем получить таблицу приближённых значений компонент скорости.

Как найти закон движения тела, т.е. таблицу приближённых значений координат x(t), y(t)? Аналогично!
Имеем
Заменив мгновенную скорость на среднюю скорость на промежутке времени , получим
Отметим, что в правых частях уравнений можно взять полусумму значений vx (vy) в точках t и t + Δt. Возможно, так точность будет выше.
Отсюда
По этим формулам мы можем вычислить таблицу приближённых значений функций x(t) и y(t) в узлах сетки.


Как теперь реализовать полученный алгоритм? Да очень просто!
Заведём массивы
vx, vy: array[0..MAX_N] of extended;
x, y: array[0..MAX_N] of extended;
Значение vx[j] равняется значению функции , для других массивов аналогично.
Теперь остаётся написать цикл, внутри которого мы будем вычислять vx[j+1] через уже вычисленное значение vx[j], и с остальными массивами то же самое. Цикл будет по j от 1 до N.
Не забудьте инициализировать начальные значения vx[0], vy[0], x[0], y[0] по формулам , x0 = 0, y0 = 0.


В Паскале и Си для вычисления синуса и косинуса имеются функции sin(x), cos(x). Обратите внимание, что эти функции принимают аргумент в радианах.

Вам необходимо построить график движения тела при k = 0 и k > 0 и сравнить полученные графики. Графики можно построить в Excel.
Отметим, что расчётные формулы настолько просты, что для вычислений можно использовать один только Excel и даже не использовать язык программирования.
Однако в дальнейшем вам нужно будет решить задачу в CATS, в которой нужно вычислить время и дальность полёта тела, где без языка программирования не обойтись.

Обратите внимание, что вы можете протестировать вашу программу и проверить ваши графики, сравнив результаты вычислений при k = 0 с точными формулами, приведёнными в разделе "Аналитическое исследование".

Поэкспериментируйте со своей программой. Убедитесь в том, что при отсутствии сопротивления воздуха (k = 0) максимальная дальность полёта при фиксированной начальной скорости достигается при угле в 45°.
А с учётом сопротивления воздуха? При каком угле достигается максимальная дальность полёта?

На рисунке представлены траектории тела при v0 = 10 м/с, α = 45°, g = 9,8 м/с2, m = 1 кг, k = 0 и 1, полученные при помощи численного моделирования при Δt = 0,01.



Здесь вы можете ознакомиться с замечательной работой 10-классников из г. Троицка, представленной на конференции "Старт в науку" в 2011 г. Работа посвящена моделированию движения теннисного шарика, брошенного под углом к горизонту (с учетом сопротивления воздуха). Применяется как численное моделирование, так и натурный эксперимент.

Таким образом, данное творческое задание позволяет познакомиться с методами математического и численного моделирования, которые активно используются на практике, но мало изучаются в школе. К примеру, данные методы использовались при реализации атомного и космического проектов в СССР в середине XX века.

Комментарии

Unknown написал(а)…
Помогите, не понимаю, откуда в решении, там, где линейная зависимость сопротивления от скорости, в выражении для Vy появился множитель (1-е^(-k/m)*t)?
Глеб Гренкин написал(а)…
Уравнение для Vy, так же как и уравнение для Vx - это линейные дифференциальные уравнения.
Чтобы было немного понятнее, запишем эти уравнения в более привычной форме:
dVx/dt + (k/m)Vx = 0,
dVy/dt + (k/m)Vy = -g.
В уравнении для Vx правая часть равна 0, такое уравнение называется однородным. Уравнение для Vy неоднородное.
Для линейных уравнений справедлив принцип линейной суперпозиции решений. А именно, допустим, что функции u(x) и v(x) подчиняются уравнениям
u'(x) + Au(x) = f(x) и v'(x) + Av(x) = g(x)
и начальным условиям u(0) = u0, v(0) = v0.
Тогда сумма z(x) = u(x) + v(x) будет удовлетворять уравнению
z'(x) + Az(x) = f(x) + g(x)
и начальному условию z(0) = u0 + v0.

Применительно к данному примеру:
Vy = U + W,
dU/dt + (k/m)U = 0, U(0) = Vy0,
dW/dt + (k/m)W = g, W(0) = 0.
Решение уравнения для U такое же, как и для уравнения для Vx, поэтому первое слагаемое в формуле для Vy такое же, как и в формуле для Vx.
Чтобы найти второе слагаемое, нужно решить уравнение для W.
Тут вы можете посмотреть стандартные методы интегрирования линейных диффуров - http://vk.com/wall-50476096_5134.
Можно ли как-то сразу написать решение для W с учётом нулевого начального условия, я не знаю.
PS Когда вы решаете линейный диффур, его не обязательно раскладывать в сумму U и W (просто применяете общую формулу - и всё). Я попытался пояснить, откуда берутся такие слагаемые.
Unknown написал(а)…
Спасибо :D
Анонимный написал(а)…
Подскажите, а что будет во время удара о землю, и каковы будут потери энергии при ударе, если его не считать упругим?
Глеб Гренкин написал(а)…
Понятия не имею.
Глеб Гренкин написал(а)…
По-моему, есть игра, когда кидают шарик, и его нужно поймать тарелкой, к которой он прилипает. Если тарелку положить на землю и шарик на неё упадёт, то он прилипнет. А если кинуть обычный теннисный шарик, то он несколько раз отскочет от земли. То есть это зависит от свойств поверхности земли и тела, например, трения.
Анонимный написал(а)…
Этот комментарий был удален автором.
Глеб Гренкин написал(а)…
Не понял, про какую конкретно формулу вы спрашиваете.
Ускорение, создаваемое силой сопротивления, направлено противоположно скорости. Отсюда
dv_x/dt = -(k/m)*v_x и dv_y/dt = -(k/m)*v_y - g.
MaxJohnson написал(а)…
Спасибо, это я понял. Я не совсем корректно задал вопрос. Дело в том, что я пытаюсь рассчитать траекторию полета мяча для голфа с учетом сопротивления воздуха. На данный момент я рассчитываю траекторию с учетом линейного сопротивления, но у меня она получается очень странная и не реалистичная. Масса мяча 0.0459 кг, начальная скорость 38.1 м/c, коэффициент k примерно равен 0.003 при придельной скорости равной 15.2 м/с. Если вы могли бы мне помочь я бы был вам очень благодарен. Заранее спасибо!
Глеб Гренкин написал(а)…
Можете прислать программу. Вдруг что-нибудь найду.
http://grenkin.github.io/email.png
Unknown написал(а)…
Этот комментарий был удален автором.
Глеб Гренкин написал(а)…
Презентация по теме: https://vk.com/wall-183726611_262
Анонимный написал(а)…
Как рассчитать время полёта тела с учётом сопротивления воздуха из этой формулы?