时间戳

首页 > 解决方案 > 如何绘制弹丸在重力、浮力和空气阻力作用下的运动？

python - 如何绘制弹丸在重力、浮力和空气阻力作用下的运动？

问题描述

我正在尝试绘制在重力、浮力和阻力作用下的质量的抛射运动图。基本上，我想在绘图上显示浮力和阻力对飞行距离、飞行时间和速度变化的影响。

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. 重力、空气阻力和浮力作用下的质量

标签： pythonnumpymatplotlibphysics

解决方案

您必须根据给定时间的力总和计算每个位置。为此，最好从随时计算净力开始，然后用它来计算加速度、速度和位置。对于下面的计算，假设浮力和重力是恒定的（这在现实中是不正确的，但在这种情况下它们的可变性的影响可以忽略不计），还假设初始位置是(0,0)尽管这可以简单地更改为任何初始位置。

F_x = tau_x 
F_y = tau_y + bouyancy + gravity

其中tau_x和分别是和方向tau_y的阻力。速度和由下式给出xyv_xv_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进行评估。dtdt * n = tn

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。

推荐阅读