INTERaction GAz-Grain

La structure carbonée

Les poussières "cibles" de l'atome d'hydrogène seront choisies comme des structures carbonées régulières. Elles sont a priori tridimensionnelles mais cette appliquette ne présentera que le cas 1D, à savoir une chaîne de n carbone fixée ou non à un bâti.

On peut s'interroger sur la pertinence d'une telle modélisation qui pourrait être qualifiée d'un peu trop simpliste. Mais cela n'est pas complètement le cas. Bien évidemment, les équations et la programmation sont plus commodes à 1D mais ces chaînes qui semblent utopiques (la nature n'aime pas trop en général ce qui n'est pas 3D, la notion de "couche" par exemple n'est qu'une vue de l'esprit, simple modèle limite des physiciens qui conduit souvent d'ailleurs à des discontinuités) ne le sont pas. Il existe dans la nature de véritables chaînes carbonées : les Cyanopolyynes [pour plus de renseignements voir http://iopscience.iop.org/0004-637X/505/1/278 (en anglais) qui discute de l'existence de telles molécules dans le milieu interstellaire.] D'autres points, liés à la modélisation 1D, sont cependant difficiles à justifier physiquement : le fait que l'atome d'hydrogène vienne percuter la chaîne exactement dans l'axe par exemple...

Quant à l'existence de 2 modèles chaînes : liée au bâti ou libre, elle se justifie assez bien. L'existence des Cyanopolyynes justifie déjà l'utilisation de la chaîne libre ; en ce qui concerne une chaîne liée à un bâti il suffit de considérer le bâti comme une très grosse molécule (son inertie très grande devant celle de l'atome d'hydrogène incident lui assurera de ne pas trop bouger pendant l'interaction, c'est-à-dire de jouer finalement le rôle d'une sorte de... bâti bien solide), la chaîne quant à elle sera tel un "cheveu" avec lequel viendra interagir l'atome d'hydrogène.

Les 2 figures qui suivent représentent ces 2 modèles. La position des atomes de carbone est repérée par les x_n tandis que l'atome d'hydrogène incident sera associé à y.

Modèle "Chaine libre"

Modèle "Chaine liée à un bâti"

Les interactions entre les atomes de carbone (de masses identiques M) sont modélisées par de simples ressorts de raideur k. La distance inter-atomes sera notée d_CC. Le potentiel d'interaction de la chaîne de ressorts aura donc pour forme $U_{CC}(x_1,...,x_n)=\frac{1}{2}k\sum_{i=1}^{n-1} (x_{i+1}-x_i-d_{cc})^2$ pour la chaîne libre -- à laquelle s'ajoutera $\frac{1}{2}k(x_1-d_{cc})^2$ si l'on considère le modèle avec bâti.

Quant aux valeurs numériques des paramètres du système on a :

M, la masse du carbone qui vaut environ 12 masses atomiques.
k, raideur du ressort qui est calculée à partir de la formule suivante : $k=4 \pi^2 M f_{vib}^2$ où f_vib vaut environ 10¹²Hz qui est la fréquence typique de vibration des liaisons inter-carbones.
d_CC, distance au repos entre les atomes qui est de l'ordre de 150 pm.

Les interactions Hydrogène-Carbone

On doit désormais s'intéresser à l'interaction que va subir l'hydrogène à l'approche de l'atome de carbone. On simplifie premièrement le modèle en supposant que l'hydrogène n'interagira qu'avec le dernier carbone de la chaîne (x_n).

Le potentiel qui assumera l'interaction sera ensuite celui de Lennard-Jones. Attractif aux longues distances (interactions classiques de type Coulomb), il est répulsif à très petite distance (résultat provenant de la physique quantique traduisant le principe d'exclusion de Pauli). Il existe donc une distance où ce potentiel est minimal r₀qui correspond à la distance d'équilibre entre un atome de carbone et un d'hydrogène. Un second paramètre E règle la profondeur de ce "puits" d'attraction.
$U_{CH}(x_n,y)=2E \left[\frac{1}{2}\left(\frac{r_0}{y-x_n}\right)^{12}-\left(\frac{r_0}{y-x_n}\right)^6\right]$

