La forma de la Tierra

Una superficie equipotencial

Un objeto situado sobre la superficie de la Tierra a una latitud λ, describe un movimiento circular de radio x=r·cosλ tal como se muestra en la figura.

Las fuerzas que actúan sobre dicho objeto desde el punto de vista de un observador no inercial que se mueve con la Tierra son:

La resultante de ambas fuerzas es la tensión T de la cuerda que sujeta la partícula de masa m.

Fuerzas conservativas. Energías potenciales

La fuerza de atracción Fg es conservativa y su energía potencial es

E g =G Mm r

De nuevo, se supone que la Tierra es aproximadamente esférica, y no se tiene en cuenta el abultamiento producido en el ecuador como consecuencia del movimiento de rotación de la Tierra.

La fuerza centrífuga depende únicamente de la distancia x al eje de rotación es por tanto, una fuerza conservativa, cuya energía potencial es

0 x m ω 2 x ·dx= E c (0) E c (x)

Tomando el nivel cero de la energía potencial en el eje de rotación x=0.

E c (x)= 1 2 m ω 2 x 2 = 1 2 m ω 2 r 2 cos 2 λ

La energía potencial total es la suma de ambas contribuciones

E p (r)= 1 2 m ω 2 r 2 cos 2 λG Mm r

La forma de la Tierra es la de una superficie equipotencial, ya que la dirección vertical o la dirección de la intensidad de la gravedad efectiva g es perpendicular a la superficie equipotencial en cada punto de dicha superficie.

Sea r=b, cuando λ=π/2, la energía potencial vale Ep=-GMm/b. Así que el valor de b determina la superficie equipotencial que describe la forma de la superficie de la Tierra.

r 3 cos 2 λ 2GM ω 2 b r+ 2GM ω 2 =0 .

Sea r=a, cuando λ=0. La raíz de la ecuación cúbica determina a en función de b, G, M y ω.

a 3 2GM ω 2 b a+ 2GM ω 2 =0

Conocido b calculamos a resolviendo la ecuación cúbica.

En la siguiente tabla, se proporcionan los datos relativos a la masa en kg, periodo de rotación en horas, y los valores de los radios a (ecuatorial) y b (polar) de algunos planetas del Sistema Solar.

Planeta Masa (kg) Periodo (h) b observado (m) a observado (m)
Tierra 0.598·1025 23.93 6.356·106 6.378·106
Marte 0.0658·1025 24.62 3.40·106 3.417·106
Júpiter 190·1025 9.9 66.93·106 71.35·106
Saturno 57·1025 10.2 54.60·106 60.40·106
Urano 9·1025 10.8 22.37·106 23.80·106
Neptuno 10·1025 15.8 21.76·106 22.20·106

Con el dato de G=6.67·10-11 Nm2/kg2, la velocidad angular de rotación de la Tierra es

ω= 2π 23.93·3600 =7.293· 10 5 rad/s

Para facilitar el cálculo en el ordenador, expresamos la distancia en unidades de 106 m. La ecuación que nos permite calcular a a partir del valor observado de b=6.356 es

a3-23594.12a+149964.23=0

Para calcular las raíces de esta ecuación cúbica sugerimos dos procedimientos:

Raíces de una ecuación cúbica

function esferoide_3
    M=0.598e25; %masa de la Tierra
    G=6.67e-11; %constante G
    w=2*pi/(23.93*3600); %velocidad angular de rotación
    b=6.356; %semieje b en 10^6 m
    p1=2*G*M*1e-18/(w^2*b);
    p0=2*G*M*1e-18/w^2;
    raiz=raices_3([1,0,-p1,p0]);
    disp(raiz)

    function x = raices_3(p)
        Q=(p(2)*p(2)-3*p(3))/9;
        R=(2*p(2)^3-9*p(2)*p(3)+27*p(4))/54;
        x=zeros(3,1); %reserva memoria para un vector de tres elementos
        if (R*R)<(Q^3)
            tetha=acos(R/sqrt(Q^3));
            x(1)=-2*sqrt(Q)*cos(tetha/3)-p(2)/3;
            x(2)=-2*sqrt(Q)*cos((tetha+2*pi)/3)-p(2)/3;
            x(3)=-2*sqrt(Q)*cos((tetha-2*pi)/3)-p(2)/3;
        else
            A=-sign(R)*nthroot(abs(R)+sqrt(R*R-Q^3),3);
            if A==0
                B=0;
            else
                B=Q/A;
            end
            x(1)=(A+B)-p(2)/3;
            x(2)=-(A+B)/2-p(2)/3+(sqrt(3)*(A-B)/2)*sqrt(-1); %mejor que i
            x(3)=-(A+B)/2-p(2)/3-(sqrt(3)*(A-B)/2)*sqrt(-1);
        end
    end
