コイルの作る磁場シミュレーション

はじめに

電流の作り出す磁場はビオザバールの法則で計算できる。

H= \frac{\mu_0}{4 \pi} \frac{I dl \times \vec{r}}{| \vec{r}^3 |}

円形コイルの中心磁場は以下のように計算できる

H = \frac{I}{2r}

円中心以外はどうなっているんだろう？
これがテーマです。

円形コイル中心の磁場計算：シミュレーション法の検証

表題の内容を計算してみる。

import math

x  = 0
y = 0


def culc_H(s_x,s_y,e_x,e_y): #電流の向きでstart-endを決めている
    n = 100
    
    #print("("+str(s_x) +"," + str(s_y) + ")" +"->"+ "("+str(e_x) +"," + str(e_y) + ")")
    
    #------計算パート-------
    #l = math.sqrt((s_x - e_x )**2 + (s_y - e_y))
    #dl = l/n
    dl_x = (e_x - s_x)/n
    dl_y = (e_y - s_y)/n
    
    H = 0
    
    for i in range(n):
        r_x = (s_x + dl_x*i) - x
        r_y = (s_y + dl_y*i) - y
        r = math.sqrt(r_x**2 + r_y**2)
        
        if r == 0:
            H = 0
            break
        
        # dl×r/|r|^3  (rの取り方が電流素から計測点への距離だったため、符号を逆にする)
        H = H - (dl_x*r_y - dl_y*r_x)/(r**3)
        #print(r_y)
        
    H = H /4 /math.pi
    return H


# n角多角形

n = 100
r = 1

dtheta = 2*math.pi/n

theta = 0
s_x = r
s_y = 0

H = 0
for i in range(n):
    theta = theta + dtheta
    e_x = r* math.cos(theta)
    e_y = r* math.sin(theta)
    dH = culc_H(s_x,s_y,e_x,e_y)
    H = H + dH
    print(dH)
    
    s_x = e_x
    s_y = e_y


print(H)

結果

\begin{eqnarray}
n = 100, r = 1 → H(0,0)/ \mu_0 = 0.500164\\
n = 100, r = 5 → H(0,0)/ \mu_0  = 0.100032\\
n = 10, r = 10 → H(0,0)/ \mu_0  = 0.517121\\
\end{eqnarray}

理論値と概ね一致している。

中心以外の磁場

import math
import numpy as np
import matplotlib.pyplot as plt

def culc_H(s_x,s_y,e_x,e_y,x,y): #電流の向きでstart-endを決めている
    n = 100
    
    #print("("+str(s_x) +"," + str(s_y) + ")" +"->"+ "("+str(e_x) +"," + str(e_y) + ")")
    
    #------計算パート-------
    #l = math.sqrt((s_x - e_x )**2 + (s_y - e_y))
    #dl = l/n
    dl_x = (e_x - s_x)/n
    dl_y = (e_y - s_y)/n
    
    H = 0
    
    for i in range(n):
        r_x = (s_x + dl_x*i) - x
        r_y = (s_y + dl_y*i) - y
        r = math.sqrt(r_x**2 + r_y**2)
        
        if abs(r) < 0.01:
            H = 0
            break
        
        # dl×r/|r|^3  (rの取り方が電流素から計測点への距離だったため、符号を逆にする)
        H = H - (dl_x*r_y - dl_y*r_x)/(r**3)
        #print(r_y)
        
    H = H /4 /math.pi
    return H

# n角多角形
def culc_H_xy(x,y):
    
    n = 20
    r = 1
    
    dtheta = 2*math.pi/n
    
    theta = 0
    s_x = r
    s_y = 0
    
    H = 0
    for i in range(n):
        theta = theta + dtheta
        e_x = r* math.cos(theta)
        e_y = r* math.sin(theta)
        dH = culc_H(s_x,s_y,e_x,e_y,x,y)
        H = H + dH
        #print(dH)
        
        s_x = e_x
        s_y = e_y
        
    return H

H = culc_H_xy(0,0)