Potentiel de Lennard-Jones pour r0 = 1 nm et E = 0,07 eV

Pour les valeurs numériques typiques on a :

E, profondeur du puits d'attraction du potentiel de Lennard-Jones qui vaut aux alentours de 0.07 eV.
r₀, distance d'équilibre pour le potentiel de Lennard-Jones, qui est à peine plus faible que d_CC c'est à dire 0.1 nm.

Les équations du mouvement

Les potentiels définis, il est désormais possible d'écrire les équations de la dynamique. On utilisera les notations et équations de la mécanique analytique. Voici la trame générale des calculs qui permettront d'aboutir à la formulation et à la résolution du problème physique qui nous intéresse.

Dédimensionnement du système
Écriture du lagrangien dédimensionné
Passage aux variables hamiltoniennes
Équations d'Hamilton
Résolution des équations d'Hamilton
Retour dans l'espace réel et redimensionnement.

Energies cinétiques, potentielles et dédimensionnement

L'énergie cinétique du système correspond à la somme des énergies cinétiques de chacun des atomes :
$T(\dot{x_1},...,\dot{x_n},\dot{y})=\frac{1}{2}M \sum_{i=1}^n \dot{x_i}^2 + \frac{1}{2}m\dot{y}^2$ où m est la masse de l'atome d'hydrogène (1 masse atomique).

Le potentiel total est la somme des potentiels élastiques de la chaîne de carbone et du potentiel de Lennard-Jones. Noter que l'on se placera pour la suite des calculs dans le cas du modèle "Chaîne libre". (Le calcul dans le cas "Chaîne avec bâti" ne diffère que très peu : voir plus haut la forme du potentiel lorsque le ressort "bâti-1^er atome" est pris en compte.)
$U(x_1,...,x_n,y)= \frac{1}{2}k\sum_{i=1}^{n-1} (x_{i+1}-x_i-d_{cc})^2+2E \left[\frac{1}{2}\left(\frac{r_0}{y-x_n}\right)^{12}-\left(\frac{r_0}{y-x_n}\right)^6\right]$

Ces énergies cinétiques et potentielles sont dimensionnées. Un léger changement de variables est souvent pratique pour ramener ces quantités à des ordres de grandeur plus facilement maniables. On utilise pour cela les données du problème comme "grandeurs typiques". Les changements concernent ici le temps et les longueurs.
$\centering X_i=\frac{x_i}{r_0} \qquad Y=\frac{y}{r_0} \qquad \tau = \sqrt{\frac{k}{M}} t$

Ces changements disent simplement que la règle unité pour mesurer les distances fera une longueur de r₀. Pour le temps, la pulsation de référence coïncidera précisément avec celle typique d'oscillation de la chaîne de carbone. Pour résumer on va manipuler des quantités sans dimension entre 10^-1 et 10.

Lagrangien

On récupère alors un Lagrangien (L=T-U) de la forme :
$L \propto \frac{1}{2} \left[ \sum_{i=1}^n \dot{X}_i^2 + \mu \dot{Y}^2 \right] - V_0\left[\frac{1}{2}\left(\frac{1}{X_n-Y}\right)^{12}-\left(\frac{1}{X_n-Y}\right)^6\right]-\frac{1}{2} \sum_{i=1}^{n-1} (X_{i+1}-X_i-D_{cc})^2$ avec $\mu = \frac{m}{M} \qquad V_0=\frac{2E}{ k r_0^2} \qquad D_{cc}=\frac{d_{cc}}{r_0}$

