Théorie des figures d'équilibre

On peut simplifier le problème encore davantage en supposant que la figure d'équilibre constitue un rotateur permanent. On définit un corps en rotation permanente comme suit[4] :

• • il est isolé dans l'espace et tourne autour d'un axe fixe ;
  • son mouvement est stationnaire dans un système d'inertie, et la densité de chaque point massique individuel reste constante au cours du mouvement ;
  • on peut négliger les effets dissipatifs, notamment les frictions ;
  • aucune force électromagnétique n'agit sur le corps.

Constatons que l'on n'a pas encore spécifié la vitesse angulaire ω. Admettons maintenant que le corps tourne autour de l'axe (O x₃) et prenons l'origine O du système d'axes au centre de masse du corps. Considérons des coordonnées cylindriques (h, λ, z) définies en termes des coordonnées cartésiennes (x₁, x₂, x₃) par la loi de transformation :

${\displaystyle {\begin{cases}h=(x_{1}^{2}+x_{2}^{2})^{\frac {1}{2}}\\\lambda =\mathrm {arctan} ({\frac {x_{2}}{x_{1}}})\\z=x_{3}\end{cases}}}$.

Ainsi, dans un système d'inertie, les composantes de la vitesse prennent la forme

${\displaystyle {\begin{cases}v_{h}=0\\v_{\lambda }=\omega h\\v_{z}=0\end{cases}}}$.

Comme la densité d'un point massique appartenant à un rotateur permanent ne varie pas le long de sa trajectoire, nous avons

${\displaystyle {\frac {\mathrm {d} \rho }{\mathrm {d} t}}={\frac {\partial \rho }{\partial t}}+v_{j}{\frac {\partial \rho }{\partial x_{j}}}=0}$.

L'équation de continuité implique donc que le mouvement est indivergentiel[5], c'est-à-dire

${\displaystyle {\frac {\partial v_{j}}{\partial x_{j}}}=0}$.

En coordonnées cylindriques, cette équation devient[6]

${\displaystyle {\frac {\partial \omega }{\partial \lambda }}=0}$,

ce qui implique qu'une rotation permanente est caractérisée par un champ de vitesse qui ne dépend pas de la variable λ[7], soit

${\displaystyle \omega =\omega (h,z)}$ .

Un cas particulier d'un rotateur permanent est un rotateur rigide, pour lequel la vitesse angulaire est constante :

${\displaystyle \omega =\mathrm {Constant} }$.

Le corps concerné tourne alors en bloc sur lui-même : on dit que la figure d'équilibre correspondante est en équilibre relatif.

Or, l'hypothèse d'une rotation permanente ne requiert pas seulement que dρ/dt = 0, mais aussi que ∂ρ/∂t = 0, d'où l'on tire au moyen de l'équation de continuité que

${\displaystyle v_{j}{\frac {\partial \rho }{\partial x_{j}}}=0}$.

Cette dernière relation montre que pour une rotation permanente les lignes de courant sont orthogonales au gradient de densité local. Dans le système de coordonnées cylindriques que nous avons adopté, cette relation s'écrit

${\displaystyle \omega {\frac {\partial \rho }{\partial \lambda }}=0}$.

En supposant que le corps tourne effectivement (ω ≠ 0) comme un rotateur permanent par rapport à un référentiel d'inertie, nous trouvons donc que

${\displaystyle \rho =\rho (h,z)}$.

On en déduit que la distribution de densité dans un rotateur permanent possède une symétrie axiale, et la même conclusion s'impose (voir ci-dessus) pour le champ de vitesse angulaire. Ainsi, ce champ ω = ω(h, z) ne possède pas de tourbillon[8] dans le plan (O x₁ x₂), donc[9]

${\displaystyle x_{2}{\frac {\partial \omega }{\partial x_{1}}}-x_{1}{\frac {\partial \omega }{\partial x_{2}}}=0}$.

En termes de coordonnées cartésiennes, le champ de vitesse peut donc s'écrire

${\displaystyle {\begin{cases}v_{1}=-\omega x_{2}\\v_{2}=\omega x_{2}\\v_{3}=0\end{cases}}}$.

