震撼！模拟器的最后一块拼图！根据线圈参数计算磁感应强度分布曲线！

m24h

5个月5天前修改于 5个月5天前 IP:上海

操作

发表于《震撼！模拟器的最后一块拼图！根据线圈参数计算磁感应强度分布曲线！》

最近正在恶补电磁学和python的数据处理见猎心喜做了一个线圈磁场模拟器从最基本的电流元叉乘距离矢量的公式开始

写了个python程序其中wire()是线圈的形状可以改写成方的多层的将线圈分成多端直线进行计算函数并不复杂输入在线圈的长度位置对应的3d坐标就行

原本是用scipy的积分器(函数B1()) 后来发现电流元是对wire()的求导然后又积分这样速度慢很多于是优化成B2() 感觉这样精度更高速度相比略快

模拟是针对整个3维空间的但是使用中发现非常靠近线圈的地方因为距离平方的倒数会相对非常大在积分或者累加等计算时候容易有较大误差这点需要注意毕竟点如果落到线圈本体上都是个奇异点

代码(请根据自己需要修改参数和线圈函数) 和效果如下 (模拟计算是需要一些时间的不是死机)

当然这只是针对静态电流的

import numpy
from scipy import integrate
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))

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()
ax1.set_title('|B| of 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)

#draw X-Z->|B|
fig, ax2=pyplot.subplots(subplot_kw={'projection':'3d'}, figsize=pyplot.figaspect(0.5))
ax2.set_title('|B| of 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)
surf=ax2.plot_surface(x, z, b.reshape(simul_xz[1], simul_xz[0]),
                      cmap=cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)

fig, ax2b=pyplot.subplots()
ax2b.set_title('|B| of X-Z plane')
ax2b.set_xlabel('X')
ax2b.set_ylabel('Z')
ax2b.plot(coils.T[0], coils.T[2], alpha=0.5)
pcm=ax2b.pcolor(x, z, b.reshape(simul_xz[1], simul_xz[0]),
                cmap=cm.coolwarm, shading='auto')
fig.colorbar(pcm, extend='max')

#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()