Пример

Пример

В качестве упражнения смоделируем полёт пушечного ядра.

Пушка совершает выстрел сферическим ядром радиуса \(r\) и массы \(m\) под углом \(\theta_0\) к горизонту с начальной скоростью \(v_0\). Смоделировать полёт снаряда. Считать, что на ядро действует лобовое сопротивление, т.е. сила сопротивления направлена против скорости движения \(\vec{v}\), её величина пропорциональна характерной площади \(S\), плотности среды \(\rho\) и квадрату скорости \(v^2\):

\[ \vec{F} = - C_{F}\frac{\rho v}{2}S \vec{v}. \]

\(C_{F}\)безразмерный аэродинамический коэффициент сопротивления, который для сферы равен \(0.47\). \(S\) считают равным площади поперечного сечения в случае сферы (\(\pi r^2\)). Обычно стандартной величиной плотности воздуха на уровне моря в соответствии с Международной стандартной атмосферой принимается значение 1.2250 кг/м³ (источник - Wiki).

Совместим начало координат с пушкой, ось \(Ox\) направим горизонтально в сторону полёта снаряда, а ось \(Oy\) вертикально.

../_images/cannonball.png

Тогда уравнения движения примут вид

\[ m \dfrac{d \vec{v}}{dt} = m\vec{g} + \vec{F}_{conp} \]

или

\[\begin{split} \begin{cases} m \dfrac{dv_x}{dt} = - C_{F}\dfrac{\rho v}{2}S v_x, \\ m \dfrac{dv_y}{dt} = - mg - C_{F}\dfrac{\rho v}{2}S v_y. \end{cases} \end{split}\]

Найдите траекторию полёта снаряда.

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)