Il est ensuite nécessaire de faire apparaître les variables physiques les plus pertinentes. Par exemple dans le potentiel de Lennard Jones, deux variables sont visibles : X_n et Y alors que c'est toujours la quantité Y - X_n dont dépend véritablement le potentiel. Chose tout à fait logique puisque ce potentiel est fonction de la distance entre les atomes C et H. Globalement on va effectuer : q_i=X_i

pour

et $q_{n+1}=Y-X_n$ . Les quantités dérivées seront transformées de la même manière.

Hamiltonien

Le précédent changement de variable n'est pas le dernier. Les variables hamiltoniennes ne sont en effet pas $(q_i,\dot{q_i})$ mais (q_i,p_i)

qu'on appelle positions, impulsions. On calcule les impulsions p_i

par la formule $p_i=\partial L/ \partial \dot{q_i}$ qui traduit le fait que l'impulsion est la variable conjuguée de la position.

On a alors un hamiltonien (sorte d'énergie mécanique dans les cas les plus simples comme ici) H=T+U de la forme :
$H \propto \frac{1}{2} \sum_{i=1}^{n-1} p_i^2 + \frac{1}{2} (p_n-p_{n+1})^2+\frac{1}{2\mu}p_{n+1}^2 + V_0\left[\frac{1}{2q_{n+1}^{12}}-\frac{1}{q_{n+1}^6}\right]+\frac{1}{2} \sum_{i=1}^{n-1} (q_{i+1}-q_i-D_{cc})^2$

Équations d'Hamilton

On peut finalement écrire l'équivalent du principe fondamental de la dynamique dans l'espace d'Hamilton, résumé par les célèbres équations d'Hamilton :
$\dot{q_i}=\frac{\partial H}{\partial p_i} \qquad \dot{p_i}=-\frac{\partial H}{\partial q_i}$

Ce qui donne, lorsqu'on les écrit pour tout i, les équations matricielles suivantes (attention elles ne sont valables que pour le modèle de la chaîne libre) :

Résolution des équations d'Hamilton

Il ne reste donc qu'à résoudre ces équations. Elles peuvent sembler un peu barbares mais en réalité elles le sont bien peu car elles sont linéaires pour la plupart des inconnues, la notation matricielle d'ailleurs le prouve et si l'on omet l'équation qui traduit l'interaction entre l'hydrogène et le carbone : $\dot{p}_{n+1}=f_{nl}(q_{n+1})$ (avec nl pour Non Linéaire) le système peut même être résolu analytiquement. Malheureusement c'est précisément cette interaction qui nous intéresse, impossible donc de l'oublier.

Le module qui va résoudre ces équations est un intégrateur numérique (son fonctionnement sera décrit dans la section algorithmique). Celui-ci étant censé fonctionner pour n'importe quelle équation différentielle, il faut encore réduire la forme des équations du mouvement à une forme ultimement épurée à savoir : $F(X)=\dot{X}$ .

En fait tout est histoire de notations, il suffit en effet de définir F et X de la sorte :

Conditions initiales

Il manque cependant une chose essentielle pour résoudre les précédentes équations différentielles. Ce sont les conditions initiales ! Mon prof de maths lors du cours dédié au théorème de Cauchy-Lipschitz, disait que pour jouer à la belote, il fallait deux choses : d'abord en connaître les règles et aussi, "et surtout" insistait-il, que chacun des 4 joueurs ait sa huitaine de cartes en main. Les règles, ce sont les équations de la physique et plus généralement les équations différentielles. La distribution des cartes, quant à elle, assure la définition d'un état initial. Ces deux choses, l'une dans l'autre, permettent ensuite à la partie de se dérouler ou autrement dit au système physique d'évoluer. (Les rigoristes les plus acharnés douteront peut-être du parallèle puisque qu'à un jeu précis donné, plusieurs stratégies peuvent s'imposer ce qui rompt ainsi l'unicité de la solution qu'est censé offrir Cauchy-Lipschitz. Mais passons, c'était sur le rôle indispensable des conditions initiales qu'il fallait s'attarder.)

