Need to continue with B export

This commit is contained in:
schoener 2021-11-12 16:38:25 +01:00
parent 67f2588f60
commit 139d79c411
11 changed files with 365 additions and 52 deletions

View File

@ -65,6 +65,7 @@ AHH_Coil.cooling(I_grad, 22.5)
#AHH_Coil.B_quick_plot(I_grad)
#AHH_Coil.B_grad_quick_plot(I_grad)
AHH_Coil.plot_raster(raster_value= 11)
AHH_Coil.plot_3d(I,50,50)
# zero = mu_it(0)
# print(f"Bz_grad({z[zero]}) = {Bz_grad[zero]} G/cm")

View File

@ -0,0 +1,95 @@
"""
Created on 11.11.21
@author: Joschka
"""
import numpy as np
import matplotlib.pyplot as plt
from src import coil_class as BC
from src import physical_constants as cs
import logging as log
log.basicConfig(level=log.DEBUG, format='%(message)s')
def main():
AHH_Coil = BC.BCoil(HH=-1, distance=69.313, radius=47, layers=8, windings=16,
wire_height=0.5, wire_width=0.5, insulation_thickness=0.068/2,
is_round=True, winding_scheme=2)
x, z = BC.BCoil.make_axis(50, 5)
I_current = 5
# print(BC.BCoil._BCoil__B_z_loop(I_current = I, r_loop = 1, z_loop = 1, r_pos = 0, z_pos = 0.1)
# - BC.BCoil._BCoil__B_z_loop(I_current = I, r_loop = 1, z_loop = -1, r_pos = 0, z_pos = 0.1))
# print(BC.BCoil._BCoil__B_z_loop(I_current=I, r_loop=1, z_loop=1, r_pos=0, z_pos= -0.1)
# - BC.BCoil._BCoil__B_z_loop(I_current=I, r_loop=1, z_loop=-1, r_pos=0, z_pos=-0.1))
step = 1
xlim_mm = 30
xlim = int(xlim_mm/step) + 1
ylim_mm = 30
ylim = int(ylim_mm/step) + 1
zlim_mm = 20
zlim = int(zlim_mm/step) + 1
B_quarter = np.zeros((xlim, ylim, zlim, 2))
#print(B_quarter)
x = np.linspace(0,xlim_mm,xlim) * 1e-3
y = np.linspace(0,ylim_mm,ylim) * 1e-3
z = np.linspace(0,zlim_mm,zlim) * 1e-3
calc_raster = AHH_Coil.full_raster(1)
rastering_value = len(calc_raster[0])
I_current /= rastering_value
for wire in range(0,AHH_Coil.get_N()):
for ii in range(0,rastering_value):
# extract position information out of raster
z_pos = calc_raster[wire, ii, 0]
r_pos = calc_raster[wire, ii, 1]
for xx in range(0,xlim):
for yy in range(0,ylim):
r = np.sqrt(x[xx]**2 + y[yy]**2)
# calculate z-component at x,y
B_quarter[xx, yy, :, 1] += BC.BCoil._BCoil__B_z_loop(I_current, r_pos, z_pos, r, z) \
+ BC.BCoil._BCoil__B_z_loop(AHH_Coil.HH * I_current, r_pos, -z_pos, r, z)
if r == 0:
continue
# calculate r-component at x, y
B_quarter[xx, yy, :, 0] += BC.BCoil._BCoil__B_r_loop(I_current, r_pos, z_pos, r, z) \
+ BC.BCoil._BCoil__B_r_loop(AHH_Coil.HH * I_current, r_pos, -z_pos, r, z)
np.save('output/B_quarter.npy',B_quarter)
#B_tot = np.zeros(((2 * xlim) - 1, (2 * ylim) - 1, (2 * zlim)-1, 2))
#AHH_Coil.B_quick_plot(5,50)
#AHH_Coil.B_grad_quick_plot(5,50)
# Bz, Bx = AHH_Coil.B_field(I_current, x, z)
# B_tot_z, B_tot_x = AHH_Coil.B_field(I_current, x, z)
if __name__ == '__main__':
main()