end
 -156.6882
  150.3213
    6.3669

La raíz buscada es a=6.3669·106 m

Alternativamente, se puede utilizar la función roots de MATLAB en vez de la función raices_3

Procedimiento iterativo

Esta ecuación se puede transformar en esta otra equivalente

a= a 3 +149964.23 23594.12

para hallar la raíz por el método de iteración partiendo de un valor inicial próximo a 6, obteniendo el valor a=6.367.

x0=6;
while(1)
    x1=(x0*x0*x0+149964.23)/23594.12;
    if abs(x1-x0)<0.001    
        break;
    end
    x0=x1;
end
disp(x0)
    6.3669

Forma del planeta

En la segunda ecuación despejamos 2GM/ω2

a 3 2GM ω 2 b a+ 2GM ω 2 =0 2GM ω 2 = a 3 b ab

Sustituimos en la primera

r 3 cos 2 λ 1 b ( a 3 b ab )r+ a 3 b ab =0

Conocidos los semiejes a y b, dado el ángulo 0≤λ<π/2 calculamos la raíz real de la ecuación cúbica en r

Representamos la forma del planeta calculado la abscisa x=r·cosλ y la ordenada y=r·sinλ de cada punto (color rojo)

Comparamos con la ecuacion de la elipse (color azul)

x 2 a 2 + y 2 b 2 =1

Vemos que coinciden. Calculando el primer cuadrante, extendemos la representación gráfica a los tres restantes

function esferoide_2
    M=0.598e25; %masa de la Tierra
    G=6.67e-11; %constante G
    w=2*pi/(23.93*3600); %velocidad angular de rotación
    b=6.356; %semieje b en 10^6 m
    p1=2*G*M*1e-18/(w^2*b);
    p0=2*G*M*1e-18/w^2;
    raiz=roots([1,0,-p1,p0]);
    a=raiz(3);
    z=@(x) b*sqrt(1-x.^2/a^2);
    fp=fplot(z,[0,a]);
    plot(fp.XData, fp.YData, 'r', fp.XData, -fp.YData,'r', -fp.XData,
 fp.YData,'r', -fp.XData, -fp.YData, 'r')
    Th=(0:89)*pi/180;
    xx=zeros(1,length(Th));
    yy=zeros(1,length(Th));
    k=1;
    for th=Th
        p1=-a^3/((a-b)*cos(th)^2);
        p0=a^3*b/((a-b)*cos(th)^2);
        raiz=roots([1,0,p1,p0]);
        xx(k)=raiz(3)*cos(th);
        yy(k)=raiz(3)*sin(th);
        k=k+1;
    end
    plot(xx,yy,'b',-xx,yy,'b', xx,-yy,'b',-xx,-yy,'b')
    hold off
    grid on
    axis equal
    xlabel('x,y')
    ylabel('z')
    title('Tierra')
end

Otra ecuación que describe la forma de la Tierra

La ecuación de la elipse es

x 2 a 2 + y 2 b 2 =1 r 2 sin 2 θ a 2 + r 2 cos 2 θ b 2 =1

Se define el coeficiente de achatamiento f

f= ab a

Para la Tierra a=6 378 137 m y b=6 356 752 m, lo que da un coeficiente f=3.35·10-3

La ecuación de la elipse se convierte en

r 2 = a 2 ( 1f ) 2 cos 2 θ+ ( 1f ) 2 sin 2 θ = a 2 ( 1f ) 2 ( 1f sin 2 θ ) 2 + f 2 sin 2 θ cos 2 θ r a = 1f 1f sin 2 θ 1 1+ f 2 sin 2 θ cos 2 θ ( 1f sin 2 θ ) 2

Teniendo en cuenta que para la Tierra f es muy pequeño

