Python中带阻力球形粒子下落：如何生成垂直速度Vy随时间变化列表及绘图

Got it, let's work through this problem of calculating the vertical velocity of a falling spherical particle with drag over time. We'll use Euler's method to discretize the given differential equation, since you provided the update formula dVy = -g*dt - (b/m)*Vy*dt — quick note: I'll adjust the sign convention so positive velocity points downward (since the particle is falling) to make the results more intuitive, but the core logic matches your equation.

We'll write a Python script that:

• Takes user input for mass (m), gravitational acceleration (g), drag coefficient constant (b), time step (∆t), and maximum time (t_max)
• Calculates velocity at each time step using the provided update rule
• Stores time and velocity values in lists for easy access
• Plots the Vy-time curve to visualize the behavior

Complete Code #

import matplotlib.pyplot as plt

# Get user input for parameters
m = float(input("Enter mass of the particle (kg): "))
g = float(input("Enter gravitational acceleration (m/s²): "))
b = float(input("Enter drag coefficient constant (b): "))
dt = float(input("Enter time step (∆t, s): "))
t_max = float(input("Enter maximum time (t_max, s): "))

# Initialize variables
t = 0.0
vy = 0.0  # Starting velocity at release (t=0)
time_list = [t]
vy_list = [vy]

# Calculate velocity at each time step
while t < t_max:
    # Update velocity using the given discrete equation
    # Adjusted sign: positive Vy is downward, so dv = g*dt - (b/m)*vy*dt
    dv = g * dt - (b / m) * vy * dt
    vy += dv
    t += dt
    
    # Append to lists
    time_list.append(t)
    vy_list.append(vy)

# Print the Vy-time values
print("\nTime (s) | Vertical Velocity (m/s)")
print("-----------------------------------")
for t_val, vy_val in zip(time_list, vy_list):
    print(f"{t_val:.2f}      | {vy_val:.4f}")

# Plot Vy vs Time
plt.figure(figsize=(10, 6))
plt.plot(time_list, vy_list, marker='o', linestyle='-', color='b')
plt.xlabel('Time (s)')
plt.ylabel('Vertical Velocity (m/s)')
plt.title('Vertical Velocity of Falling Particle with Drag Over Time')
plt.grid(True)
plt.show()

How It Works #

• Input Handling: The script prompts you to enter all required parameters directly, so you can test different values without editing the code itself.
• Euler's Method: We iterate through each time step, updating velocity using dv = g*dt - (b/m)*vy*dt — this rearranges your original equation and flips the sign to make downward velocity positive (swap g*dt to -g*dt if you prefer upward as the positive direction).
• Data Storage: Every time and velocity value is stored in lists, which we print in a clean table format and use for plotting.
• Visualization: The plot shows velocity starting at 0, increasing, and eventually leveling off at terminal velocity (when drag balances gravity — you can verify this analytically with v_terminal = (m*g)/b).

Example Output #

If you input:

• m = 0.1 kg
• g = 9.8 m/s²
• b = 0.05
• dt = 0.1 s
• t_max = 10 s

You'll see velocity approach ~19.6 m/s (the calculated terminal velocity), with the plot showing a curve that flattens as time passes.

内容的提问来源于stack exchange，提问作者Aran G