Computational Science Asked on February 19, 2021
I want to compare the field of two electrical currents and compare the resulting field to a magnetic dipol field and find magnetic momentum that minimizes the difference of the two fields.
My current approach is to calculate the field of the two currents (by superimposing the fields of the two individual currents) and then start a optimization/search algorithm that minimizes the difference of the two fields. Because I am not interested in the immediate surroundings of the dipol or wires, I mask out the area around the center (and set the values of the fields there to zero).
To compute the overall error that the search algorithm should minimize, I subtract the dipol field from the field from the two currents. I use vectorization, so I end up just substracting whole fields from each other. Feels nice, but perhaps its wrong? I assume the vecorization substracts each vector in each coordinate. I want those individual vectors to be as similar as possible and want to minimize all those individual, small errors, added up. But I fear this optimization does not take place because the fields dont look similar.
How do i compare fields effectively?
here is my code so far, so you can see all the details of how i did it.
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize
# def magnetic_dipol_field_in_point(dipol_coordinates, dipol_vector, point_coordinates):
# magnetic_permeability = 4 * np.pi * 1e-7
# direction_vector = (point_coordinates - dipol_coordinates)
# direction_vector_norm = np.linalg.norm(direction_vector, axis=0)
# direction_unit_vector = np.divide(direction_vector, direction_vector_norm, out=np.zeros_like(direction_vector),
# where=direction_vector_norm != 0)
# vector_product = np.dot(dipol_vector, direction_unit_vector)
# unit_vect_prod_ = direction_unit_vector * vector_product
# dipol_vector = dipol_vector[:, np.newaxis]
# dipol_diff = unit_vect_prod_ - dipol_vector
# b_vector = (magnetic_permeability / (4 * np.pi)) * (3 * dipol_diff)
# / np.power(direction_vector_norm, 3)
# return b_vector
def magnetic_dipol_field_in_point(dipol_coordinates, dipol_vector, point_coordinates):
mu_0 = 4 * np.pi * 1e-7
direction_vector = (point_coordinates - dipol_coordinates)
direction_vector_norm = np.linalg.norm(direction_vector, axis=0)
direction_unit_vector = (direction_vector / direction_vector_norm)
b_vector = (mu_0 / (4 * np.pi)) * (3 * direction_unit_vector * np.dot(dipol_vector.T, direction_unit_vector)
- dipol_vector) / np.power(direction_vector_norm, 3)
return b_vector
def magnetic_vector_in_point(current_coordinates, current, point_coordinates):
magnetic_permeability = 4 * np.pi * 1e-7
direction_vector = (point_coordinates - current_coordinates)
direction_vector_norm = np.linalg.norm(direction_vector, axis=0)
direction_unit_vector = np.divide(direction_vector, direction_vector_norm, out=np.zeros_like(direction_vector),
where=direction_vector_norm != 0)
# orthogonal_unit_vector = np.flip(direction_unit_vector, axis=0) * np.array([[-1], [1]])
# faster version of the above:
# https://stackoverflow.com/questions/64795722/how-to-calculate-the-orthogonal-vector-of-a-unit-vector-with-numpy/64797015#64797015
orthogonal_unit_vector = direction_unit_vector[::-1, :]
np.negative(orthogonal_unit_vector[0, :], out=orthogonal_unit_vector[0, :])
b_vector = (current * magnetic_permeability / (2 * np.pi)) * np.divide(orthogonal_unit_vector,
direction_vector_norm,
out=np.zeros_like(orthogonal_unit_vector),
where=direction_vector_norm != 0)
return b_vector
def calc_current(time, num, peak):
return peak * np.cos(num * 2 * np.pi / 3 + time) # time is in the interval 0..2pi
def mask_area(field, grid, coord, radius):
x_center = coord[0]
y_center = coord[1]
field[0] = np.where(np.sqrt((grid[0] - x_center) ** 2 + (grid[1] - y_center) ** 2) > radius, field[0], 0)
field[1] = np.where(np.sqrt((grid[0] - x_center) ** 2 + (grid[1] - y_center) ** 2) > radius, field[1], 0)
return field
def main():
size = 4
a = 2.67e-2
X, Y = np.mgrid[-size:size:.1, -size:size:.1]
grid = np.vstack([X.ravel(), Y.ravel()])
current_coordinates = np.array([[a, -a / np.sqrt(3), -a / np.sqrt(3)],
[0, a, -a]])
magnetic_fields = np.zeros([2 + 1, grid.shape[0], grid.shape[1]])
currents = np.zeros(3)
for i in range(0, 2):
currents[i] = calc_current(1 + np.pi / 2, i, 1)
currents = [5, -5, 0]
magnetic_fields[i] = magnetic_vector_in_point(current_coordinates[:, np.newaxis, i], currents[i],
grid)
magnetic_fields[i] = mask_area(magnetic_fields[i], grid, [0.00564237, 0.01335], 1)
def find_dipol_moment(x, fields, slot_cnt, end_points, koordinaten_feld):
# calculate middle between currents, where dipol should be located.
dipol_coord = np.mean(end_points[:, 0:2], axis=1)
# calculate direction of the dipol moment, it is orthogonal to the line between the currents
vect = end_points[:, 0] - end_points[:, 1]
dipol_vektor = np.flip(vect, axis=0) * np.array([-1, 1])
#
dipol_vektor = dipol_vektor[:, np.newaxis]
dipol_unit_vektor = - dipol_vektor / np.linalg.norm(dipol_vektor)
dipol_coord = dipol_coord[:, np.newaxis]
fields[slot_cnt] = magnetic_dipol_field_in_point(dipol_coord, x[0] * dipol_unit_vektor, koordinaten_feld)
fields[slot_cnt] = -mask_area(fields[slot_cnt], grid, [0.00564237, 0.01335], 1)
error_field = np.sum(magnetic_fields, axis=0)
error = np.nansum(np.linalg.norm(error_field, axis=1))
print("Dipolmoment:", x,"tFehler: ", error)
return error
res = minimize(find_dipol_moment, -5, # 1 ist mein anfangswert bei der suche
args=(magnetic_fields, 2, current_coordinates, grid))
x = res.x
print("Ergebnis: äquivalentes Dipolmoment", x)
magnetic_field = np.sum(magnetic_fields[0:2], axis=0)
dipol_field = magnetic_fields[2]
# overall_field = np.sum(magnetic_fields, axis=0)
# magnetic_field = magnetic_fields[2]
fig, ax = plt.subplots(figsize=(10, 10), dpi=180, facecolor='w', edgecolor='k')
ax.scatter(0.00564237, 0.01335, marker="o")
ax.scatter(a, 0, marker="+")
ax.scatter(-a / np.sqrt(3), a, marker="+")
q = ax.quiver(grid[0, :], grid[1, :], magnetic_field[0, :], magnetic_field[1, :], color='b')
q = ax.quiver(grid[0, :], grid[1, :], dipol_field[0, :], dipol_field[1, :], color = 'k')
# q = ax.quiver(coordiante_field[0, :], coordiante_field[1, :], overall_field.T[0, :], overall_field.T[1, :],
# color='r')
# print(str(np.nansum(np.linalg.norm(overall_field, axis=1 )) ))
ax.quiverkey(q, X=0.3, Y=1.1, U=10,
label='Quiver key, length = 10', labelpos='E')
ax.set_aspect('equal', 'box')
plt.show()
if __name__ == '__main__':
main()
We consider two thin infinite straight parallel wires at distance $2 a$ apart and carrying equal currents $I$ in opposite directions. We need to find an approximate description of the system in terms of multipole expansion moments.
Let's first make a few preliminary comments. The magnetic field is described by the magnetic vector potential $vec{A}$ that is related to the magnetic field $vec{B}$ according to $nabla times vec{A} = vec{B}$, and $nabla^2 vec{A} = - mu_0 vec{j}$, where $vec{j}$ is the electric current (using SI units).
Let's use the coordinate $z$ in the direction of the wires. For an infinite straight wire in the $hat{e}_z$ direction the vector potential is also in the $hat{e}_z$ direction. Due to the symmetry, our problem is 2D and the magnetic field is fully described by a scalar quantity $A_z(x,y)$.
For a single wire located in the 2D plane at $vec{r}_1 =(x_1,y_1)$ and carrying current $I_1$ the vector potential is found from $nabla^2 A_z = -mu_0 I_1 delta(vec{r}-vec{r}_1)$, where on the right hand side is the delta-function, and the solution is well known to be $A_z = -(mu_0/2pi) I_1 ln(|vec{r}-vec{r}_1|)$, where the distance between the observation point and the wire is $|vec{r}-vec{r}_1|$.
Further it is convenient to use the polar coordinates $(rho,theta)$, thus a wire is described by the current $I_1$ and its location ($rho_1,theta_1$). Then at the observation point $(rho,theta)$, using the cosine theorem, we have
$ A_z = -frac{mu_0}{4pi} I_1 ln(rho^2 + rho_1^2 - 2 rho rho_1 cos(theta-theta_1)) $
For the far field, assuming $rho_1 ll rho$, we expand the expression as
$ A_z = -frac{mu_0}{4pi} I_1 left[ ln(rho^2) - ln(1 - frac{2 rho_1}{rho} cos(theta-theta_1)) right] approx -frac{mu_0}{4pi}I_1 ln(rho^2) -frac{mu_0}{4pi} I_1 frac{2 rho_1}{rho} cos(theta-theta_1) $
The second term on the right-hand side is the dipole term which can be written as $ -(mu_0/2pi) vec{p}_1 cdot vec{r} / r^2$, where the dipole moment is $vec{p}_1 = I_1 vec{r}_1$
Now, using superposition, we can immediately generalize for N parallel wires at locations described in the 2D plane by vectors $vec{r}_k$ and carrying currents $I_k$. The far field at location $vec{r}$, assuming $|r| gg |r_k|$, is given by
$ A_z = -frac{mu_0}{2pi} I_{tot} ln(|r|) -frac{mu_0}{2pi} (vec{p}_{tot} cdot vec{r})/r^2 $
where the total current is $I_{tot} = sum_{k=1}^N I_k$ and the total dipole moment is $vec{p}_{tot} = sum_{k=1}^N I_k vec{r}_k$.
Of course the expansion procedure could be continued further and the exact answer would be an infinite sum of multipole moments; but for the far-away field only the lowest-order terms matter.
Back to the situation with two wires at distance $2a$ apart, the total current is zero so the first (monopole) term vanishes. The total dipole moment is $p_{tot} = 2 a I$ and is in the direction from one wire to the other one, in the plane perpendicular to the wires. The result for $A_z$ is very similar to the electric potential from a regular electric dipole; only the dependence on the distance is different.
Finally, the magnetic field is found from $A_z$ as
$ vec{B} = nabla times vec{A} = hat{e}_z times nabla A_z $
Correct answer by Maxim Umansky on February 19, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP