弹丸给出直线的Python动画

Question

目标是动画/模拟洒水器。基本思想是在一个实例中创建多个液滴，每个液滴以不同的角度离开喷头。 然而，我得到的是一条直线，而不是一个移动的点。我尝试了相同的代码，只有一到两滴/点，但我仍然得到相同的直线。起初我以为是因为我在动画中加入了拖动（在之前的试验中，当我加入拖动时，我得到了奇怪/意外的结果，所以我希望这是同样的问题）

下面是我写的代码。

编辑：dragx 和 dragy 不再设置为相同的值。初始速度在更新函数之外计算。

每次指数增加时都会重新计算阻力，因为阻力取决于速度。

现在还包括导入。

我现在还尝试使用散点图而不是线，这在原点处给了我一个点。从那以后我又回到了这条线上

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.animation as animation
from math import sin, cos
import matplotlib.axes as axess

#pressure in system
def pressure(pump):
    return 10.197*pump*10**5/998.2

#Initial velocity based on the pressure of the system
def vinit(init_head):
    return (2*9.81*init_head)**0.5
# Number of droplets per iteration
n = 2
v0 = vinit(pressure(5.8))
cd = 0.5                                             #drag coefficient of a sphere;   
rho = 1.225                                          #density of air;                 kg/m^3

fps = 20
runtime = 1*20
numframes = fps*runtime
time = np.linspace(0,runtime,numframes)

#Initialize droplets

droplets = np.zeros(numframes, dtype = [('position', float, 2),
                                ('velocity', float, 2),
                                ('force', float, 2),
                                ('angle',float, 1),
                                ('radius', float, 1)])

droplets['radius'] = np.random.normal(0.0045, 0.001, numframes)
A = 4*np.pi*droplets['radius']**2
mass = ((4*np.pi*droplets['radius']**3)/3)*1000     #mass of waterdroplet
drag = np.zeros(numframes,dtype = [('dragx', float, 1),('dragy', float, 1)])

rads = np.radians(np.random.normal(37,1,numframes))
droplets['angle'] = rads
droplets['force'][:,1] = -9.8*mass

#Initialize velocity
for i in range(0, numframes):
    droplets['velocity'][i,0] = v0*cos(droplets['angle'][i])
    droplets['velocity'][i,1] = v0*sin(droplets['angle'][i])    

#Initialize the drag
drag['dragx'] = -0.5*rho*cd*A*droplets['velocity'][:,0]
drag['dragy'] = -0.5*rho*cd*A*droplets['velocity'][:,1]


droplets['position'][:,0] = 0


#Initialize the figure
fig,ax = plt.subplots()
axess.Axes.set_frame_on(ax,True)
ax.autoscale(True)
ax.set_xlim(0)
ax.set_ylim(0)

line = ax.plot(droplets['position'][0, 0], droplets['position'][0, 1],'b.')[0]
line = ax.plot(droplets['position'][:, 0], droplets['position'][:, 1],'b.')[0]
xdata = [droplets['position'][0, 0]]
ydata = [droplets['position'][0, 1]]
ln = plt.plot([],[],'b')[0]

def initfunc():
    droplets['position'][0,:] = 0
    ln.set_data(droplets['position'][0, 0], droplets['position'][0, 1])
    return ln

def update(framenum):
    index = framenum
    cd = 0.5                                          #drag coefficient of a sphere;   
    rho = 1.225                                       #density of air;                 kg/m^3
    A = 4*np.pi*droplets['radius']**2                 #surface area of droplet;        m^2
    mass = ((4*np.pi*droplets['radius']**3)/3)*1000   #mass of waterdroplet

    #Update the drag force on the droplet
    drag['dragx'] = -0.5*rho*cd*A*droplets['velocity'][:,0]
    drag['dragy'] = -0.5*rho*cd*A*droplets['velocity'][:,1]

    droplets['force'][index,0] += drag['dragx'][index]
    droplets['force'][index,1] += drag['dragy'][index]
    #droplets['position'][0] = [0,0]
    droplets['velocity'][index,0] = droplets['velocity'][index,0] + (droplets['force'][index,0]/mass[index])*index
    droplets['velocity'][index,1] = droplets['velocity'][index,1] + (droplets['force'][index,1]/mass[index])*index
    droplets['position'][index,0] = droplets['position'][index,0] + (droplets['velocity'][index,0])*index  
    droplets['position'][index,1] = droplets['position'][index,1] + (droplets['velocity'][index,1])*index  

    xdata.append(droplets['position'][index,0])
    ydata.append(droplets['position'][index,1])
    ln.set_data(xdata,ydata)

    line.set_data(droplets['position'][:,0], droplets['position'][:,1])

    return ln



sprink = animation.FuncAnimation(fig, update,init_func = initfunc,interval= 200,  frames = numframes)
plt.show()

sprink.save('Sprinkler.mp4', fps = fps)

Answer 1

除了@JohanC 指出的错误，您使用的是表面积 (4 * pi * r^2) 而不是横截面积 (@987654329 @) 的水滴（近似为球体）。此外，您似乎正在尝试使用假设真空的运动方程，当添加空气阻力或其他非恒定力时，这些方程会失效。

我建议您阅读this answer 以获得更完整的解释来计算带空气阻力的弹道运动，但这里有一个简化版本。

计算粒子在介质中移动的位置本质上是求解first order ODE，有很多方法可以数值求解这样的方程。最简单（尽管最不稳健）的方法是将问题分解为小段，并近似从当前点到下一个点的变化。对于这个问题，假设您使用足够小的时间步长，这样的方法效果很好。

在你的情况下，我们可以做这样的事情

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
from math import sin, cos

g = 9.81
rho = 1.225
v0 = 50.0    # This can be initialized as desired, I chose 50 m/s for demonstration
C = 0.5

nframes = 100
ndrops = 10
tend = 8.0
time = np.linspace(0.1, tend, nframes)   # Runtime of 4s chosen arbitrarily
dt = time[1] - time[0]                    # change in time between steps
maxs = [0.0, 0.0]
r = [0.001, 0.0045]                       # This stores the range of radii to be used

# I have implemented a class to store the drops for brevity
class drop:

    def __init__(self, pos, vel, r):
        self.pos = pos
        self.vel = vel
        self.r = r

class sprinkler:
# I implemented this as a class to simplify updating the scatter plot
    def __init__(self):
        self.fig, self.ax = plt.subplots()
        self.drops = [None]*ndrops
        self.maxt = 0.0
        ## This step is simply to figure out how far the water drops can travel
        ## in order to set the bounds of the plot accordingly
        angles = np.linspace(0, np.pi/2, 90)    # We only need to sample from 0 to pi/2
        for angle in angles:
            m = [drop([0.0, 0.1], [v0*cos(angle), v0*sin(angle)], 0.001),
                 drop([0.0, 0.1], [v0*cos(angle), v0*sin(angle)], 0.0045)]
            for d in m:
                t = 0.0
                coef = - 0.5*C*np.pi*d.r**2*rho
                mass = 4/3*np.pi*d.r**3*1000
                while d.pos[1] > 0:
                    a = np.power(d.vel, 2) * coef * np.sign(d.vel)/mass
                    a[1] -= g
                    d.pos += (d.vel + a * dt) * dt
                    d.vel += a * dt
                    t += dt
                    if d.pos[1] > maxs[1]:
                        maxs[1] = d.pos[1]
                    if d.pos[0] > maxs[0]:
                        maxs[0] = d.pos[0]
                    if d.pos[1] < 0.0:
                        if t > self.maxt:
                            self.maxt = t
                        break
        for ii in range(ndrops):    # Make some randomly distributed water drops
            self.drops[ii] = drop([0.0, 0.0], [cos(np.random.random()*np.pi)*v0,
                             sin(np.random.random()*np.pi)*v0],
                             np.random.random()*(r[1]-r[0]) + r[0])

        anim = animation.FuncAnimation(self.fig, self.update, init_func=self.setup,
                                       interval=200, frames=nframes)
        anim.save('Sprinkler.gif', fps = 20, writer='imagemagick')

    def setup(self):
        self.scat = self.ax.scatter([d.pos[0] for d in self.drops],
                               [d.pos[1] for d in self.drops], marker='.', color='k')
        self.ax.set_xlim([-maxs[0]*1.15, maxs[0]*1.15])
        self.ax.set_ylim([0, maxs[1]*1.15])

        return self.scat,

    def update(self,frame):  # Use set_offsets to move the water drops
        self.step()          # Advance to the next 'step'
        self.scat.set_offsets(np.stack([[d.pos[0] for d in self.drops],
                                       [d.pos[1] for d in self.drops]], 1))

        return self.scat,

    def step(self):
        for ii in range(ndrops):
            coef = - 0.5*C*np.pi*self.drops[ii].r**2*rho    # Aggregated coefficient
            mass = 4/3*np.pi*self.drops[ii].r**3*1000
            a = np.power(self.drops[ii].vel, 2)*coef * np.sign(self.drops[ii].vel)/mass
            a[1] -= g
            # Approximate how much the position and velocity would change if we assume
            # a(t) does not change between t and t+dt
            self.drops[ii].pos += np.array(self.drops[ii].vel) * dt + 0.5 * a * dt**2
            self.drops[ii].vel += a * dt
            if self.drops[ii].pos[1] < 0.0:     # Check if the drop has hit the ground
                self.drops[ii].pos[1] = 0.0
                self.drops[ii].vel = [0.0, 0.0]

sprinkler()

这可能无法准确地产生您所追求的效果，但它正确地模拟了水滴的运动。

