Как построить движение снаряда под действием силы тяжести, плавучести и сопротивления воздуха?

Question

Я пытаюсь составить график движения снаряда массой под действием силы тяжести, плавучести и сопротивления. В основном, я хочу показать влияние плавучести и силы сопротивления на расстояние полета, время полета и изменение скорости на графике.

import matplotlib.pyplot as plt
import numpy as np

V_initial = 30 # m/s
theta = np.pi/6 # 30
g = 3.711
m =1
C = 0.47
r = 0.5
S = np.pi*pow(r, 2)
ro_mars = 0.0175
t_flight = 2*(V_initial*np.sin(theta)/g)
t = np.linspace(0, t_flight, 200)

# Drag force
Ft = 0.5*C*S*ro_mars*pow(V_initial, 2)

# Buoyant Force
Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))

x_loc = []
y_loc = []

for time in t:
    x = V_initial*time*np.cos(theta)
        y = V_initial*time*np.sin(theta) - (1/2)*g*pow(time, 2)
    x_loc.append(x)
    y_loc.append(y)

x_vel = []
y_vel = []
for time in t:
    vx = V_initial*np.cos(theta)
    vy = V_initial*np.sin(theta) - g*time
    x_vel.append(vx)
    y_vel.append(vy)


v_ch = [pow(i**2+ii**2, 0.5) for i in x_vel for ii in y_vel]

tau = []
for velocity in v_ch:
    Ft = 0.5*C*S*ro_mars*pow(velocity, 2)
        tau.append(Ft)

buoy = []
for velocity in v_ch:
    Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))
    buoy.append(Fb)

после этого момента я не мог понять, как построить движение снаряда под действием этих сил. Другими словами, я пытаюсь сравнить движение снаряда при трех обстоятельствах

1. Масса только под действием силы тяжести
2. Масса под действием силы тяжести и сопротивления воздуха
3. Масса под действием силы тяжести, сопротивления воздуха и плавучести

Answer 1

Используя векторную нотацию, и odeint.

import numpy as np
from scipy.integrate import odeint
import scipy.constants as SPC
import matplotlib.pyplot as plt

V_initial = 30 # m/s
theta = np.pi/6 # 30
g = 3.711
m = 1 # I assume this is your mass
C = 0.47
r = 0.5
ro_mars = 0.0175

t_flight = 2*(V_initial*np.sin(theta)/g)
t = np.linspace(0, t_flight, 200)

pos0 = [0, 0]
v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]

def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):
    x, y, x_prime, y_prime = vector

    # volume and surface
    V = np.pi * 4/3 * r**3
    S = np.pi*pow(r, 2)

    # net weight bouyancy
    if apply_bouyancy:
        Fb = (ro_mars * V - m) * g *np.array([0,1])
    else:
        Fb = -m  * g * np.array([0,1])

    # velocity vector
    v = np.array([x_prime, y_prime])

    # drag force - corrected to be updated based on current velocity
#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)
    if apply_resistance:
        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)
    else:
        Ft = np.array([0, 0])

    # resulting acceleration
    x_prime2, y_prime2 = (Fb + Ft) / m

    return x_prime, y_prime, x_prime2, y_prime2

sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))
plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')
plt.legend(loc='best')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

Обратите внимание, что я исправил вашу силу сопротивления, чтобы использовать фактическую (не начальную) скорость, я не знаю, была ли это ваша ошибка или она была умышленной.

Также, пожалуйста, ознакомьтесь с документацией для odeint, чтобы лучше понять, как превратить ODE второго порядка (как в вашей задаче) в вектор ODE первого порядка.

Чтобы удалить сопротивление воздуха или сопротивление, установите apply_bouyancy и apply_resistance на True или False, добавив их к args=(...)

Answer 2

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

F_x = tau_x 
F_y = tau_y + bouyancy + gravity

Где tau_x и tau_y - силы сопротивления в направлениях x и y соответственно. Скорости v_x и v_y тогда определяются как

v_x = v_x + (F_x / (2 * m)) * dt
v_y = v_y + (F_y / (2 * m)) * dt

Таким образом, позиции x и y, r_x и r_y в любое время t задаются суммой

r_x = r_x + v_x * dt
r_y = r_y + v_y * dt

В обоих случаях это должно быть оценено от 0 до t для некоторого dt, где dt * n = t, если n - количество шагов суммирования.

r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2
r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2

Весь расчет может быть выполнен в две строки,

r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2
r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2

За исключением того, что v_x и v_y требуют обновления на каждом временном шаге. Чтобы повторить это и вычислить позиции x и y в диапазоне времен, вы можете просто следовать приведенному ниже (отредактированному) примеру.

Следующий код включает исправления для предотвращения отрицательных позиций y, поскольку данное значение g предназначено для поверхности или Марса. Я полагаю, что это уместно - когда вы нажмете ноль y и попытаетесь продолжить движение, вы можете в конечном итоге с быстрой незапланированной разборкой, как мы, физики, называем это.

Редактировать

В ответ на отредактированный вопрос следующий пример был изменен для отображения всех трех запрошенных случаев: гравитация, гравитация плюс сопротивление и гравитация плюс сопротивление и плавучесть. Код настройки графика также был добавлен

Полный пример (отредактировано)

import numpy as np
import matplotlib.pyplot as plt

def projectile(V_initial, theta, bouyancy=True, drag=True):
    g = 9.81
    m = 1
    C = 0.47
    r = 0.5
    S = np.pi*pow(r, 2)
    ro_mars = 0.0175

    time = np.linspace(0, 100, 10000)
    tof = 0.0
    dt = time[1] - time[0]
    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))
    gravity = -g * m
    V_ix = V_initial * np.cos(theta)
    V_iy = V_initial * np.sin(theta)
    v_x = V_ix
    v_y = V_iy
    r_x = 0.0
    r_y = 0.0
    r_xs = list()
    r_ys = list()
    r_xs.append(r_x)
    r_ys.append(r_y)
    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.
    # Just make sure you use sufficiently small dt (dt is change in time between steps)
    for t in time:
        F_x = 0.0
        F_y = 0.0
        if (bouyancy == True):
            F_y = F_y + bouy
        if (drag == True):
            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)
            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)
        F_y = F_y + gravity

        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2
        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2
        v_x = v_x + (F_x / m) * dt
        v_y = v_y + (F_y / m) * dt
        if (r_y >= 0.0):
            r_xs.append(r_x)
            r_ys.append(r_y)
        else:
            tof = t
            r_xs.append(r_x)
            r_ys.append(r_y)
            break

    return r_xs, r_ys, tof

v = 30
theta = np.pi/4

fig = plt.figure(figsize=(8,4), dpi=300)
r_xs, r_ys, tof = projectile(v, theta, True, True)
plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")
r_xs, r_ys, tof = projectile(v, theta, False, True)
plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")
r_xs, r_ys, tof = projectile(v, theta, False, False)
plt.plot(r_xs, r_ys, 'k:', label="Gravity")
plt.title("Trajectory", fontsize=14)
plt.xlabel("Displacement in x-direction (m)")
plt.ylabel("Displacement in y-direction (m)")
plt.ylim(bottom=0.0)
plt.legend()
plt.show()

Обратите внимание, что это сохраняет и возвращает время пролета в переменной tof.