View File

@ -0,0 +1,54 @@
import numpy as np
import matplotlib.pyplot as plt
from src import coil_class as BC
from src import physical_constants as cs
import logging as log
def main():
AHH_Coil = BC.BCoil(HH=-1, distance=69.313, radius=47, layers=8, windings=16,
wire_height=0.5, wire_width=0.5, insulation_thickness=0.068 / 2,
is_round=True, winding_scheme=2)
step = 1
xlim_mm = 30
xlim = int(xlim_mm/step) + 1
ylim_mm = 30
ylim = int(ylim_mm/step) + 1
zlim_mm = 20
zlim = int(zlim_mm/step) + 1
x = np.linspace(0,xlim_mm,xlim) * 1e-3
y = np.linspace(0,ylim_mm,ylim) * 1e-3
z = np.linspace(0,zlim_mm,zlim) * 1e-3
B_quarter = np.load('output/B_quarter.npy')
x_m, z_m = np.meshgrid(x, z)
plt.figure(34)
plt.quiver(x_m,z_m,np.transpose(B_quarter[0, :, :, 0]), np.transpose(B_quarter[0, :, :, 1]))
plt.xlabel("x-axis [mm]")
plt.ylabel("z-axis [mm]")
plt.show()
x_m, y_m = np.meshgrid(x, y)
plt.figure(35)
plt.quiver(x_m, y_m, np.transpose(B_quarter[:, :, 0, 0]), np.transpose(B_quarter[:, :, 0, 1]))
plt.xlabel("x-axis [mm]")
plt.ylabel("y-axis [mm]")
plt.show()
# AHH_Coil.B_quick_plot(5,abs = False)
#AHH_Coil.plot_3d(5,30,20)
if __name__ == '__main__':
main()

View File

@ -0,0 +1,79 @@
"""
Created on 11.11.21
@author: Joschka
"""
import numpy as np
import matplotlib.pyplot as plt
from src import coil_class as BC
from src import physical_constants as cs
import logging as log
log.basicConfig(level=log.DEBUG, format='%(message)s')
def main():
step = 1
xlim_mm = 30
xlim = int(xlim_mm / step) + 1
ylim_mm = 30
ylim = int(ylim_mm / step) + 1
zlim_mm = 20
zlim = int(zlim_mm / step) + 1
x = np.linspace(-xlim_mm, xlim_mm, 2*xlim + 1) * 1e-3
y = np.linspace(-ylim_mm, ylim_mm, 2*ylim+1) * 1e-3
z = np.linspace(-zlim_mm, zlim_mm, 2*zlim+1) * 1e-3
B_quarter = np.load('output/B_quarter.npy')
step = 1
xlim_mm = 30
xlim = 2 * int(xlim_mm/step) + 1
ylim_mm = 30
ylim = 2 * int(ylim_mm/step) + 1
zlim_mm = 20
zlim = 2 * int(zlim_mm/step) + 1
#B_quarter = np.zeros((xlim, ylim, zlim, 2))
#print(B_quarter)
x = np.linspace(-xlim_mm,xlim_mm,xlim) * 1e-3
y = np.linspace(-ylim_mm,ylim_mm,ylim) * 1e-3
z = np.linspace(-zlim_mm,zlim_mm,zlim) * 1e-3
B_tot = np.zeros((xlim , ylim , zlim ,2))
# First quater half x>0,y>0,z>0
B_tot[xlim//2: xlim, ylim//2: ylim, zlim//2 : zlim,1] = B_quarter[:,:,:,1]
print(B_tot)
for xx in range(xlim//2, xlim):
for yy in range(ylim//2,ylim):
pass
#AHH_Coil.B_quick_plot(5,50)
#AHH_Coil.B_grad_quick_plot(5,50)
# Bz, Bx = AHH_Coil.B_field(I_current, x, z)
# B_tot_z, B_tot_x = AHH_Coil.B_field(I_current, x, z)
if __name__ == '__main__':
main()

Binary file not shown.

View File

@ -66,6 +66,7 @@ print(f"Diff B 15 mm: {Bz[30000] - Bz[15000]}, relative: {(Bz[30000] - Bz[15000]
HH_Coil.cooling(I_current,25)
HH_Coil.print_info()
#HH_Coil.B_curv_quick_plot(I_current)
#HH_Coil.plot_raster()

View File

@ -22,21 +22,21 @@ z = np.linspace(-15, 15, 30001)
# New coil
# Wire_1 = [0.5, 0.568]
Wire_1 = [0.45, 0.514]
Wire_1 = [0.5, 0.568]
#Wire_1 = [0.45, 0.514]
#I_current = 0.94
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 9, wire_height = Wire_1[0], wire_width = Wire_1[0], insulation_thickness=(Wire_1[1] - Wire_1[0]) / 2, is_round = True, winding_scheme= 2)
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8, wire_height = Wire_1[0], wire_width = Wire_1[0], insulation_thickness=(Wire_1[1] - Wire_1[0]) / 2, is_round = True, winding_scheme= 2)
I_current = 64 / HH_Coil.get_N() * 1.25
#set radius plus distance
HH_Coil.set_R_outer(50.5 - HH_Coil.get_tot_wire_width()*1e3)
HH_Coil.set_d_min(47.15)
HH_Coil.set_R_outer(50.5 - HH_Coil.get_tot_wire_width()*1e3 - 0.5)
HH_Coil.set_d_min(47.15+1)
print(f"height = {HH_Coil.get_coil_height()}")
HH_Coil.print_info()
Bz, Bx = HH_Coil.B_field(I_current, x, z, raster = 10)
Bz, Bx = HH_Coil.B_field(I_current, x, z, raster = 7)
B_tot_z, B_tot_x = HH_Coil.B_field(I_current, x, z, raster = 10)
B_tot_z, B_tot_x = HH_Coil.B_field(I_current, x, z, raster = 7)
Bz_curv = BC.BCoil.curv(Bz, z)
HH_Coil.cooling(I_current,28)

View File

@ -19,6 +19,7 @@ z = np.linspace(-lim,lim,nr_points)
def mu_it(x_pos):
step =
it = nr_points//2 + x_pos
return it

View File

@ -4,29 +4,40 @@ Created on 3.11.21
@author: Joschka
"""
import numpy as np
import matplotlib.pyplot as plt
from src import coil_class as BC
from src import physical_constants as cs
def Q_heat(flow,d_T):
V_t = 4.8 * 3.2 * 1e-6 * flow
m = cs.water_dens * V_t
Q = m * cs.water_c_p * d_T
return Q
def main():
dens = 999.78 # kg/m^3
c_p = 4190 # J/kg*K
d_T = 3
flow = 0.3#/5 #m/s
d_T = 1.5
flow = 0.3/5#/5 #m/s
#flow *=2
print(f"flow = {flow}m/s")
V_t = 4.8 * 3.2 * 1e-6 * flow
print(f"Volume rate = {V_t * 1e6} mL/s")
m = dens * V_t
Q = m * c_p * d_T
m = cs.water_dens * V_t
Q = m * cs.water_c_p * d_T
print(f"Q = {Q} J/s")
HH_Coil = BC.BCoil(HH=1, distance=51.694, radius=47.9263, layers=8, windings=8, wire_height=0.5,
#flow = np.linspace(0,5,100)
#plt.plot(flow,Q_heat(flow,3))
#plt.show()
HH_Coil = BC.BCoil(HH=1, distance=51.694, radius=47.9263, layers=8, windings=16, wire_height=0.5,
wire_width=0.5, insulation_thickness=0.034, is_round=True,
winding_scheme=2)
for I in np.arange(1,7,0.001):
for I in np.arange(1,20,0.001):
P = HH_Coil.power(I, 23)
if P > Q:
break
@ -37,6 +48,8 @@ def main():
print(f"max field = {B_max} G")
print(HH_Coil.cooling(I, 25))

View File

@ -322,30 +322,89 @@ class BCoil:
B *= 1e4 # conversion Gauss
return B
@staticmethod
def make_axis(lim,precision,two_axis=True):
"""
Gives back two arrays x and z in mm, symmetric around zero
Args:
lim: positive (and negative if lim2 = None) limit of axis
precision: nr_points per mm
two_axis: If True --> gives back two arrays with given parameter
Returns: array
"""
nr_points = (2 * lim) * precision + 1
z = np.linspace(-lim, lim, nr_points)
if two_axis:
x = np.linspace(-lim, lim, nr_points)
return x, z
@staticmethod
def mm_it(x, x_pos):
"""
Takes argument in mm returns position in x-array
Args:
x: array of interest
x_pos: position to find position of
Returns: int of position of x_pos in array x
"""
precision = len(x)//(2*np.abs(x[0]))
it = int(len(x)//2 + precision * x_pos)
return it
@staticmethod
def k_sq(R_loop, z_loc, r, z):
""" Argument for elliptic integrals"""
return (4 * R_loop * r) / ((R_loop + r) ** 2 + (z - z_loc) ** 2)
@staticmethod
def B_z_loop(I_current, R_loop, z_loc, r, z):
"""calculate z-component of B-field at position r and z for each individual loop"""
B_z = 2e-7 * I_current * 1 / ((R_loop + r) ** 2 + (z - z_loc) ** 2) ** (1 / 2) * (
sp.ellipk(BCoil.k_sq(R_loop, z_loc, r, z)) + sp.ellipe(BCoil.k_sq(R_loop, z_loc, r, z)) * (
R_loop ** 2 - r ** 2 - (z - z_loc) ** 2) / ((R_loop - r) ** 2 + (z - z_loc) ** 2))
def __B_z_loop(I_current, r_loop, z_loop, r_pos, z_pos):
"""
calculate z-component of B-field at position r and z for each individual loop
Args:
I_current: Current in [A]
r_loop: Radial position of current loop
z_loop: Axial position of current loop along symmetry axis
r_pos: radial position of calculation in [m]
z_pos: axial position of calculation in [m]
Returns: z component of B-field at position r, z in G
"""
B_z = 2e-7 * I_current * 1 / ((r_loop + r_pos) ** 2 + (z_pos - z_loop) ** 2) ** (1 / 2) * (
sp.ellipk(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) + sp.ellipe(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) * (
r_loop ** 2 - r_pos ** 2 - (z_pos - z_loop) ** 2) / ((r_loop - r_pos) ** 2 + (z_pos - z_loop) ** 2))
B_z *= 1e4 # conversion to gauss
return B_z
@staticmethod
def B_r_loop(I_current, R_loop, z_loc, r, z):
"""calculate r-component of B-field at position r and z for each individual loop"""
B_r = 2e-7 * I_current / r * (z - z_loc) / ((R_loop + r) ** 2 + (z - z_loc) ** 2) ** (1 / 2) * (
-sp.ellipk(BCoil.k_sq(R_loop, z_loc, r, z)) + sp.ellipe(BCoil.k_sq(R_loop, z_loc, r, z)) * (
R_loop ** 2 + r ** 2 + (z - z_loc) ** 2) / ((R_loop - r) ** 2 + (z - z_loc) ** 2))
def __B_r_loop(I_current, r_loop, z_loop, r_pos, z_pos):
"""
calculate z-component of B-field at position r and z for each individual loop
Args:
I_current: Current in [A]
r_loop: Radial position of current loop
z_loop: Axial position of current loop along symmetry axis
r_pos: radial position of calculation in [m]
z_pos: axial position of calculation in [m]
Returns: z component of B-field at position r, z in G
"""
B_r = 2e-7 * I_current / r_pos * (z_pos - z_loop) / ((r_loop + r_pos) ** 2 + (z_pos - z_loop) ** 2) ** (1 / 2) * (
-sp.ellipk(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) + sp.ellipe(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) * (
r_loop ** 2 + r_pos ** 2 + (z_pos - z_loop) ** 2) / ((r_loop - r_pos) ** 2 + (z_pos - z_loop) ** 2))
B_r *= 1e4 # conversion to gauss
return B_r
def B_field(self, I_current, x, z, raster=10):
def B_field_z(self,I_current):
pass
def B_field(self, I_current, x, z, raster = 7):
"""
Returns Bz along z-axis and B_x along x-axis,
HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
@ -380,11 +439,11 @@ class BCoil:
z_pos = calc_raster[wire, ii, 0]
r_pos = calc_raster[wire, ii, 1]
B_z += BCoil.B_z_loop(I_current, r_pos, z_pos, 0, z_SI) + BCoil.B_z_loop(self.HH * I_current,
r_pos, -z_pos, 0, z_SI)
B_x_pos += BCoil.B_r_loop(I_current, r_pos, z_pos, x_positive, 0) + BCoil.B_r_loop(
B_z += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, z_SI) + BCoil.__B_z_loop(self.HH * I_current, r_pos,
-z_pos, 0, z_SI)
B_x_pos += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_positive, 0) + BCoil.__B_r_loop(
self.HH * I_current, r_pos, -z_pos, x_positive, 0)
B_x_neg += BCoil.B_r_loop(I_current, r_pos, z_pos, x_negative, 0) + BCoil.B_r_loop(
B_x_neg += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_negative, 0) + BCoil.__B_r_loop(
self.HH * I_current, r_pos, -z_pos, x_negative, 0)
B_x_neg *= -1 # B_x_neg is pointing in opposite x_direction
@ -393,7 +452,7 @@ class BCoil:
# print(f"time = {end - start} s")
return B_z, B_x
def max_field(self, I_current,raster = 7):
def max_field(self, I_current, raster = 7):
B_z_0 = 0
calc_raster = self.full_raster(raster)
@ -410,10 +469,11 @@ class BCoil:
z_pos = calc_raster[wire, ii, 0]
r_pos = calc_raster[wire, ii, 1]
B_z_0 += BCoil.B_z_loop(I_current, r_pos, z_pos, 0, 0) + BCoil.B_z_loop(self.HH * I_current,
r_pos, -z_pos, 0, 0)
B_z_0 += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, 0) + BCoil.__B_z_loop(self.HH * I_current, r_pos,
-z_pos, 0, 0)
return B_z_0
def B_tot_along_axis(self, I_current, x, z, raster=10):
"""
return B_tot_z, B_tot_x
@ -453,17 +513,19 @@ class BCoil:
r_pos = calc_raster[wire, ii, 1]
# z-field along z-axis (x-Field always zero)
B_tot_z += BCoil.B_z_loop(I_current, r_pos, z_pos, 0, z_SI) + BCoil.B_z_loop(
self.HH * I_current, r_pos, -z_pos, 0, z_SI)
B_tot_z += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, z_SI) + BCoil.__B_z_loop(self.HH * I_current,
r_pos, -z_pos, 0, z_SI)
# x-field along x-axis
B_x_pos += BCoil.B_r_loop(I_current, r_pos, z_pos, x_pos, 0) + BCoil.B_r_loop(
self.HH * I_current, r_pos, -z_pos, x_pos, 0)
B_x_neg += BCoil.B_r_loop(I_current, r_pos, z_pos, x_neg, 0) + BCoil.B_r_loop(
self.HH * I_current, r_pos, -z_pos, x_neg, 0)
B_x_pos += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_pos, 0) + BCoil.__B_r_loop(self.HH * I_current,
r_pos, -z_pos, x_pos,
0)
B_x_neg += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_neg, 0) + BCoil.__B_r_loop(self.HH * I_current,
r_pos, -z_pos, x_neg,
0)
# Bz along x-axis:
B_z_x += BCoil.B_z_loop(I_current, r_pos, z_pos, x_SI, 0) + BCoil.B_z_loop(self.HH * I_current,
r_pos, -z_pos, x_SI, 0)
B_z_x += BCoil.__B_z_loop(I_current, r_pos, z_pos, x_SI, 0) + BCoil.__B_z_loop(self.HH * I_current,
r_pos, -z_pos, x_SI, 0)
B_x_neg *= -1 # B_x_neg is pointing in opposite x_direction
@ -499,18 +561,16 @@ class BCoil:
r_pos = calc_raster[wire, ii, 1]
# compute z-value of field
B[:, el, 1] += BCoil.B_z_loop(I_current, r_pos, z_pos, np.abs(x_SI[el]),
z_SI) + BCoil.B_z_loop(self.HH * I_current, r_pos, -z_pos,
np.abs(x_SI[el]), z_SI)
B[:, el, 1] += BCoil.__B_z_loop(I_current, r_pos, z_pos, np.abs(x_SI[el]), z_SI) + BCoil.__B_z_loop(
self.HH * I_current, r_pos, -z_pos, np.abs(x_SI[el]), z_SI)
# compute x-value
if x_SI[el] < 0:
B[:, el, 0] -= BCoil.B_r_loop(I_current, r_pos, z_pos, -x_SI[el],
z_SI) + BCoil.B_r_loop(self.HH * I_current, r_pos, -z_pos,
-x_SI[el], z_SI)
B[:, el, 0] -= BCoil.__B_r_loop(I_current, r_pos, z_pos, -x_SI[el], z_SI) + BCoil.__B_r_loop(
self.HH * I_current, r_pos, -z_pos, -x_SI[el], z_SI)
elif x_SI[el] > 0:
B[:, el, 0] += BCoil.B_r_loop(I_current, r_pos, z_pos, x_SI[el], z_SI) + BCoil.B_r_loop(
B[:, el, 0] += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_SI[el], z_SI) + BCoil.__B_r_loop(
self.HH * I_current, r_pos, -z_pos, x_SI[el], z_SI)
elif x_SI[el] == 0:
@ -524,18 +584,19 @@ class BCoil:
def plot_3d(self, I_current, x_lim, z_lim):
x = np.arange(-x_lim, x_lim, 5)
print(x)
z = np.arange(-z_lim, z_lim, 5)
print(z)
x_m, z_m = np.meshgrid(x, z)
if self.is_round:
B = self.B_multiple_3d(I_current, x, z, raster=3)
else:
B = self.B_multiple_3d(I_current, x, z, raster=2)
B_tot = BCoil.B_tot_3d(B)
for xx in range(0, len(x)):
for zz in range(0, len(z)):
if B_tot[xx, zz] > 15:
B[xx, zz, :] /= B_tot[xx, zz] / 15
if B_tot[zz, xx] > 15:
B[zz, xx, :] /= B_tot[zz, xx] / 15
plt.figure(33)
plt.quiver(x_m, z_m, B[:, :, 0], B[:, :, 1])
@ -659,8 +720,13 @@ class BCoil:
Power = self.power(I_current, T)
Voltage = self.resistance(T) * I_current
print("")
print("Cooling:")
print(f" current density = {j_dens} A/mm^2")
print(f" Power = {Power} W")
print(f" U = {Voltage} V")
def power(self, I_current, T):
P = self.resistance(T) * I_current ** 2

View File

@ -29,3 +29,6 @@ lambda_626 = 626.082e-9
rho_copper_20 = 1.68e-8
density_copper = 8.940e3
# water
water_dens = 999.78 # kg/m^3
water_c_p = 4190 # J/kg*K