パルス磁場シミュレーターver0.1をpythonで

実際にミニコイルに電流を流して磁場を発生するとなると、どれくらい巻くとどれくらいの磁場が出るのか、どれくらい電流が流れるのか、どれくらい発熱するのか、が気になってきます。やる前からある程度わかる方が、時間とお金の節約になります。さしあたってはサイリスタのスペックを決めるために概算が必要。

というわけでミニコイル用のシミュレータプログラムをpythonでつくってみました。というか、M田研のigor用プログラムを少し改良しつつ移植しました。改良部分は、微分方程式の解析解を使わずに、ルンゲクッタ法で数値的に解くようにしました。とりあえず動くようになりました。igorのと磁場が結構ずれるので、まだなんか間違っているかも。評価が必要です。問題点が解消したら、更新したコードをのせます。抵抗と発熱を入れると少しずれはじめるので、デバイモデル付近が怪しいのかも。

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


"""
Pulse magnet simulator v0.1
by A. Ikeda of UEC TOKYO
30 April 2022


references
[1] Example program (3) in https://qiita.com/sci_Haru/items/b8d5c9cfe64c4366630a
[2] A program in Master thesis, Takeshi Nakamura of Y. H. Matsuda lab., Univ. of Tokyo (2010)
[3] Pulse magnet for high field generation, N. Miura (2008)
"""



###input parameters ###

wire = 1.0 # wire diameter
wire_cover = 0.033 # wire insulation
n = 10 # windings
m = 10 # layer
bore = 10 #inner diameter (mm)

# resitivity is resty = resty1 * tp + resty0
resty0 = 5 * 1e-10 # (ohm m)
resty1 = 6.5 * 1e-11 # (ohm m / T)
T_D = 343 # (K) Debye temperature of Cu 

cap = 1000e-6 # (F)
v0 = 1000 #(V)
q0 = v0 * cap #(C)
tp0 = 300 #(Kelvin)

t0 = 0.0 # when simulation ends
t1 = 5e-3 # when simulation starts
N = 200
del_t = (t1 - t0) / N # time grid


F_t = 0 #thickness of the flange
z = 0 # z position
R_r = 0 #residual resistance

###imput parameters ###



### processing parameters ###

# inner radii (m)
r1 = (bore * 1e-3) / 2 
# outer radii (m) 
r2 = m * (wire + wire_cover * 2) * 1e-3 + r1 
# outer diameter (mm) 
c_r = r2 * 2 * 1e3 
# half height of coil
l1 = (wire + wire_cover * 2) * 1e-3 * n / 2 
#L calib factor
L_t = 0 
# height of coil
coil_length = l1 * 2 * 1e3 

alpha = r2 / r1 # alpha
beta = l1 / r1 # beta
mu0 = 4 * np.pi * 1e-7 # permitivity

# Nagaoka's constant
K = (1 + 0.225 * (alpha + 1) / beta + 0.64 * (alpha - 1) / (alpha + 1) + 0.42 * (alpha - 1) / beta) ** -1
# L Inductance of the coil (H)
L = mu0 * np.pi / 8 * r1 * (alpha + 1) ** 2 / beta * (n * m) ** 2 * K + L_t * 1e-6
# Crosssection of wire
S = (wire * 1e-3 / 2) ** 2 * np.pi 
#length of wire
l_length = 0 
for m0 in np.arange(m): 
    l_length += 2 * np.pi * n * (r1 + (m0 + 1 / 2) * (wire + wire_cover * 2) * 1e-3)
# position on z axis from center (m)
z0 = z * 1e-3 

### processing parameters ###




#total mol of wire
R_const = 8.31447 #J/K mol, 　Gas Constant
density = 8920 * 1e3 # Density of Copper g/m3
mol = density * l_length * S / 63.536 


def resty(tp):
    return resty0 + tp * resty1 

def r(tp):
    return resty(tp) * l_length / S

def f(q, i, tp, t):
    return -q / cap -r(tp) * i

def debye(x):
    return x ** 4 * np.exp(x) / (np.exp(x) - 1) ** 2

qpoints = np.array([])
ipoints = np.array([])
tppoints = np.array([])
tpoints = np.arange(t0, t1, del_t)


