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

文号933346

5305

gugugu2024/06/23原创磁阻研究 IP:上海

关键词

磁感应强度曲线计算python

在笔者之前写的文章《根据磁感应强度曲线来计算磁阻炮的弹丸受力》中成功计算了弹丸受力，可惜需要借用maxwell来计算磁感应强度分布曲线，这种方式由于需要maxwell的介入导致非常不方便，那么有没有一种办法可以根据线圈本身的参数直接计算磁感应强度分布曲线呢？

答案是肯定的，在github上的CoilGun.jl项目上笔者找到了这个问题的答案。计算公式如下：

不过由于笔者的电磁学能力很薄弱，所以无法跟着原文的思路进行推导，但在下面的代码验证中发现，公式应该要乘以0.5才成立。（CoilGun.jl项目的链接放着文末，有能力的大佬可以推导看看）。非常感谢托卡大佬的帮助，这个公式其实就是论文《两种载流螺线管轴线上磁场分布函数的比较》里的公式2的变形，并且确实少乘了个1/2 。

公式转化为python代码：

原文是的项目是CoilGun.jl，使用Julia代码写成的，可以用ai方便的转化成python代码。代码中的线圈参数与《根据磁感应强度曲线来计算磁阻炮的弹丸受力》中设置的一致即：线圈长度8，内径8.5，外径18，匝数100，激励电流200A。代码如下：

import math
from scipy.interpolate import interp1d  
import numpy as np
import pandas as pd  

from scipy.optimize import curve_fit  
import matplotlib.pyplot as plt  


class Coil:  
    def __init__(self, location, length, N, innerRadius, outerRadius):  
        self.length = length      #线圈长度
        self.innerRadius = innerRadius  #线圈内半径
        self.outerRadius = outerRadius  #线圈外半径
        self.N = N                      #线圈匝数
        self.μ0 = 4 * math.pi * 1e-7  # 真空磁导率   亨利每米 
  
    def BFieldCoil3(self, current, z):   #输入位置 计算磁感应强度！
        L=self.length
        a=L/2
        r2=self.outerRadius
        r1=self.innerRadius
        R=r2-r1
        def logarithm(pos):  
            return pos * math.log((math.sqrt(pos ** 2 + r2 ** 2) + r2) /  
                                  (math.sqrt(pos ** 2 + r1 ** 2) + r1))  
          
        return ( self.μ0*self.N*current/(2*L*R) ) * (logarithm(z+a)-logarithm(z-a)) 

coil = Coil(0, 8*1e-3, 100, 0.5*8.5*1e-3, 0.5*18*1e-3)  # 线圈位置0，线圈长度8 mm，线圈匝数100，线圈内半径8.5/2 mm，线圈外半径18/2 mm 
current = 200  # 电流为200A  

#计算-30至30 范围内的磁感应强度分布曲线
xx=[ pos for pos in np.arange(-30,30,0.1)]
yy=[ coil.BFieldCoil3( current, x*1e-3)   for x in xx]


#maxwell直接计算的磁感应强度分布曲线
Bx_csv = pd.read_csv('D:\pytohn_Sofe\dcp1\dcp1\mainSim\mainSim\Fb\\Bx_8_8.5_18_100_200.csv') 
Bxxx = np.array(Bx_csv.iloc[:, 3].values )   
Byyy = np.array(Bx_csv.iloc[:, 4].values )   


#画出曲线A 和 曲线B的曲线
plt.plot(xx, yy, color='yellow', label=' curve A')  
plt.scatter(Bxxx, Byyy, color='red', label=' maxwell B')  

plt.xlabel('x (mm)')  
plt.ylabel('B (T)')  
plt.legend()  
plt.show()

代码验证：

其中曲线A是代码计算的曲线，曲线B是由maxwell计算的曲线，可以看到这两条曲线的重合度非常高！之后笔者对其他参数的线圈也进行过多次验证，也表现出了极高的重合度！

结语：

现在笔者可以说是集齐了电磁炮模拟器的所有拼图了。HAcoilgun这个模拟器已经有十多年的历史了，只是非常可惜代码没有开源，我们没办法知道其是如何计算的。但是现在终于终于看到了曙光，新的模拟器即将要在kc诞生了！!!!

最后在这里我非常感谢各位大佬的帮助! 尤其是三水大佬，对我帮助非常多。 sticker

参考文献：

1: GitHub - laurium-labs/CoilGun.jl: Coil gun simulation to aid in coil gun design

2：根据磁感应强度曲线来计算磁阻炮的弹丸受力 - 科创网 (XXXXXXXXXXXX)

3: 线圈参数计算器（基于Wheeler公式计算电感量） - 科创网 (XXXXXXXXXXXX)

4: 你好，电磁炮（盒装）套件+书内购入口 - 科创网 (XXXXXXXXXXXX)

5: 行波加速器的仿真计算流程 - 科创网 (XXXXXXXXXXXX)

[修改于 5个月5天前 - 2024/06/24 18:19:39]

来自：物理高能技术 / 电磁炮

全部只看作者

已屏蔽原因：{{ notice.reason }}已屏蔽

屏蔽解除屏蔽编辑详情

~~空空如也

gugugu

作者

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

933347

1楼

在原项目里，或许是这里的公式2推导到公式3时候漏了个1/2 ？？？

引用

加载评论中，请稍候...

200字以内，仅用于支线交流，主线讨论请采用回复功能。

折叠评论

托卡马克

5个月6天前 IP:四川

933352

2楼

可以看看这篇论文

两种载流螺线管轴线上磁场分布函数的比较_肖建新.pdf 198.55KB PDF 112次下载预览

引用

加载评论中，请稍候...

200字以内，仅用于支线交流，主线讨论请采用回复功能。

折叠评论

gugugu

作者

5个月6天前 IP:上海

933353

3楼

引用托卡马克发表于2楼的内容

可以看看这篇论文

非常感谢！

那么破案了，CoilGun.jl项目中的公式其实是论文里的公式2。令x=z+a,又有L=2a。那么换元后，这俩公式就是一样了。而且CoilGun.jl项目中的公式确实少乘了1/ 2。

引用

加载评论中，请稍候...

200字以内，仅用于支线交流，主线讨论请采用回复功能。

折叠评论

三水合番

5个月6天前 IP:四川

933357

4楼

当年发的帖子都没人看了 sticker

XXXXXXXXXXXXXXXXXXXXXXXX/t/81277

引用

加载评论中，请稍候...

200字以内，仅用于支线交流，主线讨论请采用回复功能。

折叠评论

m24h

5个月6天前 IP:上海

933367

5楼

其实如果善用sympy scipy 可以做到更精确的模拟直接从安培定理开始积分也不难这正是Python的便利啊

引用

加载评论中，请稍候...

200字以内，仅用于支线交流，主线讨论请采用回复功能。

折叠评论

m24h

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

933478

6楼

最近正在恶补电磁学和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()

我不太懂电磁炮但是从模拟看出磁场最大和起伏最大的地方是靠近线圈壁的地方所以弹丸应该可以优化成空心的

另一个用途是在计算开关电源变压器磁芯时候可能会有一些指导意义吧

.........................

程序改了下增加打点显示防止误判死机增加一种新图可以看到线圈截面上的磁感应强度(有线圈阴影做对比)

引用

加载评论中，请稍候...

200字以内，仅用于支线交流，主线讨论请采用回复功能。

折叠评论

m24h

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

933517

7楼

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

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