Valeur potentiel gravitationnel

[fabio123]

Bonsoir,

dans le cadre d'un projet, j'ai été amené à résoudre sous Matlab l'équation de poisson en coordonnées cylindriques.

Ci-dessous 2 figures représentant la solution finale de cette équation dans le cas d'une galaxie (en valeur absolue) :

9o2YqLL.jpg

ZmwRO69.jpg

Pour ceux qui ont Matlab, je mets aussi le lien pour télécharger le fichier Final_Solution.fig afin qu'ils aient le zoom et la rotation :

Final_Solution.fig

J'ai 2 questions :

1*) J'aimerais savoir si les valeurs du potentiel gravitationnel que vous voyez sur les graphiques sont compatibles avec celle d'une vraie galaxie. J'ai utilisé dans programme les unités CGS (cm gramme seconde).

La galaxie représentée sur les figures fait à peu près 25 kpc de rayon (avec les bras) et qu'il y a 10240 étoiles.

Je n'ai pas trouvé de valeur standard pour l'instant.

J'ai vu en bas du lien suivant ( https://en.wikipedia.org/wiki/Gravit...merical_values ) qu'il faut une
énergie de 130 GJ/kg pour se libérer de l'attraction de la voie lactée, autrement dit son énergie potentielle est égale à 130 GJ pour un objet de 1kg.

Mais comment avoir une estimation sur le potentiel gravitationnel partant de là ?

2*) Sachant que je représente la valeur absolue de ce potentiel gravitationnel, comment se fait-il qu'au niveau des bras, les valeurs en forme de petits pics soient tournées vers le bas, c'est-à-dire vers le haut si on ne prend pas la valeur absolue ?

J'avais tendance à penser qu'une zone en légère sur-densité (comme sur les bras) accentuerait le potentiel gravitationnel, c'est-à-dire augmenterait sa valeur vers le bas (en ne prenant pas la valeur absolue, c'est-à-dire négatif comme pour sa définition de potentiel).

Cependant, j'ai l'impression sur les figures qu'il y a une légère élévation sur les bras (e qui confirmerai mon intuition) mais tout en gardant des pics vers le bas.

Tout aide est la bienvenue.
-----

[Gilgamesh]

J'ai vu en bas du lien suivant ( https://en.wikipedia.org/wiki/Gravit...merical_values ) qu'il faut une
énergie de 130 GJ/kg pour se libérer de l'attraction de la voie lactée, autrement dit son énergie potentielle est égale à 130 GJ pour un objet de 1kg.

Mais comment avoir une estimation sur le potentiel gravitationnel partant de là ?

L'énergie qu'il faut fournir à une masse unitaire pour l'envoyer à l'infini, c'est par définition, le potentiel du champs à cet endroit. Donc tu as ta réponse (à 8 kpc du centre). Vérifie que ça donne le même ordre de grandeur avec ces $*µ£% ! d'unités CGS (à masse galactique égale).

2*) Sachant que je représente la valeur absolue de ce potentiel gravitationnel, comment se fait-il qu'au niveau des bras, les valeurs en forme de petits pics soient tournées vers le bas, c'est-à-dire vers le haut si on ne prend pas la valeur absolue ?

J'avais tendance à penser qu'une zone en légère sur-densité (comme sur les bras) accentuerait le potentiel gravitationnel, c'est-à-dire augmenterait sa valeur vers le bas (en ne prenant pas la valeur absolue, c'est-à-dire négatif comme pour sa définition de potentiel).

Cependant, j'ai l'impression sur les figures qu'il y a une légère élévation sur les bras (e qui confirmerai mon intuition) mais tout en gardant des pics vers le bas.

Tout aide est la bienvenue.

Les petits pics, ce sont tes ~10 k étoiles. Cela forme des puits locaux, le potentiel est plus haut quand tu t'en approches.

[fabio123]

1*) En ce qui concerne la valeur sur le lien wikipedia que j'ai donné ci-dessus, un calcul rapide donne :

