# Méthode de Grad

À la fin du XIXe siècle on connaît l'équation de Boltzmann qui régit la dynamique du milieu gazeux à l'échelle microscopique et les équations d'Euler et de Navier-Stokes pour le niveau macroscopique. Le passage d'une échelle à l'autre constitue une partie du sixième problème de Hilbert. David Hilbert, auteur des énoncés des problèmes jugés majeurs à la fin du XIXe siècle pose les bases d'une méthode sous forme d'un développement qui porte son nom (1912). Il faudra attendre quelques années pour que Sydney Chapman et David Enskog proposent simultanément et indépendamment en 1916 et 1917 une solution à ce problème$^,$$^,$ par une méthode de perturbation consistant à définir la solution sous forme de série de fonctions de distribution en fonction d'un « petit paramètre » assimilable au nombre de Knudsen.
Harold Grad en 1949 a proposé une approche alternative consistant à chercher la solution par la méthode des moments de la fonction de distribution. L'équation de Boltzmann est multipliée par $m,m\,\mathbf {v} ,m\,\mathbf {v} \otimes \mathbf {v} ,m\,\mathbf {v} \otimes \mathbf {v} \otimes \mathbf {v} ...$ ($\mathbf {v}$ est la vitesse microscopique de l'équation de Boltzmann et $\otimes$ le produit tensoriel) et intégrée en vitesse. Dans ce type de méthode l'équation portant sur le ne moment fait apparaître le (n+1)e. Il faut donc faire une hypothèse pour « fermer » le système. Grad suppose la solution exprimée par une série tronquée de polynômes d'Hermite tensoriels. David Levermore a plus récemment (1996) proposé une fermeture qui fait appel à une propriété générale : la solution maximise l'entropie du système de fermions que sont les particules du milieu.
Ces formulations de la mécanique des fluides n'ont pas eu de développement notable car n'apportant pas d'avantage par rapport aux équations de Navier-Stokes.

## Équations d'évolution microscopique et macroscopique

On se limite à un milieu comportant une seule espèce.

### Niveau microscopique

On note $f(\mathbf {x} ,\mathbf {v} ,t)$ la fonction de distribution statistique de la vitesse $\mathbf {v}$ à l'instant $t$ au point $\mathbf {x}$ pour la particule (atome ou molécule) de masse $m$. Le nombre probable de particules dans le volume $[\mathbf {x} ,\mathbf {x} +\mathrm {d} \mathbf {x} ]$, de vitesses $[\mathbf {v} ,\mathbf {v} +\mathrm {d} \mathbf {v} ]$ à cet instant est  $f\,\mathrm {d} \mathbf {x} \,\mathrm {d} \mathbf {v}$. La distribution statistique $f$ se mesure donc en s$^3$ m$^{−6}$.
L'équation de Boltzmann s'écrit
       ${\frac {\partial f}{\partial t}}+\mathbf {v} \cdot \nabla f=\sum _{j}Q(f,f)$
où $Q$, l'opérateur (ou noyau) de collision, est un opérateur intégral quadratique décrit ci-dessous, donnant l'effet des collisions que l'on supposera élastiques pour simplifier le problème : pas d'échange entre degrés de liberté internes et translation. On exclut donc la viscosité volumique qui résulte de ce type d'échange et qui n'a d'intérêt que dans les problèmes de propagation du son.

#### La collision élastique

Les vitesses avant interaction sont $\mathbf {v} _{i}$ et $\mathbf {v} _{j}$ dans un référentiel galiléen. Ces vitesses valent $\mathbf {v} _{i}'$ et $\mathbf {v} _{j}'$ après interaction. On se place dans un système centré sur le barycentre qui a une vitesse constante du fait de la conservation de la quantité de mouvement. Dans ce système qui est donc galiléen la vitesse initiale de la particule $j$ est la vitesse relative  $\mathbf {g} _{ij}=\mathbf {v} _{i}-\mathbf {v} _{j}$. Par symétrie on peut affirmer que la trajectoire sera contenue dans le plan contenant l'origine et  $\mathbf {g} _{ij}$. On choisit un repère tel que  $\mathbf {\Omega } ={\frac {\mathbf {g} _{ij}}{||\mathbf {g} _{ij}||}}=[1,0,0]$ (voir figure). Dans ce repère la déviation est  $\theta _{ij}$, fonction du paramètre d'impact $b$, de la vitesse relative $g_{ij}$ et du potentiel d'interaction que l'on suppose ne dépendant que de la distance en les deux particules en interaction. Si cette hypothèse est rigoureuse pour l'interaction entre deux atomes, on peut la considérer utilisable pour deux molécules : le potentiel est alors un potentiel moyen statistique.
La direction de sortie d'interaction est définie par  $\mathbf {\Omega '} ={\frac {\mathbf {g'} _{ij}}{||\mathbf {g'} _{ij}||}}$. On peut calculer les vitesses finales à partir des considérations suivantes :
* la conservation de la quantité de mouvement dans l'interaction implique
       $\mathbf {v} _{i}'+\mathbf {v} _{j}'=\mathbf {v} _{i}+\mathbf {v} _{j}$
