10  Microscope à force atomique

Un microscope à force atomique (atomic force microscope ou AFM) est un appareil servant à mesurer diverses propriétés physico-chimiques de surfaces (solides ou recouvertes d’un film). Son utilisation principale est la mesure de la rugosité de surfaces, mais il est également possible de mesurer le module d’Young d’un solide en divers points de sa surface. Cette technique demande d’étudier la dynamique vibratoire de l’AFM.

Le principe de fonctionnement d’un AFM est simple: il s’agit d’une poutre console (encastrée-libre) dotée à l’extrémité libre d’une pointe qui rentre en contact avec la surface solide. La rotation de la tête par rapport à la position horizontale est mesurée à l’aide d’un miroir posé au dessus de la tête qui dévie un laser. La position du laser dévié permet de déduire la rotation, et les caractéristiques de flexion de la poutre permettent de reconstruire le déplacement vertical de la tête et la force que subit celle-ci.

Fonctionnement d’un AFM. Modifié à partir du diagrame de la page wikipedia

Dans cet exercice, on va discrétiser la poutre en 4 segments rigides reliés entre eux par des ressorts de flexion de raideur \(K\). L’inertie de la poutre est représentée par 3 masses \(m\) et une masse \(10m\) en bout de poutre qui représente la pointe. Cette dernière est liée par un ressort \(k\) au substrat. La raideur de ce ressort est considérée comme une inconnue du système puisqu’elle dépend du module d’Young du solide que l’on cherche à mesurer.

Schématisation d’un AFM discrétisé

La position verticale de chaque masse est notée \(z_i\), \(i = 1, ..., 4\), et \(z_0\) est la position verticale de l’encastrement, contrôlée par l’opérateur de la machine. Il est ainsi possible de mettre en vibration la poutre à une fréquence donnée (\(z_0(t) = Z_0\cos\Omega t\)) et de mesurer la réponse dynamique. C’est le mécanisme dont on se servira pour trouver la valeur de \(k\).

Puisqu’il est difficile de faire les calculs à la main avec 4 degrés de liberté, on utilise l’outil numérique. On se servira donc de Numpy pour le calcul matriciel et Scipy pour la diagonalisation et le calcul de matrices inverses. Les fonctions utiles à ce TD sont donc:

Notebook live

Pour faire l’exercice, vous devez vous rendre ici, depuis n’importe quel navigateur web, il n’y a pas besoin d’installer Python.

Il sera peut-être nécessaire d’aller dans le menu Help > Clear Browser Data pour recharger les notebooks.

10.1 Dynamique et modes propres

On commence par poser les équations de la dynamique du système discret.

Montrer que l’énergie cinétique et la matrice d’inertie sont respectivement:

\[ T = \frac{1}{2}m\left(\dot z_1^2 + \dot z_2^2 + \dot z_3^2 + 10\dot z_4^2\right),\quad \mathbf M = m\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 10\end{pmatrix} \]

Sous l’hypothèse des petits angles, montrer que l’énergie potentielle et la matrice de raideur sont respectivement:

\[\begin{align*} U &= \frac{1}{2}\frac{K}{L^2}\left((z_1 - z_0)^2 + (z_2 - 2 z_1 + z_0)^2 + (z_3 - 2 z_2 + z_1)^2 + (z_4 - 2 z_3 + z_2)^2\right)+ \frac{1}{2}kz_4^2, \\ \mathbf K &= \frac{K}{L^2} \begin{pmatrix} 6 & -4 & 1 & 0\\ -4 & 6 & -4 & 1\\ 1 & -4 & 5 & -2\\ 0 & 1 & -2 & 1 +\frac{kL^2}{K} \end{pmatrix} \end{align*}\]

Note

Le ressort de rotation entre les segments de poutres \((i-1, i)\) et \((i, i+1)\) ne se déforme que quand les rotations des deux segments sont différentes, soit:

