Coil simulator (mono coil)

モノコイルのシミュレーター

### Pulse magnet simulator v0.1 
### by A. Ikeda (UEC-Tokyo, 30 April 2022)

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



# # # # # # # # # # # # # # # # # # # # # # # # 
#                                             #
#              Imports                        #
#                                             #
# # # # # # # # # # # # # # # # # # # # # # # # 

import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline 
# plt.style.use('./ai.mplstyle')
# %matplotlib tk 

res_x = []
res_y = []

reserve_v = []
reserve_i = []
reserve_tp = []
reserve_t = []
reserve_bf = []

# # # # # # # # # # # # # # # # # # # # # # # # 
#                                             #
#              Input parameter                #
#                                             #
# # # # # # # # # # # # # # # # # # # # # # # # 


wire = 1. # wire diameter (mm)
wire_cover = 0.033 # wire insulation (mm)
n = 16 # windings
m = 8 # layer
bore = 8 #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_single = 800 * 1e-6 #(F)
cap_num = 3
cap = cap_single * cap_num #(F)
v0 = 2000 #(V)
q0 = v0 * cap #(C)
tp0 = 78 #(Kelvin)

t0 = 0.0 # when simulation starts
t1 = 60e-3 # when simulation ends
del_t = 5e-6 # time div


F_t = 0 #thickness of the flange
z = 0 # position of magnetic field measurement (mm)
R_r = 1 #residual resistance (mΩ)
L_r = 2.5 #residual L (μH)
# L_r = 100 #residual L (μH)



# # # # # # # # # # # # # # # # # # # # # # # # 
#                                             #
#           Processing parameter              #
#                                             #
# # # # # # # # # # # # # # # # # # # # # # # # 

# 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 
# 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_r * 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)

crosssection_cm = (wire / 2) ** 2 * np.pi / 100
length_cm = l_length * 100
weight = crosssection_cm * length_cm * 8.9

# Position of B-sensor from the field center in the z axis (m)
z0 = z * 1e-3 
#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 



# # # # # # # # # # # # # # # # # # # # # # # # 
#                                             #
#              Calculations                   #
#                                             #
# # # # # # # # # # # # # # # # # # # # # # # # 


# Definition of functions

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.trapezoid(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', [tpoints, ipoints, bf, qpoints / cap, tppoints])



maxx = np.argmax(bf)
maxt = round(tpoints[maxx] * 1000, 2)
plsx = np.argmin(np.abs(bf[maxx * 2: maxx * 3])) + maxx * 2
plswid = round(tpoints[plsx] * 1000, 2)


imax = round(ipoints[maxx] / 1000, 1)
vmax = v0
bmax = round(bf[maxx], 1)

imin = round(ipoints[int(plsx * 1.1)] / 1000, 1)
vmin = round(qpoints[int(plsx * 1.1)] / cap, 1)
bmin = round(bf[int(plsx * 1.1)], 1)

tmax = round(tppoints[plsx], 1)
tmin = tp0

efficiency = round(bmax / imax, 1)
tpos1 = maxt * 0.5
tpos2 = maxt * 1.4

llabel = [['$B$ (T)', '$V$ (V)'], ['$I$ (kA)', '$T$ (K)']]
maxs = [[bmax, vmax], [imax, tmax]]
mins = [[bmin, vmin], [imin, tmin]]

fig, axs = plt.subplots(2, 2, figsize = (7, 5), tight_layout = True)

if len(reserve_t) > 1:
    axs[0, 1].plot (reserve_t * 1000, reserve_v, lw = 0.5, color = 'red')
    axs[1, 0].plot (reserve_t * 1000, reserve_i / 1000, lw = 0.5, color = 'blue')
    axs[0, 0].plot (reserve_t * 1000, reserve_bf, lw = 0.5, color = 'green')
    axs[1, 1].plot (reserve_t * 1000, reserve_tp, lw = 0.5, color = 'black')



axs[0, 1].plot (tpoints * 1000, qpoints / cap, '.', color = 'red')
axs[0, 1].text(tpos2, vmax * 0.8, f'max: {maxt} ms')
axs[0, 1].text(tpos2, vmax * 0.5, f'pulse: {plswid} ms')
axs[1, 0].plot (tpoints * 1000, ipoints / 1000, '.', color = 'blue')
axs[1, 0].text(tpos1, imax * 0.1, f'I_max: {imax} kA')
axs[0, 0].plot (tpoints * 1000, bf, '.', color = 'green')
axs[0, 0].text(tpos1, bmax * 0.4, f'B_max: {bmax} T')
axs[0, 0].text(tpos1, bmax * 0.1, f'Efficiency: {efficiency} T/kA')
axs[1, 1].plot (tpoints * 1000, tppoints, '.', color = 'black')
axs[1, 1].text(tpos2, tmin * 1.1, f'T_max: {tmax} K')

for i in range(2):
    for j in range(2):
        axs[i, j].set_xlabel('$t$ (ms)', fontsize = 12)
        axs[i, j].set_ylabel(llabel[i][j], fontsize = 12)
        axs[i, j].set_ylim(mins[i][j] - abs(mins[i][j]) * 0.2, maxs[i][j] * 1.2)
        axs[i, j].set_xlim(-0.1, plswid * 1.1)

fig.suptitle('Pulse magnet simulator')
plt.tight_layout()
# plt.savefig('fig.pdf')
plt.show()

#### print ####

print('---Inputs---')
print(f'coil_height is {round(coil_length, 1)} mm.')
print(f'coil_diameter is {round(c_r, 1)} mm.')
print(f'bore is {round(bore, 1)} mm.')

print(f'L is {round(L * 1e6, 1)} μH.')
print(f'(L_res is {round(L_r, 1)} μH.)')
print(f'C is {round(cap * 1e6, 1)} μF.')
print(f'V_0 is {round(v0, 1)} V.')
# print(f'L_res is {round(L_res, 1)} μH.')
print(f'R_res is {round(R_r, 1)} mΩ.')

print('\n---Results---')
print(f'Total length of magnet wire is {round(l_length, 1)} m.')
print(f'Weight of magnet is {round(weight, 1)} g.')
print(f'I_max is {round(np.amax(ipoints))} A.')
print(f'B_max is {round(np.amax(bf), 1)} T.')
print(f'T_ini is {round(tp0, 1)} K.')
print(f'T_end is {round(tp, 1)} K.')
print(f'R_ini is {round(r(tp0) * 1000, 1)} mohm.')
print(f'R_end is {round(r(tp) * 1000, 1)} mohm.')
print(f'C_v  is {round(hcap(tp), 1)} J/K')
print(f'E_tot  is {round(cap * v0 ** 2 / 2)} J')


#### print ####




reserve_v = qpoints / cap
reserve_i = ipoints
reserve_tp = tppoints
reserve_t = tpoints
reserve_bf = bf