* la vitesse relative $\mathbf {g} _{ij}$ a un module constant du fait de la conservation de l'énergie, donc $g_{ij}'=g_{ij}$ ou
       $|\mathbf {v} _{i}'-\mathbf {v} _{j}'|=|\mathbf {v} _{i}-\mathbf {v} _{j}|$
Les vitesses après l'interaction sont donc
       $\mathbf {v} _{i}'(b,g)=\mathbf {v} _{i}-(\mathbf {g_{ij}} \cdot \mathbf {\Omega '} )\,\mathbf {\Omega '}$
       $\mathbf {v} _{j}'(b,g)=\mathbf {v} _{j}+(\mathbf {g_{ij}} \cdot \mathbf {\Omega '} )\,\mathbf {\Omega '}$
De plus la conservation du moment cinétique au cours de l'interaction conduit à $b'=b$.
Le système décrivant la collision est réversible. Le théorème de Liouville permet donc d'écrire
       $\mathrm {d} \mathbf {v} _{i}'\,\mathrm {d} \mathbf {v} _{j}'=\mathrm {d} \mathbf {v} _{i}\,\mathrm {d} \mathbf {v} _{j}$

#### Le noyau de collision

Le nombre probable de particules qui traversent l'aire  $2\,\pi \,b\,\mathrm {d} b$  par unité de temps est  $f(\mathbf {v} _{j})\,g_{ij}\,2\,\pi \,b\,\mathrm {d} b\,\mathrm {d} \mathbf {v} _{j}$. Elles interagissent avec le nombre probable de particules dans le volume élémentaire  $f(\mathbf {v} _{i})\,\mathrm {d} \mathbf {r} \,\mathrm {d} \mathbf {v} _{i}$. Le nombre de particules qui disparaissent de la statistique par unité de temps est  $\Theta _{ij}^{-}\,\mathrm {d} \mathbf {r} \,\mathrm {d} \mathbf {v} _{i}$ avec
       $\Theta _{ij}^{-}=2\pi \int _{\mathbf {v} }\int _{0}^{\infty }f(\mathbf {v} _{i})f(\mathbf {v} _{j})\,g_{ij}\,b\,\mathrm {d} b\,\mathrm {d} \mathbf {v} _{j}$
On compte de la même façon la quantité de particules qui apparaissent
       $\Theta _{ij}^{+}=2\pi \int _{\mathbf {v} }\int _{0}^{\infty }f(\mathbf {v} _{i}')f(\mathbf {v} _{j}')\,g_{ij}'\,b'\,\mathrm {d} b'\,\mathrm {d} \mathbf {v} _{j}'$
Compte tenu des relations données ci-dessus pour la collision l'opérateur de collision s'écrit
       $Q(f,f)=\Theta _{ij}^{+}-\Theta _{ij}^{-}=2\pi \int _{\mathbf {v} }\int _{0}^{\infty }[f(\mathbf {v} _{i}')f(\mathbf {v} _{j}')-f(\mathbf {v} _{i})f(\mathbf {v} _{j})]\,g_{ij}\,b\,\mathrm {d} b\,\mathrm {d} \mathbf {v} _{j}$
Cette équation est nommée équation de Wang Chang et Uhlenbeck.
On peut donner une formulation équivalente en introduisant la section efficace différentielle $\sigma _{ij}$ définie par
       $2\,\pi \,b\,\mathrm {d} b=2\,\pi \,\sigma _{ij}\sin \theta _{ij}\,\mathrm {d} \theta _{ij}=\sigma _{ij}\,\mathrm {d} \Omega$
d'où
       $Q(f,f)=\int _{\mathbf {v} }\int _{4\pi }[f(\mathbf {v} _{i}')f(\mathbf {v} _{j}')-f(\mathbf {v} _{i})f(\mathbf {v} _{j})]\,g_{ij}\,\sigma _{ij}\,\mathrm {d} \mathbf {\Omega }$

### Niveau macroscopique
#### Les variables

L'équation de Boltzmann décrit au niveau microscopique l'évolution de particules. Au niveau macroscopique on définit les quantités suivantes, fonctions de x et t
| - la densité particulaire | $n=\int _{\mathbf {v} }f\,\mathrm {d} \mathbf {v}$ |
| - la masse volumique | $\rho =n\,m$ |
| - la vitesse moyenne | $\mathbf {V} ={\frac {1}{n}}\int _{\mathbf {v} }\mathbf {v} f\mathrm {d} \mathbf {v}$ |
| - l'énergie interne | $e=\int _{\mathbf {v} }{\frac {1}{2}}mv^{2}f\mathrm {d} \mathbf {v}$ |

#### Les flux

Le flux de la quantité $\psi$ est par définition la quantité $\int _{\mathbf {v} }\psi f\,\mathrm {d} \mathbf {v}$, fonction de x et t
En notant $\otimes$ le produit dyadique on définit les moments d'ordre 2 (tenseur de pression) et 3 (tenseur de flux d'énergie) ainsi que les restrictions habituelles : pression scalaire et flux de chaleur vectoriel
| - tenseur de pression | $p_{ij}=\int _{\mathbf {v} }m\mathbf {v} \otimes \mathbf {v} f\,d\mathbf {v}$ | qui représente le flux de quantité de mouvement. |
| - flux du tenseur de pression | $p_{ijk}=\int _{\mathbf {v} }m\mathbf {v} \otimes \mathbf {v} \otimes \mathbf {v} f\mathrm {d} \mathbf {v}$ |
| - tenseur de flux d'énergie | $q_{ij}=\int _{\mathbf {v} }{\frac {1}{2}}mv^{2}\mathbf {v} \otimes \mathbf {v} f\mathrm {d} \mathbf {v}$ | flux du flux de chaleur |
En particulier on peut extraire les quantités habituelles
| - pression scalaire | $p={\frac {1}{3}}\int _{\mathbf {v} }mv^{2}fd\mathbf {v}$ | définie à partir de la trace du tenseur de pression. |
| - vecteur flux de chaleur | $q_{i}=\int _{\mathbf {v} }{\frac {1}{2}}mv^{2}\mathbf {v} f\mathrm {d} \mathbf {v}$ | contraction du tenseur de flux d'énergie. |
On définit la température à partir de l'équation d'état  $p=n\,k\,T$