\[ U_i = \frac{1}{2}K\Delta \theta_i^2 = \frac{1}{2}K(\theta_{i+1} - \theta_i)^2 \]

Où \(\theta_i\) est l’angle que fait le segment \((i-1, i)\) avec l’horizontale (respectivement \(\theta_{i+1}\) pour le segment \((i, i+1)\)), qu’il faut exprimer en fonction des positions \(z_i\) et \(z_{i-1}\) (respectivement \(z_{i+1}\) et \(z_i\)).

Cinématique d’une rotule

À partir de l’équation d’Euler-Lagrange, déduire que le vecteur des efforts généralisés s’écrit :

\[ \vec Q(t) = \frac{K}{L^2}\begin{pmatrix} 3 \\ -1 \\ 0 \\ 0 \end{pmatrix} z_0(t) \]

On importe les paquets nécessaire à la suite de l’exercice.

import numpy as np
import scipy as sc
import matplotlib.pyplot as plt

Implémenter les fonctions deux fonctions suivantes (vous pourrez vous aider des fonctions np.array et np.eye):

def mass_matrix(m):
    """Retourne la matrice de masse en fonction de m"""
    
    return ...

def stiffness_matrix(K, L, k):
    """Retourne la matrice de raideur en fonction de K, L et k"""

    return ...

Implémenter la fonction suivante à l’aide de sc.linalg.eigh :

def diagonalize(K_mat, M_mat):
    """Retourne les pulsation propres et la matrice des modes propres"""

    return ..., ...

Avec les valeurs de \(K, L, k\) et \(m\) ci-dessous, appelez les fonctions que vous avez implémentées pour calculer les pulsations propres et les modes propres. Vérifiez que les égalités suivantes sont vraies (à la précision machine) :

\[ \mathbf X^T \mathbf M \mathbf X = \mathbf I,\quad \mathbf X^T \mathbf K \mathbf X = \mathrm{diag}(\{\omega_i^2\}),\quad \text{où}\quad \mathrm{diag}(\{v_i\}) = \begin{pmatrix} v_1 & 0 & 0 & 0 \\ 0 & v_2 & 0 & 0\\ 0 & 0 & v_3 & 0\\ 0 & 0 & 0 & v_4 \end{pmatrix}\]

K = 1.
L = 1.
k = 0.
m = 1.
Note

On peut visuellement comparer les valeurs pour s’assurer de l’égalité, ou utiliser la fonction np.testing.assert_almost_equal(A, B) pour tester l’égalité entre les matrices A et B.

À l’aide de la méthode de Rayleigh, vérifiez la première pulsation propre.

Affichez sur un graphe la déformée des 4 modes.

10.2 Influence du substrat

10.2.1 Dimension des quantités

Dans les questions ci-dessus, on a supposé que \(k = 0\). En regardant l’expression de la matrice de raideur, l’effet de \(k\) dans cette dernière est donné par un nombre sans dimension \(kL^2/K\).

De même, on remarque en posant l’équation de la dynamique que \(\omega_0 = \sqrt{\frac{K}{L^2m}}\) est l’unité naturelle des pulsations. Choisir \(K = 1\), \(L = 1\) et \(m = 1\) revient à résoudre le problème adimensioné, et la valeur choisie de \(k\) est le nombre sans dimension \(k L^2/K\).

10.2.2 Effet de \(k\) sur les fréquences des modes

On cherche à quantifier l’influence du substrat sur les fréquences de vibration des deux premiers modes. Pour cela, on a besoin d’une fonction qui, pour une valeur de \(k\), nous donne la fréquence propre du mode désiré.

@np.vectorize
def mode_freq(k, mode_number):
    """
    Retourne la fréquence d'oscillation du mode `mode_number`
    pour une raideur `k` du substrat
    """
    K_mat = ... # en fonction de k (K = L = 1)
    M_mat = ...
    ...
    return ...