i0 = 0.0 #(A)
q, i, tp = q0, i0, tp0

# heat capacity calculation using Debye model

def hcap(tx):
    x0 = np.arange(1e-3, T_D / tx, 1e-3)
    return 9 * R_const * (tx / T_D) ** 3 * np.trapz(debye(x0), x0) * mol

for t in tpoints:
    qpoints = np.append(qpoints, q)
    ipoints = np.append(ipoints, i)
    tppoints = np.append(tppoints, tp)
    
    k1i = f(q, i, tp, t) * del_t / L
    k1q = i * del_t

    k2i = f(q + k1q / 2, i + k1i / 2, tp, t + del_t / 2) * del_t / L
    k2q = (i + k1i / 2) * del_t 

    k3i = f(q + k2q / 2, i + k2i / 2, tp, t + del_t / 2) * del_t / L
    k3q = (i + k2i / 2) * del_t 

    k4i = f(q + k3q, i + k3i, tp, t + del_t ) * del_t / L
    k4q = (i + k3i) * del_t 


    tp += r(tp) * i ** 2 * del_t / hcap(tp)
    i += (k1i + 2 * k2i + 2 * k3i + k4i) / 6
    q += (k1q + 2 * k2q + 2 * k3q + k4q) / 6

    
ipoints = -1 * ipoints
delta = z0 / r1
KBp = (beta + delta) / (alpha - 1) * np.log((alpha + (alpha ** 2 + (beta + delta) ** 2) ** 0.5) / (1 + (1 + (beta + delta) ** 2) ** 0.5))
KBm = (beta - delta) / (alpha - 1) * np.log((alpha + (alpha ** 2 + (beta - delta) ** 2) ** 0.5) / (1 + (1 + (beta - delta) ** 2) ** 0.5))
bf = mu0 * n * m * ipoints /(4 * l1) * (KBp + KBm)


np.savetxt('out.csv', bf)

fig, axs = plt.subplots(2, 2, figsize = (9, 6), constrained_layout=True)
axs[0, 1].plot (tpoints * 1000, qpoints / cap, '.', color = 'red')
axs[0, 1].set_xlabel("$t$ (ms)", fontsize = 12)
axs[0, 1].set_ylabel("$V$ (V)", fontsize = 12)
axs[1, 0].plot (tpoints * 1000, ipoints / 1000, '.', color = 'blue')
axs[1, 0].set_xlabel("$t$ (ms)", fontsize = 12)
axs[1, 0].set_ylabel("$I$ (kA)", fontsize = 12)
axs[0, 0].plot (tpoints * 1000, bf, '.', color = 'green')
axs[0, 0].set_xlabel("$t$ (ms)", fontsize = 12)
axs[0, 0].set_ylabel("$B$ (T)", fontsize = 12)
axs[1, 1].plot (tpoints * 1000, tppoints, '.', color = 'black')
axs[1, 1].set_xlabel("$t$ (ms)", fontsize = 12)
axs[1, 1].set_ylabel('$T$ (K)', fontsize = 12)

plt.savefig('fig.pdf')
plt.show()

#### print ####

print("coil_length is " + str(round(coil_length, 1)) + " mm.")
print("coil_diameter is " + str(round(c_r, 1)) + " mm.")
print("bore is " + str(round(bore, 1)) + " mm.")
print("wire length is " + str(round(l_length, 1)) + " m.")


print("L is " + str(round(L * 1e6, 1)) + " μH.")
print("C is " + str(round(cap * 1e6, 1)) + " μF.")
print("V0 is " + str(round(v0, 1)) + " V.")

print("Imax is " + str(round(np.amax(ipoints))) + " A.")
print("Bmax is " + str(round(np.amax(bf), 1)) + " T.")
print("T_ini is " + str(round(tp0, 1)) + " K.")
print("T_end is " + str(round(tp, 1)) + " K.")
print("R_ini is " + str(round(r(tp0) * 1000, 1)) + " mohm.")
print("R_end is " + str(round(r(tp) * 1000, 1)) + " mohm.")
print("C_v  is " + str(round(hcap(tp), 2)) + " J/T")

#### print ####