En vertu des relations ∂v_j/∂x_j = 0 et v_j ∂ρ/∂x_j = 0, la quantité ∂(ρ v_i v_j) / ∂_j devient simplement ρ v_j ∂v_i / ∂x_j. Dès lors, l'équation de mouvement prend la forme particulière

${\displaystyle -\rho \omega ^{2}(x_{i}-\delta _{i3}x_{3})=f_{i}-{\frac {\partial p}{\partial x_{i}}}}$.

Selon la définition d'un rotateur permanent énoncée plus haut, il est évident que le champ de force volumique f_i est simplement le champ de force engendré par l'autoattraction gravitationnelle, c'est-à-dire

${\displaystyle f_{i}=\rho g_{i}}$.

Il est bien connu que la force gravifique spécifique g_i en un point P est fournie par l'expression

${\displaystyle g_{i}(P)=G\int _{B}{\frac {\rho (Q)}{d^{2}(P,Q)}}e_{i}(P,Q)\;\mathrm {d} \tau (Q)}$,

où d(P, Q) représente la distance entre le point attiré P et un quelconque point attirant Q du corps B. La quantité e_i désigne la i^e composante d'un vecteur unitaire porté par le segment de droite qui joint P à Q, et orienté de P à Q. L'élément de volume en Q est dénoté par dτ(Q). On sait de même que la force volumique g_i dérive en chaque point de l'espace d'une fonction scalaire V, appelé potentiel gravitationnel. Ainsi,

${\displaystyle g_{i}=-{\frac {\partial V}{\partial x_{i}}}}$.

Nous définissons ici ce potentiel V, en accord avec la convention généralement adoptée par les physiciens, comme l'énergie potentielle par unité de masse d'un élément matériel situé au point P. Cette énergie potentielle résulte de l'interaction gravitationnelle de l'élément matériel en P avec tous les autres éléments matériels composant le corps B, de sorte que

${\displaystyle V(P)=-G\int _{B}{\frac {\rho (Q)}{d(P,Q)}}\;\mathrm {d} \tau (Q)}$.

Selon cette définition, le potentiel gravifique est toujours négatif, et correspond physiquement à une énergie de liaison. Il faut remarquer que les astronomes et les géodésiens, et aussi de nombreux géophysiciens, adoptent une définition différente pour le potentiel gravifique, à savoir ils l'assimilent au travail qu'un élément matériel de masse unitaire accomplit dans le champ de force gravifique g_i en passant d'un point situé à l'infini au point P. Ce travail représente l'énergie potentielle par unité de masse changée de signe, et est donc positif. La convention des physiciens adoptée ici possède l'avantage, en retenant la correspondance directe entre potentiel et énergie potentielle, de faire en sorte qu'un point massique tend toujours à mouvoir vers des régions où le potentiel est plus bas. La constante G ci-dessus est la constante de gravitation.

En fait, la symétrie axiale de V aurait pu être déduite de la symétrie axiale de ρ en se servant du principe de symétrie établi par le physicien français Pierre Curie (1859-1906). Ce principe peut s'énoncer de façon concise comme suit :

Tout effet produit par une cause déterminée doit posséder au moins les mêmes éléments de symétrie mais peut être plus symétrique que la cause.

Ou encore : toute dissymétrie rencontrée dans un effet donné doit aussi exister dans la cause qui a produit l'effet.

En corollaire, s'il existe une relation de réciprocité entre une cause et un effet, dans le sens que la cause peut tout aussi bien être considérée comme l'effet, et l'effet comme la cause, les deux — cause et effet — doivent posséder exactement les mêmes éléments de symétrie[10]. Un exemple de ce cas est fourni par une relation entre les tensions et les déformations, à savoir ce que nous avons appelé plus haut une « relation constitutive ». Il semble que l'on n'ait pas encore trouvé jusqu'à présent d'expérience ni d'observation se rapportant à des phénomènes macrophysiques qui contredise le principe de Curie, et ce dernier devrait être considéré comme une loi universelle de la Nature au même titre que les lois de la thermodynamique, par exemple[11]. Nous allons adopter ce point de vue ici pour déduire certaines propriétés des figures d'équilibre.

