Пример
Пример¶
В качестве упражнения смоделируем полёт пушечного ядра.
Пушка совершает выстрел сферическим ядром радиуса \(r\) и массы \(m\) под углом \(\theta_0\) к горизонту с начальной скоростью \(v_0\). Смоделировать полёт снаряда. Считать, что на ядро действует лобовое сопротивление, т.е. сила сопротивления направлена против скорости движения \(\vec{v}\), её величина пропорциональна характерной площади \(S\), плотности среды \(\rho\) и квадрату скорости \(v^2\):
\(C_{F}\) — безразмерный аэродинамический коэффициент сопротивления, который для сферы равен \(0.47\). \(S\) считают равным площади поперечного сечения в случае сферы (\(\pi r^2\)). Обычно стандартной величиной плотности воздуха на уровне моря в соответствии с Международной стандартной атмосферой принимается значение 1.2250 кг/м³ (источник - Wiki).
Совместим начало координат с пушкой, ось \(Ox\) направим горизонтально в сторону полёта снаряда, а ось \(Oy\) вертикально.
Тогда уравнения движения примут вид
или
Найдите траекторию полёта снаряда.
from scipy.constants import g, degree, pi
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
# params
theta0 = 45 * degree
v0 = 10 # mps
m = 1 # kg
r = 0.1 # m
v0x = v0 * np.cos(theta0)
v0y = v0 * np.sin(theta0)
# constants
C_F = 0.47
S = pi * r**2
rho = 1.2250 # kg/m^3
C_resist = C_F * S * rho / 2. # F_resist = - C_resist * ||v|| * v
def resistance_free_solution(t, v0x=v0x, v0y=v0y):
...
def plot_solution(x_nores, y_nores, x=None, y=None):
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(x_nores, y_nores, label="resistance free")
ax.set_xlabel("distance (m)")
ax.set_ylabel("height (m)")
if x is not None:
ax.plot(x, y, label="with resistance")
ax.legend()
plt.show(fig)