Il faut donc connaître ici (pour renseigner l'intégrateur numérique) position et vitesse de chacun des atomes du système à l'instant t=0.

Remarque : Le but premier de notre appliquette était de simuler la projection d'une particule d'hydrogène de vitesse modulable sur une chaîne carbonée qui n'était pas au repos. Comprendre ici que le milieu interstellaire n'est pas complètement calme. Un rayonnement électromagnétique y est toujours présent. Et celui-ci doit en principe exciter la chaîne dans un certain mode de vibration. La chaîne peut également rayonner lorsqu'elle est trop excitée. Il en résulte qu'à l'équilibre thermique (égalité entre l'absorption et l'émission de la chaîne), il est possible de définir une température typique pour la chaîne, liée directement à l'énergie de celle-ci et donc à son mode d'oscillation. Ce paramètre de température n'a malheureusement pas encore été programmé. C'est pourquoi on a supposé dans cette appliquette que la chaîne de ressorts était au repos initialement.

Chaîne libre initialement au repos : On placera bien évidemment chaque atome de carbone distant de d_ccde ses voisins. Il y a donc un seul arbitraire : la position du centre de gravité de la chaîne. Et celle-ci sera fixée à l'origine à l'instant t=0. Quant aux vitesses, elles sont toutes choisies nulles.
Chaîne liée au bâti initialement au repos : Cette fois-ci, c'est le bâti qui servira d'origine. Les atomes seront ensuite placés les uns après les autres en respectant la distance d'équilibre d_cc. Les vitesses sont encore une fois choisies toutes nulles.
Particule d'hydrogène : L'appliquette fixe initialement la particule d'hydrogène à 5r₀ du dernier atome de carbone. A cette distance le potentiel de Lennard Jones ne se fait pas encore trop sentir, de cette manière la particule incidente a un trajectoire quasi rectiligne avant de "plonger" et d'interagir avec le carbone. Pour résumer, cette distance initiale permet de voir la transition "trajectoire rectiligne uniforme - trajectoire accélérée". Quant à la vitesse v₀, l'utilisateur pourra la fixer selon son bon vouloir. Elles devront cependant être négatives pour pouvoir se rapprocher du dernier atome de carbone (sinon pas d'interaction !). Précisons de même que leur norme typique va de 500 à 1500 ms^-1.

Voici le vecteur (position, vitesse) initial pour la chaîne libre : $\left\{ (x_i(0),\dot{x_i}(0))=\left(\frac{d_{cc}}{2}[2i-(n+1)],0\right) \right\}_{i=1..n}$ .

Et pour l'hydrogène : $y(0)=y_0>d_{cc}(n-1)/2~~{\rm et}~~\dot{y}(0)=v_0<0$ .

Il reste une petite subtilité pour fournir ces données à l'intégrateur. En effet, les équations différentielles étant dans l'espace d'Hamilton, il faut passer les conditions initiales dans ce même espace. A savoir faire suivre la cheminement $(x_i,\dot{x_i})~~\rightarrow_{\rm Dedimensionnement}~~(X_i,\dot{X_i})~~\rightarrow_{\rm Lagrange}~~(q_i,\dot{q_i})~~\rightarrow_{\rm Hamilton}~~(q_i,p_i)$ aux conditions initiales précédemment explicitées. Ce qui donne, en passant les calculs :
${}^t \! X(0)=(\overbrace{\frac{x_1(0)}{r_0},...,\frac{x_n(0)}{r_0},\frac{y_0}{r_0}-\frac{x_n(0)}{r_0}}^{{}^t \! Q(0)},\overbrace{0,...,0,\mu \frac{v_0}{r_0} \sqrt{\frac{M}{k}},\mu \frac{v_0}{r_0} \sqrt{\frac{M}{k}} }^{{}^t \!P(0)})$

Ce qu'on fournit à l'intégrateur numérique est donc la fonction F et ce vecteur initial X(0). Il renverra le vecteur X(t).

Retour dans l'espace réel et redimensionnement