avec la masse de la galaxie (comptant 100 milliard d'étoiles) et r la distance depuis laquelle l'objet est amené à l'infini ( 8kpc ) :

on a donc : abs(Phi) = 6.67 1O^-8 * 10^11 * 2 10^33 / (8 kpc) = 13.34 * 10^36 / (8* 10^21) = 1.66 * 10^15 cm^2 s^-2 = 1.66 * 10^11 m^2 s^-2 = 1.66 GJ kg^-1 = 166 GJ kg^-1

ce qui correspond bien à l'ordre de grandeur donné dans la page ( > 130 GJ / kg). En unités CGS, j'obtiens donc la valeur intermédiaire ci-dessus, c'est-à-dire : 1.66 10^15 cm^2 s^-2 .

2*) En ce qui concerne ma solution finale, si je prends une masse solaire ( Ms = 2 10^33 grammes) pour chacune des 10240 étoiles, j'obtiens vers le centre la valeur suivante (en prenant un rayon de 1 kpc avec 1 kpc = 3 10^21 cm)

abs(Phi) = 6.67 10^-8 * 10240 * 2 10^33 / (1 * 3 10^21) = 4,55 10^8 cm^2 s^-2

Il faudrait donc que je retrouve cette valeur au niveau du centre de la galaxie.

Dans mon programme, je résous l'équation de Poisson en utilisant des transformée de Fourier discrètes. Je calcule d'abord la distribution de matière ( rho dans l'équation de Poisson )
comme ceci (en langage C) :

for (i=0; i< numStars; i++)
   {
     compute(node_x, node_y, node_z);
     
     density[node_x][node_y][node_z] = [node_x][node_y][node_z] + mass/((space_x*space_y*space_z)*kpc*kpc*kpc);
     
   }

avec node_x, node_y et node_z les coordonnées de la cellule sur le maillage. space_x, space_y et space_z sont la taille d'une cellule (le coté selon x, y et z) en "kpc".

le facteur "kpc" est la valeur d'un kiloparsec en centimètre (1kpc = 3.08 10^21 cm) et "mass" en grammes. J'ai donc une densité en (grammes cm^-3).

Ensuite, Je discrétise la fonction de Green selon :

for (i=0; i< numStars; i++)
   {
 dx = sqrt(pow((double) (i-Ng/2),2.0))*space_x*kpc;
 dy = sqrt(pow((double) (j-Ng/2),2.0))*space_y*kpc;
 dz = sqrt(pow((double) (k-Ng/2),2.0))*space_z*kpc;

 green_grid[i][j][k] = -1.0/(4.0*M_PI*sqrt(dx*dx + dy*dy + dz*dz));
 
   }

J'ai donc la fonction de Green en cm^-1.

J'applique à partir de là les transformées de Fourier sur la densité et fonction de Green, fais le produit des 2 (sans oublier le facteur 4piG de l'équation de Poisson) et fait une
transformée inverse pour obtenir le potentiel discrétisé

Mais malheureusement, j'obtiens des valeurs très faibles par rapport aux ordres de grandeur cités ci-dessus ( 4,55 10^8 cm^2 s^-2 en comptant 10240 étoiles d'1 masse solaire) : j'obtiens des valeurs
près du centre galactique de l'ordre de 10^-41 cm^2 s^-2 .

Je ne vois pas où est mon erreur.

Est-ce que quelqu'un pourrait me dire si j'ai fait les choses correctement au niveau des facteurs utilisées pour obtenir le potentiel en unités CGS ?

Toute aide est la bienvenue

[fabio123]

Excusez moi, pour le calcul de la densité sur le maillage, la portion de code est :

for (i=0; i< numStars; i++)
   {
     compute(node_x, node_y, node_z); // Calcul des coordonnées de la cellule
     
     density[node_x][node_y][node_z] = density[node_x][node_y][node_z] + mass/((space_x*space_y*space_z)*kpc*kpc*kpc);
     
   }