# グラフ
n = 50
a = 1.2
x = np.linspace(-1*a, a, n)
y = np.linspace(-1*a, a, n)
z = np.zeros((n,n))

for x_i in range(n):
    for y_i in range(n):
        z[y_i,x_i] = culc_H_xy(x[x_i],y[y_i])


# プロットの作成
plt.figure(figsize=(8, 6))
x,y = np.meshgrid(x,y)

sc = plt.scatter(x, y, c=z, cmap='viridis',s=100)
plt.colorbar(sc)
plt.xlabel('X axis')
plt.ylabel('Y axis')
plt.title('Z value color map on XY plane using scatter')
plt.show()



fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111,projection='3d')
#ax.set_aspect('equal')
ax.scatter(x,y,z)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

r=1の円に対してシミュレーションを実施。

ビオザバールの式で$1/r^3$があることから、$r=0$付近の計算は値が不安定になる。

円環付近は値が大きくなるが、円環の外はほぼ$H=0$、円環内部は$H=0.5$と読める。
中心からの距離と、磁場$H$をグラフにすると以下になる。

$r=0.5$で中心磁場$H=0.5$とほぼ同じ値になっている。

$r = 0 ～0.8$をグラフにすると以下になる。

四角形の磁場

import math
import numpy as np
import matplotlib.pyplot as plt

def culc_H(s_x,s_y,e_x,e_y,x,y): #電流の向きでstart-endを決めている
    n = 100
    
    #print("("+str(s_x) +"," + str(s_y) + ")" +"->"+ "("+str(e_x) +"," + str(e_y) + ")")
    
    #------計算パート-------
    #l = math.sqrt((s_x - e_x )**2 + (s_y - e_y))
    #dl = l/n
    dl_x = (e_x - s_x)/n
    dl_y = (e_y - s_y)/n
    
    H = 0
    
    for i in range(n):
        r_x = (s_x + dl_x*i) - x
        r_y = (s_y + dl_y*i) - y
        r = math.sqrt(r_x**2 + r_y**2)
        
        if abs(r) < 0.001:
            H = 0
            break
        
        # dl×r/|r|^3  (rの取り方が電流素から計測点への距離だったため、符号を逆にする)
        H = H - (dl_x*r_y - dl_y*r_x)/(r**3)
        #print(r_y)
        
    H = H /4 /math.pi
    return H

# 四角形
def culc_H_xy(x,y):
    
    # 四角形　a*bの場合
    a = 1
    b = 1
    
    H = 0
    
    H = H + culc_H(-1*a/2,b/2,-1*a/2,-1*b/2,x,y)  # 左辺
    H = H + culc_H(-1*a/2,-1*b/2,a/2,-1*b/2,x,y)  # 下辺
    H = H + culc_H(a/2,-1*b/2,a/2,b/2,x,y)        # 右辺
    H = H + culc_H(a/2,b/2,-1*a/2,b/2,x,y)        # 上辺        
    
    return H


# グラフ
n = 50
a = 1
x = np.linspace(-1*a, a, n)
y = np.linspace(-1*a, a, n)
z = np.zeros((n,n))

for x_i in range(n):
    for y_i in range(n):
        z[y_i,x_i] = culc_H_xy(x[x_i],y[y_i])


# プロットの作成
plt.figure(figsize=(8, 6))
x,y = np.meshgrid(x,y)

sc = plt.scatter(x, y, c=z, cmap='viridis',s=100)
plt.colorbar(sc)
plt.xlabel('X axis')
plt.ylabel('Y axis')
plt.title('Z value color map on XY plane using scatter')
plt.show()

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111,projection='3d')
#ax.set_aspect('equal')
ax.scatter(x,y,z)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

$1 \times 1$の長方形のコイルに対し実施した。相変わらず、コイル付近の値は値が大きく出てしまっている。

四角形の内部は約$H=1.0$に落ち着いている。

結論

コイル内部の磁界は思ったより平坦だった。コイル線付近のみ値が大きくなっている。