Pour des valeurs de \(kL^2/K\) de 0 à 20 (np.linspace), affichez sur un graphe les fréquences des deux premiers modes (à l’aide de la fonction ci-dessus). Sur quelle plage de valeur les modes sont-ils sensibles à une variation de \(k\) ?

10.3 Réponse libre

Exprimer les coordonnées modales \(p_i(t)\) solutions des équations modales:

\[\ddot p_i + \omega_i^2 p_i = 0\]

Déterminez les constantes inconnues à partir des conditions initiales suivantes :

\[ \mathbf z(0) = 0,\quad\dot{\mathbf{z}}(0) = \mathbf v = -10^{-2}\begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} \]

Calculer \(z_{i}\) en fonction du temps sur deux périodes du mode 1.

def displacement(t, X, omega, B):
    """Retourne un tableau (t.size * 4) contenant les valeurs de z_i en fonction de t"""
    # Étape 1. Calculer les coordonnés modales p(t)
    # Étape 2. Repasser dans la bases des degrés de liberté z(t)
    return ...

Calculer et afficher l’angle \(\theta_4\) (l’angle du mirroir réfléchissant le laser) en fonction du temps.

Bonus traitement de signal : afficher le spectre du signal ci-dessus. Vérifiez que les fréquences dominantes sont bien les fréquences propres des 4 modes.

10.4 Réponse forcée

On reprend le vecteur des efforts généralisés:

\[ \mathbf Q(t) = \frac{K}{L^2}\begin{pmatrix} 3 \\ -1 \\ 0 \\ 0 \end{pmatrix} z_0(t) = \mathbf E z_0(t) \]

Pour un mouvement harmonique de fréquence \(\Omega\) (\(z_0(t) = Z_0\cos(\Omega t)\)), calculez l’amplitude de la réponse forcée des coordonnées modales.

E = K / L**2 * np.array([3, -1, 0, 0])

def forced_response_modal(Omega, Z0, E, omega, X):
    """Retourne l'amplitude de la réponse forcée pour les modes"""
    return ...

On repasse ensuite dans la base des degrés de liberté, calculez l’amplitude de la réponse forcée pour les \(z_i\).

def forced_response_displacement(Omega, Z0, E, omega, X):
    """Retourne l'amplitude de la réponse forcée pour les z_i"""
    # Il faut appeler la fonction ci-dessus !
    return ...

À l’aide des fonctions ci-dessus, calculez la fonction de transfert de \(\theta_4\) en fonction de \(k\) et de \(\Omega\).

@np.vectorize
def transfer_theta_4(Omega, k, E):
    """Pour une valeur de Ω et de k, calcule l'amplitude d'oscillation de la rotation du mirroir"""
    omega, X = ... # en fonction de k
    Z = ...
    theta_4 = ...
    return theta_4

Pour les valeurs de \(\Omega\) et \(k\) ci-dessous, affichez la fonction de transfert de \(\theta_4\).

Omega_values = np.geomspace(1e-2, 1e1, 300)
k_values = [0, 0.05, 5, 10]

Chaque masse a un seul degré de liberté de translation, l’énergie cinétique est triviale et \(M_{ij} = \frac{\partial^2 T}{\partial \dot z_i\partial \dot z_j}\).

L’énergie des ressorts 1, 2, 3 peut être écrite avec la formule ci-dessus, il faut juste exprimer l’angle \(\theta_i\) dans l’hypothèse des petits angles:

\[ \theta_i = \frac{z_i - z_{i-1}}{L},\quad\text{et}\quad \theta_{i+1}-\theta_i = \frac{1}{L}(z_{i+1}-2z_i+z_{i-1}) \]

L’énergie potentielle du ressort 0 (celui de l’encastrement) dépend uniquement de l’angle \(\theta_1\), et l’énergie du ressort du substrat dépend uniquement de \(z_4\), on a donc:

