![0I@C1SG_J$7Y6_9F6$`7]HD.png](https://img.kechuang.org:81/r/355045?c=resource)
在原项目里，或许是这里的公式2推导到公式3时候漏了个1/2 ？？？
90268
%7B%22isLastPage%22%3Atrue%2C%22notes%22%3A%5B%5D%2C%22pid%22%3A%22933346%22%2C%22tid%22%3A%2290268%22%2C%22mainForumsId%22%3A%5B%22367%22%5D%2C%22categoriesId%22%3A%5B%5D%2C%22tcId%22%3A%5B%5D%7D
%7B%22isEditMode%22%3Afalse%7D
震撼！模拟器的最后一块拼图！ 根据线圈参数计算 磁感应强度分布曲线！
关键词
磁感应强度曲线计算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()代码验证：
![03Y3EZ}BE0TR2OHWU{DI])2.png](https://img.kechuang.org:81/r/355044?c=resource)
其中曲线A是代码计算的曲线，曲线B是由maxwell计算的曲线，可以看到这两条曲线的重合度非常高！ 之后笔者对其他参数的线圈也进行过多次验证，也表现出了极高的重合度！
结语：
现在笔者可以说是集齐了电磁炮模拟器的所有拼图了。HAcoilgun这个模拟器已经有十多年的历史了，只是非常可惜代码没有开源，我们没办法知道其是如何计算的。 但是现在终于终于看到了曙光，新的模拟器即将要在kc诞生了！!!!
最后在这里我非常感谢各位大佬的帮助! 尤其是三水大佬，对我帮助非常多。
参考文献：
1: GitHub - laurium-labs/CoilGun.jl: Coil gun simulation to aid in coil gun design
2：根据磁感应强度曲线来计算磁阻炮的弹丸受力 - 科创网 (XXXXXXXXXXXX)
3: 线圈参数计算器（基于Wheeler公式计算电感量） - 科创网 (XXXXXXXXXXXX)
[修改于 5个月3天前 - 2024/06/24 18:19:39]
可以看看这篇论文
引用
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。
引用
1
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。
其实如果善用sympy scipy 可以做到更精确的模拟 直接从安培定理开始积分也不难 这正是Python的便利啊
引用
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。
最近正在恶补电磁学和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字以内，仅用于支线交流，主线讨论请采用回复功能。
增加了梯度的模拟 也许对受力分析有用
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分量 应该更合适 我这里就不改了 需要的可以注意一下
引用
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。