Tout d'abord, nous notons que la force d'attraction gravifique est centrale, et possède donc le plus haut degré de symétrie, celui de la sphère. Le principe de Curie requiert alors que l'arrangement des masses produit par la seule action de l'attraction gravifique doit correspondre à la symétrie sphérique : une quelconque configuration moins symétrique est exclue si l'on veut avoir équilibre. Une preuve mathématique du fait que la sphère est l'unique figure d'équilibre possible d'un corps fluide isolé au repos fut donnée par Lyapunov[12] en 1884.

Ensuite, considérons l'action combinée d'une force d'attraction gravifique et d'une force axifuge. Cette dernière introduit une certaine asymétrie par rapport à la sphère par le fait que la force axifuge possède la symétrie d'un cylindre de révolution. Le principe de Curie requiert donc que la distribution d'équilibre des masses possède la même symétrie axiale, puisqu'il existe une relation de réciprocité entre la distribution de densité et le potentiel. C'est le cas étudié plus haut en termes mathématiques.

Un point particulier doit retenir toute notre attention à ce stade. En effet, afin de compléter la formulation du problème, nous devons inclure une équation d'état appropriée. Sous forme générale, celle-ci peut s'écrire

${\displaystyle p=p(\rho ,T,\varphi _{1},...\varphi _{N})}$.

La lettre p désigne la pression, ρ est la densité, T est la température, et l'ensemble des variables φ₁ à φ_N décrit la composition chimique et minéralogique. Par définition, une structure baroclinique (ou, simplement, un barocline) est un système pour lequel une relation physique de cette forme existe. Toutefois, dans certaines circonstances spéciales qui sont approximativement remplies dans le cas de l'intérieur de la Terre[13], la pression peut simplement être reliée à la densité, à savoir

${\displaystyle p=p(\rho )}$.

Si tel est le cas, on dira qu'il s'agit d'une équation d'état barotropique et la structure correspondante est un barotrope. En première approximation, nous pouvons considérer la Terre comme un barotrope.

La distinction fondamentale entre modèles barotropiques et modèles barocliniques se trouve dans leurs stratifications respectives. Il est clair que les surfaces d'égale densité (surfaces isopycniques) et les surfaces d'égale pression (surfaces isobares) coïncident dans un barotrope. Par contre, dans les baroclines les surfaces d'égale densité coupent les surfaces d'égale pression, à moins que des conditions très spécifiques ne soient remplies.

Tenant compte de la symétrie axiale de ρ et V, l'équation de mouvement d'un rotateur permanent devient en coordonnées cylindriques

${\displaystyle {\frac {1}{\rho }}{\frac {\partial \rho }{\partial h}}=-{\frac {\partial V}{\partial h}}+\omega ^{2}h}$,

${\displaystyle {\frac {1}{\rho }}{\frac {\partial p}{\partial z}}=-{\frac {\partial V}{\partial z}}}$,

${\displaystyle {\frac {1}{\rho h}}{\frac {\partial p}{\partial \lambda }}=0}$.

Le fait que ∂p/∂λ = 0 établit la symétrie de révolution du champ de pression, c'est-à-dire

${\displaystyle p=p(h,z)}$.

On définit les composantes de la pesanteur en coordonnées cylindriques par les relations

${\displaystyle g_{h}=-{\frac {\partial V}{\partial h}}+h\omega ^{2}}$,

${\displaystyle g_{z}=-{\frac {\partial V}{\partial z}}}$,

${\displaystyle g_{\lambda }=0}$.

Dès lors, dans un corps en rotation permanente, la pesanteur est partout orthogonale aux isobares, comme on le constate en se souvenant de la définition du gradient. Constatons que cette propriété est très générale : elle est valable tout aussi bien pour les baroclines que pour les barotropes, et qu'elle est valable pour une rotation uniforme comme pour une rotation différentielle, pour autant que le champ de vitesse angulaire conserve une symétrie axiale. En général, on définit le potentiel axifuge comme

${\displaystyle Z=-{\frac {1}{2}}\omega ^{2}h^{2}}$,

