Para los cálculos reflejados en las figuras de abajo se ha usado \(r=50.000\ m\), \(D=R-2r\), y como valores de las constantes :
radio terrestre \(R=6.371.000\ m\)
constante de gravitación universal \(G=6\text{,}67 \ 10^{-11} m^3/(kg\ s^2)\)
masa de la Tierra \(M=5\text{,}97 \ 10^{24} kg\)
velocidad angular de la Tierra \(\omega=2 \pi/(23\text{,}9345\times 60\times 60) rad/s\).
La ventaja de semejante modelo es que el potencial gravitacional \(V\) es fácil de calcular: se pone el potencial de la esfera terrestre de radio \(R\) y masa \(M\), y se le resta el que crearía la cavidad si estuviera llena, con una masa \(M (r/R)^3 \), con lo que en un punto fuera de la esfera terrestre tenemos, siendo \(m_r= {\left(\frac{ r }{R }\right)} ^3\) , es
\[V(P)=GM \left(\frac{ -1 }{d }+ \frac{ m_r }{e } \right)\] La superficie equipotencial que pasa a 200 metros sobre el polo norte, el punto en que \(\theta=90º\) y \(d=R+200\), es la mostrada en color azul en las gráficas de la siguiente figura, donde la ordenada es la altura sobre la superficie terrestre, \(d-R\).
En la gráfica de la derecha se ha considerado el potencial efectivo \(W = V + C\), que incluye el potencial centrífugo \[ C(P) = -\frac{ 1 }{2 }(\omega \cdot d \cdot cos\theta )^2 \] La curva equipotencial de W se muestra en color rojo.
Empezemos excluyendo de momento el efecto centrífugo. Si toda la esfera terrestre estuviese rellena homogéneamente la simetría provocaría que la superficie equipotencial pasando a 200m sobre el Polo Norte fuese una esfera perfecta. Debido a la existencia de ese nada despreciable hueco esférico de 50km de radio, centrado a una profundidad de 100km, la superficie equipotencial azul se hunde hasta casi 200m cerca de la oquedad, cuando \(\theta \approx 0º\), y sobresale muy poquito cuando \(\theta \approx 180º\). Lo que parece claro físicamente es que cuando \(\theta \approx 0º\) también deberá ser menor la gravedad \(g\) debido a la ausencia de masa en el hueco. En este caso por tanto un hundimiento de la superficie equipotencial conlleva disminución de \(g\).
Si ahora incluimos el potencial centrífugo \(C\), su influencia es muy marcada, y provoca una gran elevación de la superficie equipotencial de \(W=V+C\) en el ecuador ( \(\theta\approx 0º \ \text{ó} \ \theta \approx 180º\) ) respecto a los polos ( \(\theta\approx 90º \) ). La intuición física nos hace esperar ahora que esa elevación vaya aparejada con una disminución de \(g\), pues cuanto mayor sea la fuerza centrífuga más se contraresta la atracción gravitacional.
Bueno, en un caso es el hundimiento y en el otro la elevación de la superficie equipotencial lo que se corresponde con la disminución de \(g\). Lo que ilustra esta aparente contradicción es que \(g\) no va ligado a subidas o bajadas de la superficie equipotencial, sino a diferencia de distancias entre superficies equipotenciales contiguas. Para calcular aproximadamente \(g\) se buscan dos superficies equipotenciales próximas, por ejemplo las que pasan a 200m y a 300m sobre el Polo Norte. Los valores de \(g\) sin y con el efecto centrífugo se muestran en la siguiente figura. En adelante se expresarán siempre en \(m/s^2\).
Vemos que \(g\) disminuye tanto por el efecto de falta de masa de la oquedad, como por el efecto centrífugo, que es bastante más acusado. Este es el responsable en \(\theta = 180º\) de una disminución de 0,068 en \(g\) comparado con solamente 0,019 en \(\theta = 0º\) debido a la oquedad. Claro que la disminución centrífuga sale casi del doble de lo que uno esperaría considerando el efecto centrífugo sobre la superficie terrestre, \( \omega^2 R \approx 0,034 \). Menudo susto, un factor dos, las cifras que no cuadran, a ver dónde está el error. Bueno, en este caso no hay un 2 extra que sobre, es que \(g\) se calcula sobre la superficie equipotencial, que en el ecuador se eleva 11060 m, con lo que la atracción gravitacional allí se reduce respecto a la superficie terrestre en
\[
\frac{GM}{(R+200)^2} - \frac{GM}{(R+11060+200)^2} \approx 0\text{,}034
\]
¡Vaya! 0,034. De modo que la disminución de 0,068 desde el polo al ecuador es una contribución equitativa de efecto centrífugo "directo" y alejamiento del centro. Es curioso que ese factor 2 no es casual. Olvidando la oquedad, atendiendo al efecto centrífugo exclusivamente, el valor de \(g\) en \(\theta=90º\) viene dado por \(GM/d^2\), y para \(\theta=0º\) viene dado por \(GM/{d'}^2-\omega^2 d'\). Si se restan ambos valores teniendo en cuenta que estamos sobre la misma superficie equipotencial, y que por tanto
\[
- \frac{GM}{d} = - \frac{GM}{d'} - \frac{1}{2} \omega^2 {d'}^2
\]
usamos la expresión de \(1/d\) en términos de \(1/d'\) y ...
\[
\frac{GM}{d^2} - \left(\frac{GM}{{d'}^2} - \omega^2 d' \right) = \frac{\omega^4 {d'}^4}{4GM} + \omega^2 d' + \omega^2 d' \approx 2 \omega^2 d'
\]
puesto que \( \omega^4 {d'}^4 / 4GM = 0\text{,}000029\).
La superficie terrestre tiene cierta plasticidad, y la forma abombada en el ecuador se adapta al efecto centrífugo. En el modelo sencillo utilizado la diferencia entre ecuador y polo norte para la superficie equipotencial es de 11.060m. En realidad la Tierra está muy bien descrita por un elipsoide de radio ecuatorial de 6.378.137 m y radio polar de 6.356.752 m, que supone una diferencia de 21.385 m. Ese elipsoide se diferencia del geoide en pocos cientos de metros. Una vez que hay redistribución de masa el propio potencial gravitacional cambia y no se calcula tan sencillamente.
En cuento a los efectos gravitacionales por exceso o defecto de masa, ya se ve que una oquedad tan inmensa e irreal como la supuesta conlleva un efecto relativamente "débil", de 0,019. No es de extrañar que en geofísica las anomalías de la gravedad a lo largo del terreno se midan en mGal, miligales, siendo un Gal una peculiar designación de \(cm/s^2 = 10^{-2} m/s^2\). Es decir, \(1mGal = 10^{-5}m/s^2\), así de pequeña es la variación local de \(g\).
En el sitio técnico de la misión GOCE hay una enorme cantidad de recursos, en inglés, como el de la gravedad en detalle donde nos presentan entre otras unidades el Gal. El apartado de cantidades del campo de gravedad (Gravity Field Quantities) contiene unas definiciones muy útiles, algunas de las cuales traduzco aquí:
Potencial gravitacional (V): Potencial generado por la atracción de las masas.
Potencial de gravedad (W): Suma del potencial gravitacional (V) y del potencial centrífugo (C) de la Tierra en rotación. Las diferencias entre dos puntos pueden observarse por nivelación.
Superficie equipotencial: Una superficie en la que W es constante. Sus puntos pueden determinarse localmente con mareógrafos, que definen el nivel medio del mar a escala regional.
Geoide: Superficie equipotencial que aproxima el nivel global medio del mar, es decir, una red global de mareógrafos y de señales de nivelación, tras sustracción de los componentes dinámicos. Puede considerarse como un hipotético océano en reposo.
Gravedad: La magnitud, g, del gradiente de W en la superficie terrestre y de V en el espacio. Puede observarse mediante una técnica absoluta (p.ej. experimento de caida libre) o relativa (como una diferencia) mediante un gravímetro de muelle.
He usado esta notación de V, W y C en las anteriores figuras.
En cuanto al código de Octave para sacar semejantes gráficas, primero va esto
#unidades mks global G = 6.67e-11; #cte gravitacion unidades m^3/(kg s^2) global MT=5.97e24; #masa de la Tierra en kg #producto G por MT, unidad de longitud metros global GMT = G*MT; global RT = 6371000; #radio terrestre medio en m global omegaT = (2*pi/(23.9345*60*60));#vel. angular Tierra tomando periodo rot. 23.9345h function val = Vd(x,m,D,theta) global GMT; val = ((-1/x)+(m/sqrt(x^2+D^2-2*x*D*cos(theta))))*GMT; endfunction function val = fdz(xp,m,D,theta) global GMT; C=Vd(xp,m,D,pi/2); fun = @(x) ((-1/x)+(m/sqrt(x^2+D^2-2*x*D*cos(theta))))*GMT - C; [val, fval, info] = fzero(fun,xp+0.000001*rand); endfunction function val = fdzcf(xp,m,D,theta) global GMT; global omegaT; C=Vd(xp,m,D,pi/2); fun = @(x) (((-1/x)+(m/sqrt(x^2+D^2-2*x*D*cos(theta))))*GMT)+ (-0.5*(omegaT*x*cos(theta))^2) - C; [val, fval, info] = fzero(fun,xp+0.000001*rand); endfunction R=RT; r=50000; D=RT-2*r; m=(r/R)^3; theta=(0:180) *pi/180; res1=zeros(1,181); x1=R+200; C1PN=Vd(x1,m,D,pi/2); for k=0:180 res1(k+1)=fdz(x1,m,D,theta(k+1)); endfor res2=zeros(1,181); x2=x1+100; C2PN=Vd(x2,m,D,pi/2); for k=0:180 res2(k+1)=fdz(x2,m,D,theta(k+1)); endfor cfres1=zeros(1,181); for k=0:180 cfres1(k+1)=fdzcf(x1,m,D,theta(k+1)); endfor cfres2=zeros(1,181); for k=0:180 cfres2(k+1)=fdzcf(x2,m,D,theta(k+1)); endfor
y después
subplot(1,2,1)
plot(theta*180/pi,res1-R,'b',90,200,'xk')
set(text(55,190,"equipotencial a 200m"),"color","black")
set(text(55,180,"sobre el Polo Norte"),"color","black")
set(text(55,170,"V (sin efecto centrífugo)"),"color","blue")
set (gca (), "xlim", [0, 180]);
set (gca (), "ylim", [0, 210]);
xlabel(strcat('\theta (grados)'));
ylabel(strcat('d-R (metros)'));
subplot(1,2,2)
plot(theta*180/pi,res1-R,'b',theta*180/pi,cfres1-R,'r',90,200,'xk')
set(text(50,11000,"equipotencial a 200m"),"color","black")
set(text(50,10500,"sobre el Polo Norte"),"color","black")
set(text(5,500,"V (sin efecto centrífugo)"),"color","blue")
set(text(25,10000,"W = V + C (con efecto centrífugo)"),"color","red")
set (gca (), "xlim", [0, 180]);
set (gca (), "ylim", [0, 11500]);
xlabel(strcat('\theta (grados)'));
y para el cálculo y gráfica de \(g\)
#calcula g para caso sin cf
gg=zeros(1,181);
alfa=zeros(1,181);
for k=2:180
alfa(k)=asin((res1(k+1)-res1(k-1))./(res1(k-1)*2));
gg(k)=(C2PN-C1PN)./((res2(k)-res1(k))*cos(alfa(k)));
endfor
gg(1)=(C2PN-C1PN)./(res2(1)-res1(1));
gg(181)=(C2PN-C1PN)./(res2(181)-res1(181));
#calcula g para caso con cf
ggcf=zeros(1,181);
alfa=zeros(1,181);
for k=2:180
alfa(k)=asin((cfres1(k+1)-cfres1(k-1))./(cfres1(k-1)*2));
ggcf(k)=(C2PN-C1PN)./((cfres2(k)-cfres1(k))*cos(alfa(k)));
endfor
ggcf(1)=(C2PN-C1PN)./(cfres2(1)-cfres1(1));
ggcf(181)=(C2PN-C1PN)./(cfres2(181)-cfres1(181));
subplot(1,1,1)
plot(theta*180/pi,gg,'b',theta*180/pi,ggcf,'r',90,gg(91),'xk')
set(text(120,9.813,"V (sin efecto centrífugo)"),"color","blue")
set(text(120,9.795,"W=V+C (con efecto centrífugo)"),"color","red")
title('gravedad en equipotencial a 200m sobre el Polo Norte');
set (gca (), "xlim", [0, 180]);
set (gca (), "ylim", [9.72, 9.82]);
xlabel(strcat('\theta (grados)'));
ylabel(strcat('g (m/s^2)'));
y ... se acabó. No siento los miligales.
