Nuage d’ions

On considère un nuage de \(N\) ions de masse \(m\) et charge \(q\). On note \(\vec{r_i}\) et \(\vec{v_i}\) la position et la vitesse du ième ion. Les ions sont dans un piège électrostatique. De plus, ils interagissent entre eux par la force de Coulomb.

La force électrostatique dérive d’un potentiel électrostatique :

\[V(\vec{r}) = \frac12 k_x r_x^2 +\frac12 k_y r_y^2 + \frac12 k_z r_z^2\]

La force de Coulomb s’écrira sous la forme:

\[\vec{f_i}(\vec{r_i}, \vec{r_j}) = \kappa \frac{q^2}{|\vec{r_i} - \vec{r_j}|^3}(\vec{r_i} - \vec{r_j})\]

Pour simplifier et avoir une meilleure représentation graphique, on suprimera totalement la dimension \(z\).

On utilisera des unités adimensionnées avec des constantes de l’ordre de 1. Par défaut on prendra : \(m=9\), \(q=1\), \(k_x=1\), \(k_y=1.3\), \(\kappa=1\).

L’état des \(N\) particules à un instant donné est représenté par 4 tableaux numpy de taille \(N\): r_x, r_y, v_x, v_y. Toutes les fonctions seront écrites avec ces variables. On ne regroupera les 4 tableaux que dans la fonction f(t, y) qui sera utilisée par solve_ivp.

Le tableau y est défini en rassemblant les 4 tableaux numpy. On définit alors les fonctions join et split (voir ci-dessous)

Les paramètres seront des constantes globales.

  1. Écrire la fonction force_piege(r_x, r_y) qui renvoie la force dérivant du potentiel (f_x et f_y).

  2. Écrire la fonction force_coulomb(r_x, r_y) qui renvoie \(f_x\) et \(f_y\), la force de Coulomb. On ne cherchera pas à éviter les boucles for dans cette fonction.

  3. Écrire la fonction f(t, y) qui définit la dynamique du problème.

  4. On condidère un nuage avec une distribution initiale donnée par une loi de Maxwell-Boltzmann (avec \(k_B=1\) et \(T=1\)) pour les particules sans interaction. Calculer jusqu’au temps \(T=20\) l’évolution de la position des particules.

  5. Vérifier que l’énergie totale est conservée. On utilisera les fonctions ci-dessous.

  6. Faire une animation (voir exemple de code ci dessous)

  7. On rajoute une force de dissipation, selon l’axe \(x\) :

    \[ \vec F_{i,x} = -\alpha v_{i, x} \]

    Simuler l’experience en prenant \(\alpha = 0.1\) et \(N=20\). Que se passe-t-il au temps long ?

  8. Utiliser le decorateur jit de numba pour accélérer le code. Faire une comparaison de vitesse.

  9. On considère un nuage contenant deux types de particules (masse \(m_1\) et \(m_2\)) et de même charge \(q=1\). Seuls les particules de type \(1\) sont soumises à la force dissipative. Simuler l’expérience en prenant les mêmes paramètres qu’à la question 7. On prendra par exemple \(m_1=9\), \(N_1=20\) et \(m_2=2\) et \(N_2=5\)

# Paramètres
m = 9; q = 1
k_x, k_y = 1, 1.3
kappa=1
k_B = 1
T = 1

Quelques fonctions pour vous aider :

  • Fonction pour passer de 4 tableaux à un seul et réciproquement

def join(r_x, r_y, v_x, v_y):
    return np.concatenate((r_x, r_y, v_x, v_y))

def split(y):
    N = len(y)//4
    return y[:N], y[N:2*N], y[2*N:3*N], y[3*N:4*N]

  • Calcul de l’énergie

def energie_cinetique(r_x, r_y, v_x, v_y):
    return np.sum(.5*m*v_x**2 + .5*m*v_y**2)

def energie_piege(r_x, r_y):
    return np.sum(k_x*r_x**2/2 + k_y*r_y**2/2)

def energie_coulomb(r_x, r_y):
    N = len(r_x)
    total = 0
    for i in range(N-1):
        for j in range(i+1, N):
            d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2
            total += kappa*q**2/np.sqrt(d2)
    return total

def energie_totale(r_x, r_y, v_x, v_y):
    return (energie_coulomb(r_x, r_y) +
            energie_piege(r_x, r_y) +
            energie_cinetique(r_x, r_y, v_x, v_y))

from IPython.display import HTML
from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots()
l, = ax.plot(split(res.y[:,0])[0], split(res.y[:,0])[1], 'o')
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)

def animate(i):
    r_x, r_y, v_x, v_y = split(res.y[:,i])
    l.set_data(r_x, r_y)

ani = FuncAnimation(fig, animate, frames=len(res.t))


HTML(ani.to_jshtml())