mais il est évident que cette définition n'est valable que pour une rotation uniforme. Si tel est le cas, nous pouvons employer un système d'axes fixé dans le corps et tournant avec lui. Si chaque particule du corps est au repos par rapport à ce référentiel qui n'est pas d'inertie, on dit que le corps est en équilibre hydrostatique relatif ou, plus simplement, en équilibre hydrostatique. Le plus souvent, dans les théories géodésiques ou géophysiques, on admet que tel est le cas pour la Terre.

Au contraire, si la vitesse angulaire n'est pas constante d'un point à l'autre, il ne peut exister qu'un équilibre hydrodynamique. Une rotation différentielle de ce type joue probablement un rôle important dans les étoiles et les planètes géantes, mais semble exclue dans les planètes telluriques à cause d'un trop fort couplage visqueux ou électromagnétique entre les différentes strates. Les études des figures d'équilibre générales donnent souvent le plus d'importance à la rotation zonale, que l'on nomme circulation méridienne ou courants de von Zeipel, ou courants d'Eddington, ou encore courants de Jardetzky[14]. Toutefois, l'étude des courants de convection ne fait généralement pas partie de la théorie des figures d'équilibre.

Les relations ci-dessus définissant les composantes de la pesanteur montrent qu'un potentiel axifuge Z, et donc une surface de niveau

${\displaystyle U=V+Z={\rm {Constante}}}$

peut exister si et seulement si ω ne dépend pas de z, c'est-à-dire si la vitesse angulaire est constante sur des surfaces cylindriques centrées sur l'axe de rotation, autrement dit si ω = ω(h). Dans ce cas plus général que la rotation uniforme, on définit le potentiel axifuge par