## Équations des moments de la fonction de distribution

L'équation de Boltzmann est multipliée par $m$ puis par  $1,\mathbf {v} ,\mathbf {v} \otimes \mathbf {v} ,\mathbf {v} \otimes \mathbf {v} \otimes \mathbf {v} ...$  et intégrée en vitesse. On obtient le système suivant en se restreignant aux quantités  $\rho ,\mathbf {v} ,{\mathsf {P}},\mathbf {q}$, soit 1 + 3 + 6 + 3 = 13 quantités indépendantes (« méthode des 13 moments »)
       ${\frac {\partial \rho }{\partial t}}+{\frac {\partial (\rho v_{i})}{\partial x_{i}}}=0$
       ${\frac {\partial \rho v_{i}}{\partial t}}+{\frac {\partial (\rho v_{i}v_{j}+p_{ij})}{\partial x_{j}}}=0$
       ${\frac {\partial p_{ij}}{\partial t}}+{\frac {\partial (p_{ijk}+p_{ij}v_{k})}{\partial x_{k}}}+p_{ki}{\frac {\partial v_{j}}{\partial x_{k}}}+p_{kj}{\frac {\partial v_{i}}{\partial x_{k}}}=P_{ij}$
       ${\frac {\partial q_{i}}{\partial t}}+{\frac {\partial (q_{ij}+q_{i}v_{j})}{\partial x_{j}}}+p_{ijk}{\frac {\partial v_{j}}{\partial x_{k}}}+q_{j}{\frac {\partial v_{i}}{\partial x_{j}}}-{\frac {p_{ki}}{\rho }}{\frac {\partial p_{kj}}{\partial x_{j}}}-{\frac {1}{2}}{\frac {p_{rr}}{\rho }}{\frac {\partial p_{ij}}{\partial x_{j}}}=Q_{i}$