Il restera alors simplement à faire revenir ce vecteur X(t) qui sera dans l'espace d'Hamilton, dans l'espace réel. Il n'y aura qu'à suivre exactement le chemin inverse de tout à l'heure. Connaissant donc X c'est à dire les q_i et les p_i à chaque instant, on en déduira les positions et vitesses de chaque atome grâce aux formules :
$x_i=q_i r_0 {\rm~pour~} i=1..n \qquad y=r_0(q_{n+1}+q_n)$ pour les positions
$\dot{x_i}=r_0 \sqrt{\frac{k}{M}} p_i {\rm~pour~} i=1..n-1 \qquad \dot{x_n}=r_0 \sqrt{\frac{k}{M}}(p_n-p_{n+1}) \qquad \dot{y}=r_0 \sqrt{\frac{k}{M}} \frac{1}{\mu}p_{n+1}$ pour les vitesses.

Intégrateur numérique

L'intégrateur numérique est un outil qui sert à résoudre une équation différentielle donnée. Dans notre cas (étude de la dynamique d'un système à plusieurs corps) la mise en équations du problème aboutit à un système d'équations différentielles couplées. Cependant, si elle ne contient pas de termes en y', on peut toujours transformer une équation différentielle du second ordre en deux équations différentielles du premier ordre via l'introduction d'une variable intermédiaire : $y = a\cdot y '' \Leftrightarrow \{ z=y' \quad {\rm et} \quad y=a\cdot z'\}$ .

Il suffit donc que l'intégrateur numérique soit capable de résoudre un système d'équations différentielles du premier ordre pour pouvoir résoudre n'importe quel système d'équations différentielles linéaires.

On veut donc résoudre un système du type : $\dot{X} = F(X)$ où X est une grandeur vectorielle.

Méthode de Runge-Kutta

La méthode de Runge-Kutta peut être vue comme une amélioration de la méthode d'Euler. Dans la méthode d'Euler, on résout numériquement une équation différentielle de type y'(x)=f(y(x))

de la façon suivante : on connaît une condition initiale (x_0, y(x_0))

. On peut donc calculer f(y(x_0))

qui n'est autre que y'(x_0)

. En choisissant un pas d'intégration h on peut donc déterminer approximativement y(x_0+h)

:
$y(x_0+h) \approx y(x_0) + hy'(x_0)$

Connaissant y(x_0+h)

on peut ensuite évaluer y'(x_0+h)

, puis calculer y(x_0+2h)

grâce à la formule précédente et ainsi de suite...

Bien qu'assez simple, cette méthode accumule une erreur assez importante rapidement car l'erreur lors du calcul de y(x_0+h)

est d'ordre h^2

. On cherche donc à évaluer de façon plus précise y(x_0+h)

. La méthode de Runge-Kutta d'ordre 2 consiste à évaluer y(x_0+h/2)

par une méthode d'Euler et à utiliser ensuite cette valeur pour calculer la dérivée y'(x_0+h/2)

. On calcule ensuite y(x_0+h)

par :
$y(x_0+h) \approx y(x_0) + hf(y(x_0)+h/2)$

On peut prouver que l'erreur commise sur y(x_0+h)

dans ce cas est d'ordre h^3

, on peut par ailleurs se convaincre du fait que cette méthode est plus efficace en traçant un graphe (avec une fonction y(x)

convexe par exemple). Bien que l'erreur commise sur une itération soit d'ordre h^3

, cette méthode s'appelle RK2 car l'erreur commise sur toute la durée d'intégration est h^2