\[\begin{align*} U & = \frac{1}{2}K(\theta_1^2 + (\theta_2 - \theta_1)^2 + (\theta_3 - \theta_2)^2 + (\theta_4 - \theta_3)^2) + \frac{1}{2}kz_4^2\\ & = \frac{1}{2}\frac{K}{L^2}\left((z_1 - z_0)^2 + (z_2 - 2 z_1 + z_0)^2 + (z_3 - 2 z_2 + z_1)^2 + (z_4 - 2 z_3 + z_2)^2\right) \\ &+ \frac{1}{2}kz_4^2 \end{align*}\]

Attention, \(z_0\) n’est pas un degré de liberté et n’intervient pas dans le calcul de \(\mathbf K\), que l’on détermine par dérivée:

\[ K_{ij} = \frac{\partial^2 U}{\partial z_i\partial z_j}. \]

On pourra vérifier que tous les termes diagonaux de \(\mathbf K\) sont positifs et que \(\mathbf K\) est symétrique. Si ce n’est pas le cas il y a une erreur de calcul.

En calculant la première dérivée de \(U\), on fait apparaître des termes en \(z_0\). En isolant ces termes dans l’équation de la dynamique \(\mathbf M \ddot{\vec z} + \mathbf K{\vec z} = \vec Q(t)\) on peut identifier \(\vec Q\). La procédure est identique à l’exercice du salaire de la peur.

def mass_matrix(m):
    """Retourne la matrice de masse en fonction de m"""
    M_mat = m * np.eye(4)
    M_mat[3, 3] *= 10
    return M_mat

def stiffness_matrix(K, L, k):
    """Retourne la matrice de raideur en fonction de K, L et k"""
    K_mat = K/L**2 * np.array([
        [6, -4, 1, 0],
        [-4, 6, -4, 1],
        [1, -4, 5, -2],
        [0, 1, -2, 1 + k * L**2 / K],
    ])

    return K_mat
def diagonalize(K_mat, M_mat):
    """Retourne les pulsation propres et la matrice des modes propres"""
    omega2, X = sc.linalg.eigh(K_mat, M_mat)
    return np.sqrt(omega2), X
K_mat = stiffness_matrix(K, L, k)
M_mat = mass_matrix(m)

omega, X = diagonalize(K_mat, M_mat)

np.testing.assert_almost_equal(X.T @ M_mat @ X, np.eye(4))
np.testing.assert_almost_equal(X.T @ K_mat @ X, np.diag(omega**2))
# S'il n'y a pas d'erreur alors tout va bien :)

Pour utiliser le quotient de Rayleigh, il faut la déformation statique de la structure due à une distribution de forces proportionnelles aux masses.

\[\mathbf f = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 10\end{pmatrix}\]

On calcule ensuite la solution de \(\mathbf K\mathbf z_\mathrm{stat} = \mathbf f\), et on applique le quotient de Rayleigh:

\[\omega_\mathrm{R}^2 = \frac{\langle \mathbf z_\mathrm{stat}, \mathbf K\mathbf z_\mathrm{stat}\rangle}{\langle \mathbf z_\mathrm{stat}, \mathbf M\mathbf z_\mathrm{stat}\rangle}\]

f = np.array([1, 1, 1, 10])
z_stat = sc.linalg.solve(K_mat, f)

omega_R = np.sqrt(
    np.dot(z_stat, K_mat @ z_stat)
    / np.dot(z_stat, M_mat @ z_stat)
)

print(f"""Rayleigh: ω = {omega_R:.3e}
Exact:    ω = {omega[0]:.3e}
Erreur:   {np.abs(omega_R - omega[0])/omega[0]*100:.3f}%""")
Rayleigh: ω = 5.609e-02
Exact:    ω = 5.608e-02
Erreur:   0.023%
fig, ax = plt.subplots()

position = np.arange(5) # entiers de 0 à 4

for mode_num, (ω, x) in enumerate(zip(omega, X.T)):
    x = np.concatenate([[0], x]) # on rajoute z_0 = 0 pour l'encastrement
    ax.plot(position, x,
            marker='o', markerfacecolor='w',
            label=f'Mode {mode_num+1} (ω = {ω:.2e})')

ax.set(xticks=position, xlabel='Position', ylabel='Déplacement vertical')
ax.grid()
ax.legend()
plt.show()

On remarque que les modes 2, 3, 4 ne font presque pas osciller la pointe de l’AFM comparé au mode fondamental.

@np.vectorize
def mode_freq(k, mode_number):
    """
    Retourne la fréquence d'oscillation du mode `mode_number`
    pour une raideur `k` du substrat
    """
    K_mat = stiffness_matrix(1, 1, k)
    M_mat = mass_matrix(1)
    omega, _ = diagonalize(K_mat, M_mat)
    return omega[mode_number] / (2 * np.pi)
fig, ax = plt.subplots()

k_vals = np.linspace(0, 20, 100)
freqs_1 = mode_freq(k_vals, 0)
freqs_2 = mode_freq(k_vals, 1)

ax.plot(k_vals, freqs_1, label='Mode 1')
ax.plot(k_vals, freqs_2, label='Mode 2')
ax.grid()
ax.set(xlabel='Raideur $k L^2 / K$',
       ylabel='Fréquence $f / \\sqrt{K/(L^2m)}$')
ax.legend()
plt.show()

Le mode 1 est le plus sensible entre les valeurs de 0 et 5, au-delà une variation de \(k\) n’affecte pas la fréquence. Inversement, le mode 2 est relativement insensible aux faibles valeurs de \(k\), et la fréquence est affectée au-delà de 5. Le dispositif est donc efficace pour mesurer des valeurs de \(kL^2/K\) entre 0 et 20 (au moins).

\[p_i(t) = A_i\cos(\omega_it) + B_i \sin(\omega_i t)\]

On a d’abord \[ \mathbf z(0) = 0 \Leftrightarrow \mathbf{X}\mathbf p(0) = 0 \Leftrightarrow \mathbf p(0) = 0 \Leftrightarrow A_i = 0\,\forall i\]

Puis

\[ \dot{\mathbf{z}}(0) = 0 \Leftrightarrow \mathbf X \mathbf{\dot p}(0) = \mathbf v \Leftrightarrow \mathbf{\dot p}(0) = \mathbf X^{-1}\mathbf v\]

Puisque \(A_i = 0\), on peut exprimer \(\mathbf{\dot p}(0) = \mathrm{diag}(\{\omega_i\})\mathbf B\), on obtient ainsi:

\[ \mathbf B = \mathrm{diag}(\{1/\omega_i\})\mathbf X^{-1}\mathbf v \]

v = -1e-2 * np.array([1., 1., 1., 1.])
B = np.diag(1 / omega) @ sc.linalg.inv(X) @ v
def displacement(t, X, omega, B):
    """Retourne un tableau (t.size * 4) contenant les valeurs de z_i en fonction de t"""

    # lignes = temps, colonnes = mode
    coord_modales = np.zeros([t.size, omega.size]) 
    for i in range(omega.size):
        # on selectionne la "i"-eme colonne du tableau de p(t)
        coord_modales[:, i] = B[i] * np.sin(omega[i] * t)

    z = np.zeros_like(coord_modales)
    for i in range(t.size):
        # on repasse dans la base des DDL pour chaque temps
        z[i, :] = X @ coord_modales[i, :]
    return z
tmax = 2 * 2 * np.pi / omega[0]
Ns = 500
t = np.linspace(0, tmax, Ns)

Z = displacement(t, X, omega, B)
for i, z in enumerate(Z.T):
    plt.plot(t, z, label=f'$z_{i+1}$')
plt.grid()
plt.gca().set(xlabel='Temps', ylabel='Déplacement')
plt.legend()
plt.show()

On a \(\theta_4 = \frac{z_4 - z_3}{L}\).