[A voir en vidéo sur Futura]

[fabio123]

le problème semble être résolu en exprimant aussi la constante de gravitation G en unités CGS.

[Gilgamesh]

Pourquoi ces unités CGS sont encore en usage, dans un domaine où les ordres de grandeurs sont la seule chose qui compte et peuvent se calculer de tête mais s'étendent sur des dizaines d'unités et font appel à des unités complexes et des constantes dimensionnées. Toute chose qui devrait favoriser hautement le SI à l'exclusion de tout autre système fantaisiste.

C'est un des grands mystères de l'univers.

[Calvert]

Pourquoi ces unités CGS sont encore en usage

Une très longue tradition qu'il est impossible de changer sans qu'il y ait une volonté forte derrière. Mais ça n'arrivera pas pendant les prochaines années.

[Gilgamesh]

Une très longue tradition qu'il est impossible de changer sans qu'il y ait une volonté forte derrière. Mais ça n'arrivera pas pendant les prochaines années.

Oui, je me souviens de ta réponse il y a quelques années... Mais tradition ou pas, chacun travaille bien avec les unités qui lui conviennent et j'imagine qu'au moment de publier, personne ne va pas se faire retoquer son article parce qu'il exprime l'énergie d'une supernova en Joule plutôt qu'en erg... Donc je me dis que ça ne gène pas pour travailler au quotidien car sinon le CGS disparaitrait faute d'usagers.

[fabio123]

Je reviens vers vous car j'ai un problème pour la validation des résultats que j'obtiens avec Matlab en utilisant 2 méthodes différentes.

Ces 2 méthodes utilisent un fichier contenant la distribution des particules (des étoiles) avec des distances en kiloparsec (1 parsec = 3.26 années-lumière).

Il y a en tout 10240 particules.

La première méthode est le calcul du potentiel gravitationnel simplement avec la formule suivante :



Partant de la distribution de départ (fichier contenant 3 colonnes pour les coordonnées x,y,z), le calcul se fait de la manière suivante en unités CGS (le potentiel a la dimension de cm^2.s^-2) :

G=6.67428*10^-8; %% cm^3 g^-1 s^-2
kpc=3.085678*10^21; %% 1kpc = 3.08 10^21 cm
mass=2.0*10^33; %% 1 solar mass = 2 10^33 g

% Loading Positions
data=load('Positions_Disk.dat');

% Number of star coordinates : 10240 with Positions_Disk.dat
numStars=numel(data(:,1));

% Number of r values
numPoints=64;

% Min and Max pour Grid
x0=linspace(min(data(:,1)),max(data(:,1)),numPoints);
y0=linspace(min(data(:,2)),max(data(:,2)),numPoints);
z0=linspace(min(data(:,3)),max(data(:,3)),numPoints);

[x y z]=meshgrid(x0,y0,z0);

%% Direct computing of gravitational potential
   for j=1:numPoints
      for k=1:numPoints
         for l=1:numPoints
            Phi(j,k,l)=0.0;
            for i=1:numStars
               Phi(j,k,l)=Phi(j,k,l)-G*mass/((sqrt((x(j,k,l)-data(i,1))^2+(y(j,k,l)-data(i,2))^2+(z(j,k,l)-data(i,3))^2))*kpc);
            end
            Value(j,k,l)=abs(Phi(j,k,l));
         end
      end
   end

% find z index corresponding to z =0
index_z=find(abs(z0)==min(abs(z0)));

% Plot potential
figure(1);
surf(x(:,:,index_z),y(:,:,index_z),Value(:,:,index_z));

Voici la figure que j'obtiens avec cette méthode (potentiel pour z = 0) :

Matlab_Potential_Direct_Computation.jpg

Maintenant, toujours en partant du fichier de distribution, je calcule le potentiel via les transformées de Fourier rapides et la fonction de Green.