On voit donc l'utilité de la méthode de Runge-Kutta : on double le nombre de calculs nécessaires (puisqu'il faut évaluer deux fois la dérivée) mais on gagne un ordre sur la précision. Il est ensuite possible de faire plus d'évaluations intermédiaires de la dérivée afin d'obtenir des erreurs encore plus faibles. Par exemple la méthode RK4, où :
$y(x_0+h) \approx y(x_0) + \frac{h}{6}(k_1 + k_2 + k_3 + k_4)$
$k_1 = f(y(x_0)) \quad ; \quad k_2 = f(y(x_0)+\frac{h}{2}k_1)$
$k_3 = f(y(x_0)+\frac{h}{2}k_2) \quad ; \quad k_4 = f(y(x_0)+hk_3)$

permet d'obtenir une erreur en h^5

sur le pas, d'où la dénomination RK4. C'est cette méthode qui est utilisée dans le code.

Pas variable

On peut maintenant intégrer des systèmes d'équations différentielles avec une bonne précision mais un problème persiste : le pas d'intégration

doit être fixé à l'avance. Or, pour une trajectoire X(t)

donnée, on va accumuler une erreur importante quand les dérivées n-ièmes ont une grande valeur, et au contraire on va commettre peu d'erreur quand les dérivées n-ièmes sont petites. Pour chaque itération de l'intégration on va chercher le pas

engendrant une erreur fixée à l'avance par l'utilisateur.

Le principe est donc le suivant (on considère une fonction de la variable temporelle x(t)

solution de l'équation différentielle $\dot{x}(t) = f(x(t))$ : à partir d'une position initiale x(t_0)

on calcule x(t_0+h)

en faisant un pas de taille

(on note

) et on calcule x(t_0+h)

en faisant deux pas de taille

(on note

). L'erreur commise sur x_1(t_0+h)

est

où C est une constante et l'erreur commise sur x_2(t_0 + h)

est

. On a donc $|x_1(t_0 + h)-x_2(t_0 + h)| \approx \frac{15Ch^5}{16} \approx Ch^5 \equiv E$ .

Si maintenant l'utilisateur à fixé une erreur relative sur la précision $\epsilon_r$ , l'erreur absolue vaut $E_0 = | \epsilon_r x_2(t_0 + h)| = Ch_0^5$ . h_0

est donc le pas permettant d'obtenir l'erreur relative $\epsilon_r$ . On obtient donc au final :
$h_0 = 0.95h\left(\frac{E}{E_0}\right)^\frac{1}{5}$

Le 0.95 permet de se laisser une marge de sécurité pour être sûr d'avoir une précision suffisante, puisque |x_1(t_0 + h)-x_2(t_0 + h)|

n'est pas exactement égal à

Une fois

calculé on refait toute la procédure jusqu'à ce que l'erreur commise avec h_0

soit inférieure à l'erreur souhaitée. On peut ensuite calculer
$x_3(t_0 + h) = \frac{2^4x_2(t_0 + h)-x_1(t_0 + h)}{2^4-1}$

ce qui permet d'avoir une erreur sur le pas en h^6

, et donc globalement un RK5.

Méthodes JAVA

L'implémentation de cet intégrateur est contenue dans une classe JAVA "intégrateur". Cette classe contient une méthode rk4 qui calcule une itération à la méthode RK4, une méthode rk4pv qui calcule une itération avec le pas optimal pour l'erreur demandé, et une classe rk4dt qui effectue une intégration sur une période T en tronquant le pas variable de sorte à ce que l'intégration finisse exactement au temps T. Ceci est nécessaire pour pouvoir afficher les points des solutions à intervalles réguliers.

INTERaction GAz-Grain

Table des matières

Champ

Introduction

Objectifs

Explications

Modèles

La structure carbonée

Les interactions Hydrogène-Carbone

Les équations du mouvement

Energies cinétiques, potentielles et dédimensionnement

Lagrangien

Hamiltonien

Équations d'Hamilton

Résolution des équations d'Hamilton

Conditions initiales

Retour dans l'espace réel et redimensionnement

Intégrateur numérique

Méthode de Runge-Kutta

Pas variable

Mode d'emploi de l'applet

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Exemple d'utilisation

Liste des paramètres d'entrée de l'applet

Liste des paramètres de sortie de l'applet