# Numerical simulation of magnetic dipole in a inhomogeneous magnetic field

Goal:

I want to use Python to illustrate how a magnetic dipole with magnetic moment m2 moves in a non-homogeneous magnetic field in a 2D-Plane. This field is generated by another magnetic dipole with moment m1. However, I get results that are not plausible.

For the calculation of the trajectory of m2 I use the equations from wiki - force between magnets

Theoretical Background:

Given is a magnetic dipole with constant moment m1 and constant position r1!. m1 generates a magnetic flux density

$${\displaystyle \mathbf {B}_1 (\mathbf {r} )={\frac {\mu _{0}}{4\pi }}{\frac {3\mathbf {\hat {r}} (\mathbf {\hat {r}} \cdot \mathbf {m} )-\mathbf {m} }{|\mathbf {r} |^{3}}}\qquad (1)}$$

where r is a vector from the center of the magnetic dipole to the location where the magnetic field is measured. In this non-homogeneous magnetic field a magnetic dipole with moment m2 and position r2 is placed. Its magnetic field can be neglected B2 = 0.

Only the dynamics of m2 shall be considered. According to the linked Wikipedia page, m2 will be subjected to a torque D and a force F. The torque is neglected D=0. The magnetic moment m2 is therefore constant.

Due to the acting force, the position vector r_2(t) depends on time. Its time evolution is obtained from the solution of the equation of motions, which is my objective:

$$\frac{d^2}{dt^2}\mathbf{r}_2 (t) = \binom{\dot{r}_{2x}}{\dot{r}_{2y}} = \mathbf{F}$$

According to Wikipedia, m1 exerts the following force on m2: (r = r2 - r1, r = ||r||_2):

\begin{align} \mathbf{F} &= \nabla |\mathbf{m}_2\cdot \mathbf{B}_1| \qquad(3) \\\\ &= {\frac {1}{ r^{5}}}\left[(\mathbf {m} _{1}\cdot \mathbf {r} )\mathbf {m} _{2}+(\mathbf {m} _{2}\cdot \mathbf {r} )\mathbf {m} _{1}+(\mathbf {m} _{1}\cdot \mathbf {m} _{2})\mathbf {r} -{\frac {5(\mathbf {m} _{1}\cdot \mathbf {r} )(\mathbf {m} _{2}\cdot \mathbf {r} )}{r^{2}}}\mathbf {r} \right] \qquad (4) \\ &= \frac{\mathbf{m}_1 \mathbf{r}}{r^5} \mathbf{m}_2 + \frac{\mathbf{m}_2\mathbf{r}}{r^5}+\frac{\mathbf{m}_1\mathbf{m}_2}{r^5}\mathbf{r} - \frac{5(\mathbf{m}_1 \mathbf{r})(\mathbf{m}_1\mathbf{r})}{r^7}\mathbf{r}\\\\ &= \mathbf{c}_1 + \mathbf{c}_2 + c_3 (\mathbf{r}_2 - \mathbf{r}_1) - c_4 (\mathbf{r}_2 - \mathbf{r}_1) \qquad (5)\\\\ &= \mathbf{c}_1 + \mathbf{c}_2 + \mathbf{r}_1 (c_4 - c_3) + \mathbf{r}_2 (c_3 - c_4)\qquad (6)\\\\ &= \mathbf{c}_1 + \mathbf{c}_2 + \mathbf{c} + d\cdot\mathbf{r}_2, \quad \mathbf{c}_1, \mathbf{c}_2, \mathbf{c}\in \mathbb{R}^2, d\in\mathbb{R} \qquad (7) \end{align}

The introduced variables c1,c2, c, d should be self-explanatory. Now we can define 4 ODEs and solve them numerically with the Runge-Kutta method:

\begin{align} \dot{r}_{2x} &= v_{2x} \qquad (8)\\\\ \dot{r}_{2y} &= v_{2y} \qquad (9)\\\\ \dot{v}_{2x} &= c_{1x} + c_{2x} + c_x + d\cdot r_{2x} \qquad (10)\\\\ \dot{v}_{2y} &= c_{1y} + c_{2y} + c_y + d\cdot r_{2y} \qquad (11)\\\\ \end{align}

What to expect?

An attractive force should act on the magnet m2, so that