通过简单地修改update函数（create为简单起见添加方法):

    def update(self,frame):  # Use set_offsets to move the water drops
        # Create between 0 and ndrops new drops, but only until tend-maxt
        if time[frame] < tend - self.maxt*1.1:
            self.create(np.random.randint(0, ndrops))
        self.step()          # Advance to the next 'step'
        self.scat.set_offsets(np.stack([[d.pos[0] for d in self.drops],
                                        [d.pos[1] for d in self.drops]], 1))

        return self.scat,

    def create(self, i):
        for ii in range(i):
            self.drops.append(drop([0.0, 0.0], [cos(np.random.random()*np.pi)*v0,
                             sin(np.random.random()*np.pi)*v0],
                             np.random.random()*(r[1]-r[0]) + r[0]))

编辑 - 生成第二个动画的整个代码

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
from math import sin, cos

g = 9.81
rho = 1.225
v0 = 50.0    # This can be initialized as desired, I chose 50 m/s for demonstration
C = 0.5

nframes = 100
ndrops = 10
tend = 8.0
time = np.linspace(0.1, tend, nframes)   # Runtime of 4s chosen arbitrarily
dt = time[1] - time[0]                    # change in time between steps
maxs = [0.0, 0.0]
r = [0.001, 0.0045]                       # This stores the range of radii to be used

# I have implemented a class to store the drops for brevity
class drop:

    def __init__(self, pos, vel, r):
        self.pos = pos
        self.vel = vel
        self.r = r

class sprinkler:
# I implemented this as a class to simplify updating the scatter plot
    def __init__(self):
        self.fig, self.ax = plt.subplots()
        self.drops = [None]*ndrops
        self.maxt = 0.0
        ## This step is simply to figure out how far the water drops can travel
        ## in order to set the bounds of the plot accordingly
        angles = np.linspace(0, np.pi/2, 90)    # We only need to sample from 0 to pi/2
        for angle in angles:
            m = [drop([0.0, 0.1], [v0*cos(angle), v0*sin(angle)], 0.001),
                 drop([0.0, 0.1], [v0*cos(angle), v0*sin(angle)], 0.0045)]
            for d in m:
                t = 0.0
                coef = - 0.5*C*np.pi*d.r**2*rho
                mass = 4/3*np.pi*d.r**3*1000
                while d.pos[1] > 0:
                    a = np.power(d.vel, 2) * coef * np.sign(d.vel)/mass
                    a[1] -= g
                    d.pos += (d.vel + a * dt) * dt
                    d.vel += a * dt
                    t += dt
                    if d.pos[1] > maxs[1]:
                        maxs[1] = d.pos[1]
                    if d.pos[0] > maxs[0]:
                        maxs[0] = d.pos[0]
                    if d.pos[1] < 0.0:
                        if t > self.maxt:
                            self.maxt = t
                        break
        for ii in range(ndrops):    # Make some randomly distributed water drops
            self.drops[ii] = drop([0.0, 0.0], [cos(np.random.random()*np.pi)*v0,
                             sin(np.random.random()*np.pi)*v0],
                             np.random.random()*(r[1]-r[0]) + r[0])

        anim = animation.FuncAnimation(self.fig, self.update, init_func=self.setup,
                                       interval=200, frames=nframes)
        anim.save('Sprinkler.gif', fps = 20, writer='imagemagick')

    def setup(self):
        self.scat = self.ax.scatter([d.pos[0] for d in self.drops],
                                    [d.pos[1] for d in self.drops], marker='.', color='k')
        self.ax.set_xlim([-maxs[0]*1.15, maxs[0]*1.15])
        self.ax.set_ylim([0, maxs[1]*1.15])

        return self.scat,

    def update(self,frame):  # Use set_offsets to move the water drops
        # Create between 0 and ndrops new drops, but only until tend-maxt
        if time[frame] < tend - self.maxt*1.1:
            self.create(np.random.randint(0, ndrops))
        self.step()          # Advance to the next 'step'
        self.scat.set_offsets(np.stack([[d.pos[0] for d in self.drops],
                                        [d.pos[1] for d in self.drops]], 1))

        return self.scat,

    def create(self, i):
        for ii in range(i):
            self.drops.append(drop([0.0, 0.0], [cos(np.random.random()*np.pi)*v0,
                             sin(np.random.random()*np.pi)*v0],
                             np.random.random()*(r[1]-r[0]) + r[0]))
    def step(self):
        for ii in range(len(self.drops)):
            coef = - 0.5*C*np.pi*self.drops[ii].r**2*rho    # Aggregated coefficient
            mass = 4/3*np.pi*self.drops[ii].r**3*1000
            a = np.power(self.drops[ii].vel, 2) * coef * np.sign(self.drops[ii].vel)/mass
            a[1] -= g
            # Approximate how much the position and velocity would change if we assume
            # a(t) does not change between t and t+dt
            self.drops[ii].pos += np.array(self.drops[ii].vel) * dt + 0.5 * a * dt**2
            self.drops[ii].vel += a * dt
            if self.drops[ii].pos[1] < 0.0:     # Check if the drop has hit the ground
                self.drops[ii].pos[1] = 0.0
                self.drops[ii].vel = [0.0, 0.0]

sprinkler()