${\displaystyle Z=-\int _{0}^{h}\omega ^{2}(h')h'\;{\rm {d}}h'}$.

Nous constatons, en particulier, que pour une rotation uniforme, cette relation se réduit bien à la relation précédente comme il se doit.

Une conclusion importante qu'on peut tirer de l'existence d'un potentiel de rotation, et donc de l'existence d'un potentiel de pesanteur, est

${\displaystyle U(h,z)=V(h,z)-\int _{0}^{h}\omega ^{2}(h')h'\;{\rm {d}}h'}$,

de sorte que

${\displaystyle g_{h}=-{\frac {\partial U}{\partial h}}}$,

${\displaystyle g_{z}=-{\frac {\partial U}{\partial z}}}$,

${\displaystyle g_{\lambda }=0}$.

On donne à cette conclusion souvent le nom de théorème de Wavre-Poincaré[15], et on peut l'énoncer sous la forme suivante : Pour un rotateur permanent, l'une quelconque des propriétés qui suivent entraînent automatiquement les trois autres, à savoir

1. la vitesse angulaire est constante sur des cylindres centrés sur l'axe de rotation ;
2. la pesanteur dérive d'un potentiel ;
3. la pesanteur est normal aux surfaces d'égale densité ;
4. les surfaces d'égale densité coïncident avec les surfaces d'égale pression[16]. En effet, selon les équations pour le gradient de la pression et pour les composantes de la pesanteur fournies plus haut, on peut écrire pour la différentielle totale

${\displaystyle {\rm {d}}p={\frac {\partial p}{\partial h}}{\rm {d}}h+{\frac {\partial p}{\partial z}}{\rm {d}}z+{\frac {\partial p}{\partial \lambda }}{\rm {d}}\lambda }$

l'équation suivante :

${\displaystyle {\rm {d}}p=\rho \left(g_{h}{\rm {d}}h+g_{z}{\rm {d}}z\right)}$,

donc

${\displaystyle {\rm {d}}U=-{\frac {{\rm {d}}p}{\rho }}}$.

Un déplacement arbitraire sur une surface de niveau fournit dU = 0, impliquant donc que dp = 0. Les surfaces isobares coïncident donc avec les surfaces de niveau, et p = p(U) ou U = U(p). Nous trouvons ainsi à cause de la relation sur dU que la densité est aussi constante sur surface de niveau. En conclusion, les isobares, les isopycniques et les surfaces de niveau coïncident toutes s'il existe un potentiel de pesanteur, et la force par unité de masse due à la pesanteur est normale à ces surfaces. En particulier, elle est normale aux isopycniques. Réciproquement, considérons un corps pour lequel les isobares et les isopycniques coïncident. On peut alors définir une fonction U(p), en termes de laquelle la relation donnant la différentielle de la pression devient

${\displaystyle {\rm {d}}U=-g_{h}{\rm {d}}h-g_{z}{\rm {d}}z}$,

montrant que dU est une différentielle totale exacte. Par conséquent, les relations g_h(h,z) = – ∂U/∂h, g_z(h,z) = – ∂U/∂z, g_λ(h,z) = 0 doivent s'appliquer, et le champ de pesanteur dérive d'un potentiel. Finalement, supposons que la pesanteur est partout normale aux surfaces d'égale densité, soit gx∇ρ = 0. Comme g = – ∇U = ρ⁻¹ ∇p, la coïncidence des surfaces isobares et isopycniques est établie, et tous les énoncés du théorème de Wavre-Poincaré sont démontrés.

Il est important de se rendre compte que la condition ω = ω(h), ou la condition équivalente ∂ω/∂z = 0, ne requiert aucune connaissance préalable de l'équation d'état. Selon le théorème de Wavre-Poincaré, les isopycniques d'un barocline coïncideront donc avec ses isobares si ∂ω/∂z = 0. Toutefois, il est évident que de tels baroclines particuliers devraient être appelés de façon plus appropriée des pseudo-barotropes, parce qu'ils partagent la plupart de leurs propriétés avec les barotropes véritables, et se distinguent clairement des baroclines véritables.

Cette constatation nous permet de revenir aux considérations de symétrie. Tout d'abord, le principe de Curie nous apprend que les barotropes sont au moins aussi symétriques que les baroclines véritables, mais peuvent être plus symétriques.

Le principe de Curie (sans symétrie brisée) permet d'inférer une autre propriété de symétrie très importante des figures d'équilibre : un quelconque rotateur permanent barotropique ou pseudo-barotropique doit posséder un élément de symétrie supplémentaire, à savoir un plan équatorial, c'est-à-dire un plan normal à l'axe de rotation qui passe par le centre de masse. Cet énoncé constitue le théorème de Lichtenstein.

Mais cette preuve physique du théorème de Lichtenstein dans toute sa généralité exigerait l'unicité de solution, et la démonstration mathématique rigoureuse n'est pas triviale du tout[17].

Le cas d'un barocline véritable est quelque peu plus compliqué, mais Dive a démontré en 1930 qu'un barocline en rotation différentielle est symétrique par rapport à un plan équatorial si la vitesse angulaire est partout une fonction biunivoque de la densité et de la distance à l'axe de rotation[18].

Pour un barotrope quelconque en rotation lente, et en particulier pour les modèles de Terre, les surfaces de niveau ne diffèrent que légèrement de la forme sphérique. Il est donc naturel, compte tenu des éléments de symétrie que doivent posséder ces surfaces, de rechercher l'équation décrivant la surface d'une couche quelconque sous la forme

${\displaystyle r(s,\theta )=s\left[1+\sum _{n=0}^{\infty }s_{2n}(s)P_{2n}(cos\theta )\right]}$.

La variable s désigne ici un paramètre judicieusement choisi devant permettre de reconnaître de quelle couche il s'agit — pour simplifier, on dira que s représente le nom de la couche considérée. Ainsi, on peut choisir s comme rayon équatorial a, ou comme rayon polaire c, de la couche en question. Dans la suite, nous choisissons s comme rayon équivolumétrique moyen de la couche, autrement dit, comme rayon de la sphère qui possède le même volume que celui contenu sous la couche.

De la sorte, le problème de trouver la forme d'équilibre est ramené au problème de déterminer un ensemble infini de fonctions s_2n qu'on appelle les « fonctions de figure ». Afin de pouvoir utiliser la formule précédente en pratique, il est nécessaire que l'ordre de grandeur des fonctions de figure décroisse suffisamment vite avec l'ordre n pour assurer une convergence raisonnablement rapide de la série figurant dans le membre de droite. Pour les modèles de Terre en particulier, tel est bien le cas comme il sera montré plus loin.

Si nous tronquons la série à n = 1, nous sommes en présence de l'approximation du premier ordre qui correspond au cas classique de la théorie de Clairaut. Si nous tronquons la série à n = 2, nous obtenons l'approximation du second ordre qui correspond à la théorie de Darwin–de Sitter. Pour la Terre, l'approximation d'ordre 2 correspond à négliger des effets d'environ 10^-5 ; comme le rayon moyen de la Terre est de 6 371 km, on commet donc des erreurs qui peuvent se chiffrer typiquement à quelques mètres. De telles erreurs sont incompatibles avec la précision des mesures géodésiques atteinte à l'heure actuelle, qui estiment le rayon de la Terre à quelques centimètres près. Depuis la fin des années 1970, divers auteurs[20] ont considéré pour le cas de la figure hydrostatique de la Terre des approximations d'ordre 3 et en ont fourni les formules appropriées. En réalité, pour autant qu'on ne désire pas faire des calculs effectifs, il n'est guère plus difficile de développer la théorie des figures d'équilibre à un ordre général n = N que de se limiter à l'ordre n = 2 ou n = 3, par exemple.

Si nous négligeons les termes en ${\displaystyle {\mathcal {O}}(m^{2})}$ ou plus petits, les équations de figure se réduisent à la seule équation

${\displaystyle F_{2}=-s_{2}S_{0}+S_{2}+T_{2}-{\frac {1}{3}}m=0}$,

et l'aplatissement devient

${\displaystyle f=-{\frac {3}{2}}s_{2}}$.

Combinant les deux, on obtient l'équation intégro-différentielle de Clairaut

${\displaystyle fS_{0}+{\frac {3}{2}}(S_{2}+T_{2})-{\frac {1}{2}}m=0}$.

Les quantités S₂ et T₂ ne doivent être évaluées qu'à l'ordre ${\displaystyle {\mathcal {O}}(m)}$, et S₀ à l'ordre ${\displaystyle {\mathcal {O}}(m)}$. Dès lors, avec ${\displaystyle x=s/s_{1}}$ et ${\displaystyle \delta =\rho /{\bar {\rho }}}$, on trouve

${\displaystyle S_{2}(x)=-{\frac {2}{5}}{\frac {1}{x^{5}}}\int _{y=0}^{x}\delta (y){\rm {d}}(y^{5}f)}$,

${\displaystyle T_{2}(x)=-{\frac {2}{5}}\int _{y=x}^{1}\delta (y){\rm {d}}f}$.

Sous cette forme, l'équation de Clairaut est une équation intégro-différentielle linéaire pour l'aplatissement des strates internes de la Terre. Elle peut être résolue itérativement par la méthode générale indiquée plus haut pour les fonctions de figure. Toutefois, l'approche classique consiste à dériver l'équation de Clairaut sous la forme intégro-différentielle par rapport à x et éliminer les quantités m et T₂ en faisant usage des deux dernières relations, obtenant ainsi

${\displaystyle {\frac {5}{2}}S_{2}=\left(f-{\frac {1}{3}}x{\frac {{\rm {d}}f}{{\rm {d}}x}}\right)S_{0}}$,

puis à multiplier ce résultat intermédiaire par x⁵ et dériver encore une fois par rapport à x. On obtient alors l'équation différentielle de Clairaut

${\displaystyle {\frac {{\rm {d}}^{2}f}{{\rm {d}}x^{2}}}+6{\frac {\delta }{x\;S_{0}}}{\frac {{\rm {d}}f}{{\rm {d}}x}}+{\frac {6}{x^{2}}}\left({\frac {\delta }{S_{0}}}-1\right)f=0}$.

On peut intégrer numériquement cette dernière équation différentielle ordinaire en imposant que f reste fini en x = 0 (au centre) et que f satisfasse en x = 1 (à la surface) à la condition

${\displaystyle \left.{\frac {{\rm {d}}f}{{\rm {d}}x}}\right|_{1}+2f(1)={\frac {5}{2}}m}$.

On aboutit à cette condition en multipliant l'équation intégro-différentielle de Clairaut ci-dessus par x⁵, en posant ensuite x = 1 et en notant que S₀(1) = 1.

Tout au long du XIX^e siècle, en l'absence des données fournies plus tard par la sismologie, les astronomes essayaient de déterminer la densité des couches internes de la Terre en imposant deux contraintes à diverses fonctions d'essai ρ = ρ(x) : (1) la fonction d'essai, par intégration sur le volume, devait fournir la valeur observée M de la masse de la Terre ; (2) en intégrant l'équation de Clairaut, la fonction d'essai devait fournir la valeur observée f de l'aplatissement géométrique en surface. Or, avant l'avènement des ordinateurs et, surtout, des microordinateurs personnels dans la deuxième moitié du XX^e siècle, la solution de l'équation différentielle de Clairaut s'obtenait par une intégration numérique qui devait se faire à la main et qui était fort fastidieuse.

Pour cette raison, l'élégante méthode de solution proposée en 1885[21] par l'astronome français d'origine prussienne Rodolphe Radau (1835–1911), qui simplifie énormément l'intégration de l'équation de Clairaut, devint la méthode de résolution standard et reste jusqu'à présent d'une grande valeur dans les travaux théoriques ayant recours à l'aplatissement calculé au premier ordre[22].

La méthode de Radau consiste à transformer l'équation de Clairaut en introduisant la fonction sans dimension dite « paramètre (ou variable) de Radau »

${\displaystyle \eta (x)={\frac {x}{f}}{\frac {{\rm {d}}f}{{\rm {d}}x}}}$,

au moyen de laquelle on obtient l'équation de Clairaut-Radau

${\displaystyle {\frac {{\rm {d}}x^{5}{\sqrt {1+\eta }}\;S_{0}}{{\rm {d}}x}}=5x^{4}S_{0}\;F(\eta )}$,

avec

${\displaystyle F(\eta )={\frac {1+{\frac {\eta }{2}}-{\frac {\eta ^{2}}{10}}}{\sqrt {1+\eta }}}}$.

• R. Aris (1962). Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Prentice-Hall, Englewood Cliffs (Corrected reprint, Dover Publications, 1962). (ISBN 0-486-66110-5).
• C. Denis (1989). The Hydrostatic Figure of the Earth, dans Physics and Evolution of the Earth's Interior (volume 4), édité par R. Teisseyre : Gravity and Low-Frequency Geodynamics, PWN-Polish Scientific Publishers, Warszawa & Elsevier, Amsterdam, pp. 111–186. (ISBN 83-01-07986-X).
• R. Dugas (1950). Histoire de la Mécanique, Éditions du Griffon, Neuchatel & Éditions Dunod, Paris.
• W.S. Jardetzky (1958). Theories of Figures of Celestial Bodies, Interscience Publishers, New York.
• P. Lanzano (1982). Deformations of an Elastic Earth, Academic Press, New York. (ISBN 0-12-436620-1).
• W.A. Magnizki, W.W. Browar & B.P. Schimbirew (1964). Theorie der Figur der Erde, VEB Verlag für Bauwesen, Berlin.
• H. Moritz (1990). The Figure of the Earth (Theoretical Geodesy and the Earth's Interior), Wichmann Verlag, Karlsruhe. (ISBN 3-87907-220-5).
• I. Todhunter (1873). A History of the Mathematical Theories of Attraction and the Figure of the Earth (From the time of Newton to that of Laplace), 2 volumes, Macmillan and Company, London. [Fac-similé des deux volumes réunis en un seul publié par Dover Publications, New York, 1962.]
• C. Truesdell (1968). Essays on the History of Mechanics, Springer-Verlag, Heidelberg.
• V.N. Zharkov & V.P. Trubitsyn (1978). Physics of Planetary Interiors (Translated, edited, and with additional material by W.H. Hubbard), Pachart Publishing House, Tucson. (ISBN 0-912918-15-2).

• Concepts de base en théorie des milieux continus
• Équations fondamentales de la mécanique des milieux continus
• Figure de la Terre
• Champ de pesanteur

• Portail de la géodésie et de la géophysique
• Portail de l’astronomie