$$t \gg 1: \mathbf{r}_2 = \mathbf{r}_1$$

My Results

m2 is repelled instead of being attracted.

Probably I made a mistake in the implementation. But I have already checked it several times and found no error... Does anyone of you know what went wrong

Python 3 Code

import numpy as np
import matplotlib.pyplot as plt

"""
Important quantites that will influence the dynamics

m1:  Magnetic Moment of the dipole that creates a constant magnetic field
m2:  Magnetic Moment of a dust particle
r1:  Position of m1: r1(t) = const
r2:  Position of m2: r2(t) will depend on time
"""

# Constants
m1 = np.array([0, 1], float)
m2 = np.array([0, -1], float)
r1 = np.array([0, 0], float)

# Initial conditions of the of m2 R = [r2_x, r2_y, v2_x, v2_y]
R = np.array([0, 3, 0, 0], float)

def distance(r1, r2, n):
"""
computes the distance ||.||_2 between position r1 and r2
"""

return np.linalg.norm(r1 - r2) ** (n)

def pre_factors(r1, r2, m1, m2):
"""
Computes the prefactors in the equation of motion
d^2 r2 / dt^2 = F = c1 + c2 + c + r2 * d   r2,c1,c2,c € R^2, d € R
"""

r = r2 - r1
l = distance(r1, r2, 5)

c1 = ( np.dot(m1, r) * m2 ) / l
c2 = ( np.dot(m2, r) * m1 ) / l

c3 = np.dot(m1, m2) / l
c4 = ( 5 * np.dot(m1, r) * np.dot(m2, r) ) / distance(r1, r2, 7)

c = (c4 - c3) * r2
d = c3 - c4

return c1, c2, c, d

def f(r1, R, m1, m2):
"""
now we can compute the righthand side of the four 1st order ODEs:
d/dt r2_x = v2_x
d/dt r2_y = v2_y
d/dt v2_x = c1_x + c2_x + c_x + d*r2_x
d/dt v2_y = c1_y + c2_y + c_y + d*r2_y
"""

r2x = R[0]
r2y = R[1]
v2x = R[2]
v2y = R[3]

c1, c2, c, d = pre_factors(r1, R[:2], m1, m2)

f_r2x = v2x
f_r2y = v2y
f_v2x = c1[0] + c2[0] + c[0] + d * r2x
f_v2y = c1[1] + c2[1] + c[1] + d * r2y

return np.array([f_r2x, f_r2y, f_v2x, f_v2y], float)

# Constants for the Simulation
a = 0.0  # t_0
b = 20.0  # t_end
N = 1000 # number of iterations
h = (b-a) / N # time step, e.g. h = (100 - 0) / 1000 = 0.1s

# Create lists to store the computed values
tpoints = np.arange(a,b+h,h)

xpoints, ypoints = [], []
vxpoints, vypoints = [], []

for dt in tpoints:
#print(R)
xpoints.append(R[0])
ypoints.append(R[1])
vxpoints.append(R[2])
vypoints.append(R[3])

# Runge-Kutta-4th Order Method
k1 = dt * f(r1, R, m1, m2)
k2 = dt * f(r1, R + 0.5 * k1, m1, m2)
k3 = dt * f(r1, R + 0.5 * k2, m1, m2)
k4 = dt * f(r1, R + k3, m1, m2)
R += (k1 + 2*k2 * 2*k3 + k4)

plt.style.use('seaborn')
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 3))
axes[0].plot(tpoints, xpoints, label = '$$r_{x_2}(t)$$')
axes[0].plot(tpoints, ypoints, label = '$$r_{y_2}(t)$$')

axes[0].legend()

axes[1].plot(tpoints, vxpoints,'r', label = '$$v_{x_2}(t)$$', alpha = 0.5)
axes[1].plot(tpoints, vypoints, 'k',label = '$$v_{y_2}(t)$$', alpha = 0.5)
axes[1].legend()

fig.suptitle("Results for: $$r_1 = (0,0), r_2 = (0,3), m_1 = (0,1), m_2 = (0,-1)$$")
#plt.savefig("2.pdf" , format='pdf',  bbox_inches="tight")
fig.tight_layout()

• Why do you think the two dipoles should be attracted to each other? It depends on their orientation, in general the force can be attraction, repulsion, or zero. Feb 5, 2021 at 0:04