1000テスラ理論の会
時之栖という宿泊施設でしたが、とっても素晴らしい場所でした。少年や社会人のサッカーの合宿で来ている人がたくさんいたように思います。家族でも行ってみたいと思いました。
実際にミニコイルに電流を流して磁場を発生するとなると、どれくらい巻くとどれくらいの磁場が出るのか、どれくらい電流が流れるのか、どれくらい発熱するのか、が気になってきます。やる前からある程度わかる方が、時間とお金の節約になります。さしあたってはサイリスタのスペックを決めるために概算が必要。
というわけでミニコイル用のシミュレータプログラムを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 ####
コメント
コメントを投稿