mirror_angle = (Z[:, 3] - Z[:, 2]) / L
plt.plot(t, mirror_angle * 360 / (2 * np.pi))
plt.grid()
plt.gca().set(xlabel='Temps', ylabel='Angle du mirroir en degrés')
plt.show()

plt.psd(mirror_angle, Fs=Ns/tmax, NFFT=Ns)
for ω in omega:
    plt.axvline(ω / (2 * np.pi), color='k', ls='--', zorder=-1)
plt.show()

On doit effectuer un changement de base de l’équation de la dynamique pour calculer une solution particulière dans la base des coordonnées modales:

\[\begin{align*} \mathbf M \ddot{\mathbf z} + \mathbf K\mathbf z = \mathbf Q & \Leftrightarrow \mathbf X^T \mathbf M \mathbf X \mathbf{\ddot p} + \mathbf X^T \mathbf K \mathbf X \mathbf{\ddot p} = \mathbf X^T\mathbf Q \\ & \Leftrightarrow \ddot{\mathbf p} + \mathrm{diag}(\{\omega_i^2\})\mathbf p = \mathbf X^T \mathbf Q \end{align*}\]

Puisque \(z_0(t) = Z_0\cos(\Omega t)\), on cherche des solutions de la forme \(p_i(t) = P_i\cos(\Omega t)\). En remplaçant dans l’équation ci-dessus, on obtient:

\[\begin{align*} -\Omega^2\mathbf P \cos(\Omega t) + \mathrm{diag}(\{\omega_i^2\})\mathbf P\cos(\Omega t) & = \mathbf X^T \mathbf E Z_0\cos(\Omega t) \\ \Rightarrow (\mathrm{diag}(\{\omega_i^2\}) - \Omega^2\mathbf I)\mathbf P & = \mathbf X^T \mathbf E Z_0 \\ \Rightarrow \mathbf P & = (\mathrm{diag}(\{\omega_i^2\}) - \Omega^2\mathbf I)^{-1} \mathbf X^T \mathbf E Z_0\\ \Rightarrow \mathbf P & = \mathrm{diag}\left(\left\{\frac{1}{\omega_i^2 - \Omega^2}\right\}\right) \mathbf X^T \mathbf E Z_0 \\ \end{align*}\]

def forced_response_modal(Omega, Z0, E, omega, X):
    """Retourne l'amplitude de la réponse forcée pour les modes"""
    Q_modal = X.T @ E * Z0
    amplitude_modal = np.zeros(omega.size)
    for i in range(omega.size):
        amplitude_modal[i] = Q_modal[i] / (omega[i]**2 - Omega**2)
    return amplitude_modal

L’amplitude s’exprime \(\mathbf Z = \mathbf X \mathbf P\).

def forced_response_displacement(Omega, Z0, E, omega, X):
    p = forced_response_modal(Omega, Z0, E, omega, X)
    return X @ p
@np.vectorize
def transfer_theta_4(Omega, k):
    """Pour une valeur de Ω et de k, calcule l'amplitude d'oscillation de la rotation du mirroir"""
    omega, X = diagonalize(stiffness_matrix(1, 1, k), mass_matrix(1))
    Z = forced_response_displacement(Omega, 1, E, omega, X)
    theta_4 = np.abs(Z[3] - Z[2]) / L
    return theta_4
for k in k_values:
    T = transfer_theta_4(Omega_values, k)
    plt.loglog(Omega_values, T, label=f'k = {k}')

plt.gca().set(xlabel='$\\Omega$', ylabel='Transfert $\\theta_4$')
plt.legend()
plt.gcf().set_size_inches(4, 4)
plt.show()

On remarque que les fréquences des résonances des modes 3 et 4 ne sont pas affectées par la valeur de \(k\), il faut donc observer les pics de résonance des deux premiers modes.

On remarque aussi une fréquence caractéristique d’excitation pour laquelle le miroir ne vibre pas. Pour les faibles valeurs de \(k\) cette fréquence est plus petite que la première résonance. En augmentant \(k\) elle finit par dépasser la première résonance.