リベンジ・・・磁場分布シミュレータv0.2をpythonで

v0.2です。 v0.1はちょっと間違っていたので、直しました。

電流のループに沿って周回積分をするようにしました。現実のコイルに即して。

さて、直したことで劇的によくなったかというと・・・？

ぱっと見は同じ・・・。元もそんなに悪くない。

でもまあ、よく見ると、コイルの外側の磁場がちゃんと消えるようになっているなどの変化があります。より正しそうですね。

軸上の磁場分布をプロットしてみるともう少し明らかです。

ちゃんと計算した方では、スプリットギャップ内での磁場の落ちが急峻です。

スプリットはやっぱり厳しいんですね。

元のプログラムではその辺が楽観的過ぎるので、そこはよくなかった。

というわけで新しくしたことは意味があったと言うことです。勉強になりました。

たくさん積分しているけれど、コードは割と短いです。

このプログラムは、巻き数とレイヤー数しか指定できないですが、

改良すれば、レイヤーごとに任意の巻き数での磁場分布もシミュレートできるとおもいます。やってみてね！

ちなみに、長さの次元はワイヤーの直径でスケールされています。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


"""
Magnetic field distribution simulator v0.2

A. Ikeda of UEC Tokyo
5 May 2022

References
[1]https://krajit.github.io/sympy/vectorFields/vectorFields.html
[2]Yosuke Nagaoka, Butsuri Nyumon 3, Iwanami Shoten

"""


### basic parameters of the split coil ###
bore = 20
wind = 10
layer = 10
gap = 20 # gap = 0 for mono coil

hgap = gap / 2
hbore = bore / 2



### basic parameters for calculations ###
ran = 30 # range of z-y plane to calculate
div = 30 # how many points to divide the plane
mg = np.linspace(-ran, ran, div)
z, y = np.meshgrid(mg, mg) # B field will be calculated on z-y plane



### vector definition ###

# The coil axis is directing in z direction. 
# We calculate B fields on the y-z plane.

# vector from wire positions to a field position, so to speak, (r - r')
def wpxx(rw, q, x0):
    return rw * np.cos(q) - x0
def wpyy(rw, q, y0):
    return rw * np.sin(q) - y0
def wpzz(zw, q, z0):
    return zw - z0

# norm of the vector above, so to speak, |r - r'|
def norm(rw, zw, q, x0, y0, z0):
    return np.sqrt(wpxx(rw, q, x0) ** 2 + wpyy(rw, q, y0) ** 2 + wpzz(zw, q, z0) ** 2)

# vector of current, so to speak, i(r')
def ix(i0, q):
    return -i0 * np.sin(q)
def iy(i0, q):
    return i0 * np.cos(q)
def iz(i0, q):
    return 0

### vector definition ###


### Biot-Savart's law for calculating B_z and B_y by a tiny part of magnet wire ###

def bz_wireat(rw, zw, q, x0, y0, z0, i0):
    return norm(rw, zw, q, x0, y0, z0) ** -3 * (ix(i0, q) * wpyy(rw, q, y0) - iy(i0, q) * wpxx(rw, q, x0))

def by_wireat(rw, zw, q, x0, y0, z0, i0):
    return norm(rw, zw, q, x0, y0, z0) ** -3 * (iz(i0, q) * wpxx(rw, q, x0) - ix(i0, q) * wpzz(zw, q, z0))



# division of 2pi integration of a coil
qqdiv = 36
qdiv = 2 * np.pi / qqdiv
qq = 0
u = z - z
v = y - y


### calculation of B_z and B_y on the plane (z, y, x = 0) 
### with the integration in 2pi of a current loop
### and further integration for winding and layers of the coil

for rr in np.arange(layer):
    for zz in np.arange(wind):
        for qq in np.arange(0, 2 * np.pi, qdiv):
            u +=  bz_wireat(rw = hbore + rr, zw = zz + hgap, q = qq, x0 = 0, y0 = y, z0 = z, i0 = 1) * qdiv
            v +=  by_wireat(rw = hbore + rr, zw = zz + hgap, q = qq, x0 = 0, y0 = y, z0 = z, i0 = 1) * qdiv

for rr in np.arange(layer):
    for zz in np.arange(wind):
        for qq in np.arange(0, 2 * np.pi, qdiv):
            u +=  bz_wireat(rw = hbore + rr, zw = -(zz + hgap), q = qq, x0 = 0, y0 = y, z0 = z, i0 = 1) * qdiv
            v +=  by_wireat(rw = hbore + rr, zw = -(zz + hgap), q = qq, x0 = 0, y0 = y, z0 = z, i0 = 1) * qdiv
            
            
            
plt.figure(figsize = (10,10))
plt.quiver(z,y,u,v) 
plt.savefig('save0.png')
plt.show()

plt.plot(z[0], np.abs(u[round(div/2)]))
np.savetxt('x0.csv', z[0])
np.savetxt('y0.csv', np.abs(u[round(div/2)]))
plt.savefig('save1.png')
plt.show()