增加了梯度的模拟 也许对受力分析有用

import numpy
from scipy import integrate, optimize
from matplotlib import pyplot, cm

numpy.seterr(all='warn')

#define of electromagnetic parameters
em_u=1 # times of u0
em_I=200 # Amps, current from wire_range[0] to wire_range[1]

#define of simulation cube 
simul_cube=[[-0.01, 0.01], [-0.01,0.01], [-0.01,0.01]] # coordinate range of simulation
simul_mf=[6, 6, 6] # number of simulating points of X/Y/Z axises
simul_xz=[30,30]
simul_z=200

#define of conducting wire
wire_range=[-10., 10.] # range of `l` for wire()
wire_segments=2000
      
def wire(l):   # wire shape
  r=0.007
  n=100
  h=0.008
  a=2*numpy.pi*n/(wire_range[1]-wire_range[0])
  b=h/(wire_range[1]-wire_range[0])
  return numpy.array([numpy.cos(a*l)*r, 
                      numpy.sin(a*l)*r, 
                      b*l])

coils=numpy.vectorize(wire, signature='()->(n)')(numpy.linspace(*wire_range, wire_segments))
coils_y0l=optimize.root(lambda l:wire(l)[1], numpy.linspace(*wire_range, 1000)).x
coils_y0l=numpy.extract(coils_y0l>=wire_range[0], coils_y0l)
coils_y0l=numpy.extract(coils_y0l<=wire_range[1], coils_y0l)
coils_xz=wire(coils_y0l).take([0,2],0)

def dB1(l, p):
  dl=(wire_range[1]-wire_range[0])/wire_segments
  w=(wire(l+dl/2)-wire(l-dl/2))/dl
  r=p-wire(l)
  return numpy.cross(w, r)/numpy.linalg.norm(r)**3

def B1(p):
  print('.', end='', flush=True)
  return 1e-7*em_I*integrate.quad_vec(dB1, *wire_range, args=(p,))[0]

def B2(p):
  print('.', end='', flush=True)
  ret=numpy.array([0.,0.,0.])
  a=p-coils[0]
  for t in coils[1:]:
    b=p-t
    ret+=numpy.cross(a, b)/numpy.linalg.norm((a+b)/2)**3
    a=b
  return 1e-7*em_I*ret

B=numpy.vectorize(B2, signature='(n)->(n)')

#draw Z->|B|
fig, ax1=pyplot.subplots()
fig.subplots_adjust(right=0.85)
ax1.set_title('|B| on Z axis')
ax1.set_xlim(*simul_cube[2])
ax1.set_xlabel('Z')
ax1.set_ylabel('|B|')
zs=numpy.linspace(*simul_cube[2], simul_z)
xyz=numpy.r_[numpy.zeros_like(zs), numpy.zeros_like(zs),zs].reshape(3,-1)
bz=numpy.linalg.norm(B(xyz.T), axis=1)
ax1.plot(zs, bz, color='black', label='|B|')
grad=numpy.gradient(bz)
ax1b=ax1.twinx()
ax1b.set_ylabel('Grad', color='blue')
ax1b.plot(zs, grad, color='blue', label='Gradient')

#draw X-Z->|B|
fig, ax2=pyplot.subplots(subplot_kw={'projection':'3d'}, figsize=pyplot.figaspect(0.5))
ax2.set_title('|B| in X-Z plane')
ax2.set_xlim(*simul_cube[0])
ax2.set_ylim(*simul_cube[2])
ax2.set_xlabel('X')
ax2.set_ylabel('Z')
ax2.set_zlabel('|B|')
x,z=numpy.meshgrid(numpy.linspace(*simul_cube[0], simul_xz[0]),
                   numpy.linspace(*simul_cube[2], simul_xz[1]))
b=numpy.linalg.norm(B(numpy.r_[x,numpy.zeros_like(x),z].reshape(3,-1).T), axis=1).reshape(simul_xz[1], simul_xz[0])
surf=ax2.plot_surface(x, z, b,
                      cmap=cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)

fig, ax2b=pyplot.subplots()
ax2b.set_title('|B| in X-Z plane')
ax2b.set_xlabel('X')
ax2b.set_ylabel('Z')
pcm=ax2b.pcolor(x, z, b,
                cmap=cm.coolwarm, shading='auto')
fig.colorbar(pcm, extend='max')
ax2b.scatter(*coils_xz, s=2, c='k', alpha=0.2)

#draw gradient
grad_x=numpy.gradient(b, axis=1)
grad_z=numpy.gradient(b, axis=0)
grad=numpy.sqrt(grad_x**2+grad_z**2)
fig, ax2c=pyplot.subplots()
ax2c.set_title('gradient of |B| in X-Z plane')
ax2c.set_xlabel('X')
ax2c.set_ylabel('Z')
pcm=ax2c.pcolor(x, z, grad,
                cmap=cm.coolwarm, shading='auto')
fig.colorbar(pcm, extend='max')
ax2c.scatter(*coils_xz, s=2, c='k', alpha=0.2)
ax2c.quiver(x.reshape(-1), z.reshape(-1), grad_x.reshape(-1), grad_z.reshape(-1))

#draw magnetic field
fig, ax3=pyplot.subplots(subplot_kw={'projection':'3d'}, figsize=pyplot.figaspect(0.5))
ax3.set_title('Magnetic field')
ax3.set_xlim(*simul_cube[0])
ax3.set_ylim(*simul_cube[1])
ax3.set_zlim(*simul_cube[2])
ax3.set_xlabel('X')
ax3.set_ylabel('Y')
ax3.set_zlabel('Z')
ax3.plot(*coils.T, alpha=0.5)
xyz=numpy.r_[numpy.meshgrid(numpy.linspace(*simul_cube[0], simul_mf[0]),
                            numpy.linspace(*simul_cube[1], simul_mf[1]),
                            numpy.linspace(*simul_cube[2], simul_mf[2]))
            ].reshape(3,-1)
mf=B(xyz.T)
ax3.quiver(*xyz, *mf.T, length=numpy.abs(simul_cube[0][1]-simul_cube[0][0])/simul_mf[0], normalize=True)

#show all
pyplot.show()

感觉相比于线圈中心部分 靠近线圈壁的受力要大很多 从加速的角度看 空心柱形弹丸会非常有最终速度上的优势

------

忽然想起 上图颜色对应的是梯度的模 如果改成Z分量 应该更合适 我这里就不改了 需要的可以注意一下