目标是动画/模拟洒水器。基本思想是在一个实例中产生多个液滴，每个液滴以不同的角度离开喷头。 然而，我得到的不是一个移动的点，而是一条静态的直线。我已经尝试过只用一两滴/点的相同代码，但我仍然得到相同的直线。起初我以为这是因为我在动画中加入了拖动(在之前的试验中，当我加入拖动时，我得到了奇怪/意外的结果，所以我希望这是同样的问题)
下面是我写的代码。
编辑: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)
最佳答案
除了 @JohanC 指出的错误之外，您还使用表面积 (4 * pi * r^2) 而不是交叉水滴(近似于球体)的截面积 (pi * r^2)。此外，您似乎正在尝试使用假设真空的运动方程，当添加空气阻力或其他非常量力时，这些方程就会失效。
我建议您阅读 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()
这可能不会产生您想要的结果，但它可以正确模拟水滴的运动。
将 Drops 作为一个类实现并将它们存储在列表中，就像我一样，只需修改 update 函数即可轻松创建更多 Droplet 来模拟上述的许多“迭代”(为简单起见添加了 >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()
关于弹丸给出直线的Python动画，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/58944881/