2
$\begingroup$

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

enter image description here enter image description here

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()
$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$ – Maxim Umansky Feb 5 at 0:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.