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

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

В разделе "Аналитическое исследование" изложена теория. Этот раздел можно пропустить, но он должен быть, в основном, понятным для вас, потому что большую часть из этого вы проходили в школе.
В разделе "Численное исследование" содержится описание алгоритма, который необходимо реализовать на компьютере. Алгоритм простой и краткий, поэтому все должны справиться.

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

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

Запишем второй закон Ньютона: .
Ускорение в каждый момент времени есть (мгновенная) скорость изменения скорости, то есть производная от скорости по времени: .

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

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

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

Из следует, что (равноускоренное движение).
Так как (r - радиус-вектор), то .
Отсюда .
Эта формула есть не что иное, как знакомая вам формула закона движения тела при равноускоренном движении.
Так как , то .
Учитывая, что и , получаем из последнего векторного равенства скалярные равенства:

Проанализируем полученные формулы.
Найдём время полёта тела. Приравняв y к нулю, получим

Дальность полёта равна значению координаты x в момент времени t₀:

Из этой формулы следует, что максимальная дальность полёта достигается при .
Теперь найдём уравнение трактории тела. Для этого выразим t через x

   

и подставим полученное выражение для t в равенство для y.

Полученная функция y(x) -- квадратичная функция, её графиком является парабола, ветви которой направлены вниз.
Про движение тела, брошенного под углом к горизонту (без учёта сопротивления воздуха), рассказывается в этом видеоролике.

Теперь рассмотрим второй случай: .

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

Мы получили два линейных дифференциальных уравнения.
Первое уравнение имеет решение

в чём можно убедиться, подставив данную функцию в уравнение для v_x и в начальное условие .
Здесь e = 2,718281828459... -- число Эйлера.
Второе уравнение имеет решение

Так как , , то при наличии сопротивления воздуха движение тела стремится к равномерному, в отличие от случая 1, когда скорость неограниченно увеличивается.
В следующем видеоролике говорится, что парашютист сначала движется ускоренно, а потом начинает двигаться равномерно (даже до раскрытия парашюта).

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

Отсюда

Нам осталось рассмотреть случай 3, когда .
Второй закон Ньютона имеет вид
, или .
В скалярном виде это уравнение имеет вид:

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

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

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

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

Скорость тела подчиняется следующим уравнениям:

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

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

Наша задача -- приближённо вычислить значения в узлах сетки.

Заменим в уравнениях ускорение (мгновенную скорость изменения скорости) на среднюю скорость изменения скорости, рассматривая движение тела на промежутке времени :

Теперь подставим полученные аппроксимации в наши уравнения.

Отсюда

Полученные формулы позволяют нам вычислить значения функций в следующем узле сетки, если известны значения этих функций в предыдущем узле сетки.

При помощи описанного метода мы можем получить таблицу приближённых значений компонент скорости.

Как найти закон движения тела, т.е. таблицу приближённых значений координат x(t), y(t)? Аналогично!
Имеем

Заменив мгновенную скорость на среднюю скорость на промежутке времени , получим

Отметим, что в правых частях уравнений можно взять полусумму значений v_x (v_y) в точках 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] по формулам , x₀ = 0, y₀ = 0.


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

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

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

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

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

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

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