Интегрирование уравнения материальной движения точки и верификация решения
Численное решение задачи баллистики с учётом сопротивления воздуха и верификация полученного решения в среде Python, используя теорему об изменении кинетической энергии.
import numpy as np
from scipy.integrate import solve_ivp
from collections import namedtuple
import matplotlib.pyplot as plt
import scipy.integrate as integrate
Параметры модели
Параметры модели задаются при помощи именованного кортежа
p = namedtuple("params", "m Cx rho Sm g")
# Масса
p.m = 5;
# Ускорение свободного падения
p.g = 9.807
# Аэродинамический коэффициент лобового сопротивления
p.Cx = 1.0;
# Характерная площадь
p.Sm = np.pi*0.2**2/4
# Плотность воздуха
p.rho = 1.2
Функция правых частей
Функция правых частей описывает дифференциальное уравнение второго порядка (как систему двух уравнений первого порядка) движения материальной точки под действием силы тяжести в однородном поле и аэродинамической силы сопротивления, направленной в противоположную сторону скорости: \[m \frac{d \vec{v}}{dt} = - m \vec{g} - \frac{\vec{v}}{|\vec v|} C_x S_m q.\]
def dydt(t, q, p):
r = q[0:2];
v = q[2:4];
# Скоростной напор
dynamic_pressure = p.rho*np.dot(v,v)*0.5;
# Аэродинамическая сила
Fa = -v/np.sqrt(np.dot(v,v))*p.Cx*p.Sm*dynamic_pressure;
# Сила тяжести
G = np.array([0,-p.m*p.g]);
# Координатный столбец главного вектора сил, действующего на тело
F = G + Fa;
# Ускорение
a = F/p.m;
return np.hstack([v,a])
Функция для остановки процесса интегрирования при h = 0
Функция-“детектор”, передаваемая в интегратор (параметр events), для определения времени достижения нулевой высоты и остановки процесса интегрирования.
def event_h_eq_0(t, q):
# Возвращаем высоту (контролируемая функция)
return q[1]
# функция-детектор срабатывает при условии h = 0 при движении "вниз" (при убывании контролируемой функции)
event_h_eq_0.direction = -1
# функция-детектор останавливает процесс
event_h_eq_0.terminal = True
Начальные условия и запуск численного интегрирования
v0 = 30
phi0 = np.deg2rad(45)
r0_vec = [0,0]
v0_vec = [v0*np.cos(phi0),v0*np.sin(phi0)]
q0 = np.hstack([r0_vec,v0_vec])
sol = solve_ivp(lambda t,q: dydt(t,q,p), [0, 5], q0, method='RK45', rtol = 1e-8, max_step = 0.5, events = [event_h_eq_0])
Вид решения
message: A termination event occurred.
success: True
status: 1
t: [ 0.000e+00 1.445e-02 1.590e-01 5.280e-01 9.195e-01
1.343e+00 1.790e+00 2.244e+00 2.703e+00 3.187e+00
3.687e+00 4.053e+00]
y: [[ 0.000e+00 3.064e-01 ... 6.725e+01 7.295e+01]
[ 0.000e+00 3.053e-01 ... 6.322e+00 -8.882e-15]
[ 2.121e+01 2.118e+01 ... 1.585e+01 1.535e+01]
[ 2.121e+01 2.104e+01 ... -1.577e+01 -1.880e+01]]
sol: None
t_events: [array([ 4.053e+00])]
y_events: [array([[ 7.295e+01, -8.882e-15, 1.535e+01, -1.880e+01]])]
nfev: 68
njev: 0
nlu: 0
Решение
Траектория
plt.plot(sol.y[0],sol.y[1],'.-'); plt.grid(ls=':'); plt.xlabel('x, м'); plt.ylabel('y, м');
Верификация
Теорема об изменении кинетической энергии \[T_k - T_0 = A_a + A_g,\]
\(A_a\) - работа аэродинамической силы сопротивления, \(A_g = 0\) т.к. изменение высоты равно нулю. Работа силы сопротивления: \[dA_a = \vec{F} \cdot \vec{V} dt, \quad A_a = \int_0^{t_k} \vec{F} \cdot \vec{V} dt\]
v = sol.y[2:4].T
v
array([[ 21.21320344, 21.21320344],
[ 21.17864132, 21.03701046],
[ 20.84522028, 19.29945152],
[ 20.0850104 , 15.04180359],
[ 19.39798133, 10.75264161],
[ 18.76020505, 6.30836944],
[ 18.17157011, 1.79943502],
[ 17.62244665, -2.6357364 ],
[ 17.07976753, -6.99151935],
[ 16.49182395, -11.41599742],
[ 15.84673062, -15.77466726],
[ 15.34602733, -18.80247627]])
Определим модуль скорости:
v_norm = np.sqrt(np.sum(v*v,axis=1))
v_norm
array([30. , 29.85114164, 28.40760528, 25.09309662, 22.17884085,
19.79244347, 18.26044705, 17.81846605, 18.45534615, 20.05754856,
22.35976294, 24.27001584])
Скоростной напор:
dynamic_pressure = p.rho*(v_norm**2)*0.5
Единичный вектор направления скорости точки:
ev = v/np.reshape(v_norm,[-1,1])
ev
array([[ 0.70710678, 0.70710678],
[ 0.70947509, 0.70473052],
[ 0.73379013, 0.67937622],
[ 0.80041976, 0.59943991],
[ 0.87461655, 0.48481531],
[ 0.94784684, 0.31872616],
[ 0.99513282, 0.09854277],
[ 0.98899909, -0.14792162],
[ 0.92546449, -0.37883437],
[ 0.8222253 , -0.56916215],
[ 0.7087164 , -0.70549349],
[ 0.63230397, -0.7747204 ]])
Аэродинамическая сила
Fa = -ev*p.Cx*p.Sm*np.reshape(dynamic_pressure,[-1,1])
Мощность аэродинамической силы
dW = np.sum(Fa*v,axis=1)
plt.figure(figsize=[8,4])
plt.plot(sol.t,dW,'.-')
plt.ylim([-550,0])
plt.fill_between(sol.t,dW, alpha=0.3)
plt.xlabel('t, c')
plt.ylabel('Мощность, Вт')
Работа аэродинамической силы
Work_aero = integrate.simps(dW,x=sol.t)
Work_aero
-777.4154150947423
Начальная кинетическая энергия
T0 = v_norm[ 0]**2*p.m*0.5
Конечная кинетическая энергия
TK = v_norm[-1]**2*p.m*0.5
Изменение кинетической энергии
delta_T = TK - T0
delta_T
-777.4158275842642
Погрешность
Относительная погрешность \[\varepsilon = \frac{\Delta T - A_a}{\Delta T} 100 \%\]
100*(delta_T-Work_aero)/delta_T
5.305905890102044e-05
Абсолютная погрешность \[a = \Delta T - A_a\]
delta_T-Work_aero
-0.0004124895218637903