où les quantités  $P_{ij}$  et  $Q_{i}$  correspondent au noyau de collision
       $P_{ij}=2\pi \int _{v_{i}}\int _{v_{j}}\int _{0}^{\infty }m\,(v'_{i}v'_{j}-v_{i}v_{j})f(v_{i})f(v_{j})\,g\,b\,\mathrm {d} b\,\mathrm {d} v_{i}\mathrm {d} v_{j}$
       $Q_{i}=2\pi \int _{v_{i}}\int _{v_{j}}\int _{0}^{\infty }{\frac {1}{2}}m\,(v'^{2}v'_{j}-v^{2}v_{j})f(v_{i})f(v_{j})\,g\,b\,\mathrm {d} b\,\mathrm {d} v_{i}\mathrm {d} v_{j}$
Comme dans tout système aux moments, celui-ci n'est pas fermé : il n'existe pas d'équation pour p$_{ijk}$. Il faut donc « fermer » le système par une hypothèse ad hoc.

## Méthode de Grad

Dans cette méthode on cherche la fonction de distribution sous forme d'un développement en série de polynômes d'Hermite tensoriels en négligeant une partie des termes du second ordre  $a_{ijk}H_{ijk}$ .
       $f=f^{(0)}\left(aH+a_{i}H_{i}+{\frac {1}{2}}a_{ij}H_{ij}+{\frac {1}{10}}a_{rri}H_{ssi}\right)$
où f$^{(0)}$ est la distribution maxwellienne
       $f^{(0)}=n\left({\frac {m}{kT}}\right)^{\frac {3}{2}}\omega (\xi )\,,\qquad \omega (\xi )={\frac {e^{-{\frac {\xi ^{2}}{2}}}}{(2\pi )^{\frac {3}{2}}}}$
cette expression étant applicable à toutes les composantes
       $\xi _{i}={\sqrt {\frac {m}{kT}}}v_{i}$
Les polynômes d'Hermite sont définis par
       $H_{i_{1}i_{2}...i_{N}}={\frac {(-1)^{N}}{\omega }}{\frac {\partial ^{N}\omega }{\partial \xi _{1}\partial \xi _{2}...\partial \xi _{N}}}$
Ils sont orthogonaux avec le poids ω
       $\int \omega H_{i_{1}i_{2}...i_{N}}H_{j_{1}j_{2}...j_{M}}=\delta _{NM}(\delta _{i_{1}j_{1}}\delta _{i_{2}j_{2}}...\delta _{i_{N}j_{N}}+...)$
le terme au second membre comportant toutes les permutations possibles sur j.
Les premiers polynômes sont
       ${\begin{array}{lcl}H(\xi )&=&1\\[0.6em]H_{i}(\xi )&=&\xi _{i}\\[0.6em]H_{ij}(\xi )&=&\xi _{i}\xi _{j}-\delta _{ij}\\[0.6em]H_{ijk}(\xi )&=&\xi _{i}\xi _{j}\xi _{k}-(\xi _{i}\delta _{jk}+\xi _{j}\delta _{ik}+\xi _{k}\delta _{ij})\end{array}}$
En reportant le développement dans les équations du système on obtient
       $a_{1}=1\,,\qquad a_{i}=0\,,\qquad a_{rr}=0\,,\qquad a_{ij}={\frac {p_{ij}}{p}}\,,\qquad a_{rri}={\frac {2q_{i}}{p}}{\sqrt {\frac {m}{kT}}}$
On définit  $p'_{ij}=p_{ij}-{\frac {1}{3}}p_{rr}\delta _{ij}$  qui est le déviateur du tenseur de pression (tenseur de cisaillement).
La fonction de distribution de Grad est
       $f=f^{(0)}\left\{1+{\frac {2\beta ^{2}}{\rho }}\left[p'_{ij}v_{i}v_{j}+{\frac {4}{5}}q_{i}v_{i}\left(\beta v^{2}-{\frac {5}{2}}\right)\right]\right\}\,,\qquad \beta ={\frac {m}{2kT}}$
On peut alors calculer les moments d'ordre élevé (fermeture du système) à partir des relations de définition
       ${\begin{array}{rcl}p_{ijk}&=&{\frac {2}{5}}\left(q_{i}\delta _{jk}+q_{j}\delta _{ik}+q_{k}\delta _{ij}\right)\\[0.6em]q_{ij}&=&{\frac {5p^{2}}{2\rho }}\delta _{ij}+{\frac {7p}{2\rho }}p'_{ij}\end{array}}$
À partir de cette solution on peut calculer les coefficients de transport : viscosité dynamique et conductivité à partir des intégrales de collision. Ces coefficients sont les mêmes que ceux issus de la méthode de Chapman-Enskog.