D'après Fonction de Green et Transformée de Fourier, on a :

%% Compute the spacing of all cells
Ng=256;
Max_x=max(data(:,1));
Min_x=min(data(:,1));

Max_y=max(data(:,2));
Min_y=min(data(:,2));

Max_z=max(data(:,3));
Min_z=min(data(:,3));

if (Max_x>abs(Min_x))
  space_x=Max_x/(Ng/2);
else
  space_x=abs(Min_x)/(Ng/2);
end

if (Max_y>abs(Min_y))
  space_y=Max_y/(Ng/2);
else
  space_y=abs(Min_y)/(Ng/2);
end

if (Max_z>abs(Min_z))
  space_z=Max_z/(Ng/2);
else
  space_z=abs(Min_z)/(Ng/2);
end

%% Computing discretized Green's function
for i=1:Ng
   for j=1:Ng
      for k=1:Ng
         
         %% Second solution : without space_x, space_y and space_z factors 
         dx = sqrt((i-Ng/2)*(i-Ng/2))*space_x;
         dy = sqrt((j-Ng/2)*(j-Ng/2))*space_y;
         dz = sqrt((k-Ng/2)*(k-Ng/2))*space_z;
         
         %% Units of Green's function : kpc 
         green_grid(i,j,k) = -1.0/(4.0*pi*sqrt(dx*dx + dy*dy + dz*dz));

      end
   end
end

%% For avoiding divergence at Ng/2
green_grid(Ng/2,Ng/2,Ng/2) = 1.0;

%% FFT of Green : FFT(G(r) = G(r) [G(r)] [kpc] => [FFT(G(r)] = value_fft_Green [kpc]^2
fft_green=fftn(green_grid);

%% FFT of density : [FFT(density)] = value_fft_Density [density(x,y,z)].[kpc^3] = value_fft_Density [g].[kpc^-3].[kpc^3] = [g]
fft_density=fftn(density);

%% Product in Fourier space : [G] = value_G [g^-1].[cm^3].[s^-2] = value_G [g^-1].[kpc^3].[s^-2] /(kpc*kpc*kpc)

%% => FFT(Phi) = value_fft_Green * value_fft_Density * value_G [kpc^2].[g].[g^-1].[kpc^3].[s^-2] /(kpc*kpc*kpc)
%%             = value_fft_Green * value_fft_Density * value G [kpc^5].[s^-2] /(kpc*kpc*kpc)
product_fft=4*pi*G*fft_green.*fft_density*factor_fft;

%% Backward of product above

%% Inverse_FFT = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [kpc^-1].[kpc^5].[s^-2] /(kpc*kpc*kpc)
%%       = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [kpc^4].[s^-2] /(kpc*kpc*kpc)
%%       = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [cm^2].[s^-2] * kpc*kpc/(kpc*kpc*kpc)
%%       = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [cm^2].[s^-2] /kpc

fft_potential=real(ifftn(product_fft))/kpc;

%% Shift final solution
fft_potential_shift=fftshift(fft_potential);

%% Plot of potential (absolute value)
figure(2);
surf(x1,y1,abs(fft_potential_shift(:,:,index_z)));

Voici la figure obtenue ( potentiel pour z = 0) :

Matlab_Potential_FFT_Computation.jpg

Comme vous pouvez le voir, il y a toujours une différence assez importante entre les résultats de la première méthode et cette figure (à peu près un facteur 10^4).

D'après ce que j'ai vu, il y a 3 pistes :

1*) Dans la définition de la Fonction de Green, je me suis aperçu qu'elle a la dimension d'une distance même si elle est définie par

Est-ce que je peux simplement l'implémenter comme je l'ai fait ci-dessus par :

%% Computing discretized Green's function
for i=1:Ng
   for j=1:Ng
      for k=1:Ng
         
         %% Second solution : without space_x, space_y and space_z factors 
         dx = sqrt((i-Ng/2)*(i-Ng/2))*space_x;
         dy = sqrt((j-Ng/2)*(j-Ng/2))*space_y;
         dz = sqrt((k-Ng/2)*(k-Ng/2))*space_z;
         
         %% Units of Green's function : distance but cm ? kpc ?
         green_grid(i,j,k) = -1.0/(4.0*pi*sqrt(dx*dx + dy*dy + dz*dz));

      end
   end
end

2*) Est-ce que cet écart de 10^4 résulterait d'un facteur de normalisation dans la FFT du terme de droite dans l'équation de Poisson (4*pi*G*rho), dans la FFT de la fonction de Green et enfin dans la FFT Inverse du produit des 2 FFT précédentes ?