r a 1f 1f sin 2 θ ( 1 1 2 f 2 sin 2 θ cos 2 θ ( 1f sin 2 θ ) 2 ) r a 1f 1f sin 2 θ ( 1f )( 1+f sin 2 θ+... )1f( 1 sin 2 θ )=1f cos 2 θ

>> syms x;
>> taylor(1/sqrt(1+x^2),x)
ans =(3*x^4)/8 - x^2/2 + 1
>> taylor(1/(1-x),x)
ans =x^5 + x^4 + x^3 + x^2 + x + 1

La definción del polinomio P2(x) de Legendre

P 2 ( cosθ )= 1 2 ( 3 cos 2 θ1 ) r a 1 f 3 2 3 f P 2 ( cosθ )

El radio medio de la esfera equivalente, se obtiene igualando el volumen del elipsoide de revolución de semiejes a y b con el volumen de una esfera de radio Rt

4 3 π R t 3 = 4 3 π a 2 b R t 3 = a 3 ( 1f ) R t =a ( 1f ) 1/3 a( 1 f 3 )

La ecuación de la forma del planeta se escribe

ra( 1 f 3 ) 2 3 af P 2 ( cosθ ) r R t ( 1 2 3 f P 2 ( cosθ ) ) r R t ( 1+ 1 3 ff cos 2 θ )

Esta es la primera aproximación a la forma de la Tierra. La segunda ecuación, expresada en términos del polinomio de Legendre, nos será útil para calcular la energía potencial producida por un esferoide homogéneo

El semieje mayor a, se obtiene para θ=π/2 y el b, para θ=0

b R t ( 1+ 1 3 ff )= R t ( 1 2 3 f ) a R t ( 1+ 1 3 f )

Comparamos la ecuación de la elipse con la expresión polar aproximada, para f=0.1 y 0.3, tomando Rt=1

x 2 a 2 + y 2 b 2 =1 r R t ( 1+ 1 3 ff cos 2 θ ),{ x=rsinθ y=rcosθ f= ab a , R t = ( a 2 b ) 1/3

Rt=1;
f=0.1;
r=@(th) Rt*(1-f*cos(th).^2+f/3);
hold on
fplot(@(th) r(th).*sin(th), @(th) r(th).*cos(th),[0,2*pi] )
a=Rt*(1+f/3);
b=Rt*(1-2*f/3);
z=@(x) b*sqrt(1-x.^2/a^2);
fp=fplot(z,[0,a]);
plot(fp.XData, fp.YData, 'r', fp.XData, -fp.YData, 'r', -fp.XData,
 fp.YData, 'r', -fp.XData, -fp.YData, 'r')
hold off
grid on
axis equal
xlabel('x,y')
ylabel('z')
title('esferoide')

Para f=0.3, hay diferencia

Para f=0.1, casi coinciden

En el caso de la Tierra a=6 378 137 m y b=6 356 752 m

a=6378137/1e6;
b=6356752/1e6; %en unidades 10^6 m
Rt=(a^2*b)^(1/3); %radio medio de la esfera equivalente
f=(a-b)/a; %achatamiento
r=@(th) Rt*(1-f*cos(th).^2+f/3);
hold on
fplot(@(th) r(th).*sin(th), @(th) r(th).*cos(th),[0,2*pi] ) %aproximada
%circunfeencia de radio Rt
fplot(@(th) Rt*sin(th), @(th) Rt*cos(th),[0,2*pi],'lineStyle','--' )
%circunfeencia de radio Rt
a=Rt*(1+f/3);
b=Rt*(1-2*f/3);
z=@(x) b*sqrt(1-x.^2/a^2);
fp=fplot(z,[0,a]);
plot(fp.XData, fp.YData, 'r', fp.XData, -fp.YData, 'r', -fp.XData, 
fp.YData, 'r', -fp.XData, -fp.YData, 'r') %elipse
hold off
grid on
axis equal
xlabel('x,y, 10^6 m')
ylabel('z, 10^6 m')
title('Tierra')

La Tierra se desvía muy poco de la forma esférica. Utilizando varias veces la herramienta Zoom podremos ver la diferencia entre la circunferencia de radio Rt y la elipse de semiejes a y b

Referencias

Bolemon. Shape of the rotating planets and the Sun: A calculation for elementary mechanics. Am. J. Phys. 44 (11) November 1976, pp. 1125-1128.

William Lowrie. A Student’s Guide to Geophysical Equations. Cambridge University Press, (2011), pp. 86-90