Est-ce que dans Matlab, la FFT inverse d'un produit de 2 FFT doit être divisé par le nombre de points des 2 signaux de départ pour garantir la validité des résultats obtenus ?

3*) Est-ce que le calcul se fait bien en respectant les unités CGS ? En effet, le potentiel final obtenu a la dimension d'une vitesse, soit [cm^2][s^-2]

Mon problème se situe dans l'analyse dimensionnelle de la deuxième méthode (avec les FFT); on a :

%% FFT of Green : FFT(G(r) = G(r) [G(r)] [kpc] => [FFT(G(r)] = value_fft_Green [kpc]^2
fft_green=fftn(green_grid);

%% FFT of density : [FFT(density)] = value_fft_Density [density(x,y,z)].[kpc^3] = value_fft_Density [g].[kpc^-3].[kpc^3] = [g]
fft_density=fftn(density);

%% Product in Fourier space : [G] = value_G [g^-1].[cm^3].[s^-2] = value_G [g^-1].[kpc^3].[s^-2] /(kpc*kpc*kpc)

%% => FFT(Phi) = value_fft_Green * value_fft_Density * value_G [kpc^2].[g].[g^-1].[kpc^3].[s^-2] /(kpc*kpc*kpc)
%%             = value_fft_Green * value_fft_Density * value G [kpc^5].[s^-2] /(kpc*kpc*kpc)
product_fft=4*pi*G*fft_green.*fft_density*factor_fft;

%% Backward of product above

%% Inverse_FFT = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [kpc^-1].[kpc^5].[s^-2] /(kpc*kpc*kpc)
%%       = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [kpc^4].[s^-2] /(kpc*kpc*kpc)
%%       = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [cm^2].[s^-2] * kpc*kpc/(kpc*kpc*kpc)
%%       = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [cm^2].[s^-2] /kpc

[/CODE]

Est-ce que j'ai le droit d'écrire, d'un point de vue dimensionnel :

%% Inverse_FFT = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [kpc^4].[s^-2] /(kpc*kpc*kpc)
%%                      = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [cm^2].[s^-2] * kpc*kpc/(kpc*kpc*kpc)
%%                      = Inverse_FFT_value (value_fft_Green * value_fft_Density * value G) [cm^2].[s^-2] /kpc

??

ce qui ferait que je n'aurais qu'à diviser par le facteur kpc au moment de faire la FFT inverse du produit des 2 FFT précédentes.

Si quelqu'un pouvait m'éclairer sur ces différentes pistes afin de retrouver, avec la méthode utilisant les FFT, les résultats de la première méthode ...

[stefjm]

Pourquoi ces unités CGS sont encore en usage, dans un domaine où les ordres de grandeurs sont la seule chose qui compte et peuvent se calculer de tête mais s'étendent sur des dizaines d'unités et font appel à des unités complexes et des constantes dimensionnées. Toute chose qui devrait favoriser hautement le SI à l'exclusion de tout autre système fantaisiste.

C'est un des grands mystères de l'univers.

Une très longue tradition qu'il est impossible de changer sans qu'il y ait une volonté forte derrière. Mais ça n'arrivera pas pendant les prochaines années.

C'est fait sciemment pour compartimenter le savoir.