Modelo PREM del interior de la Tierra
Preparación de los datos
Visitamos la página web, titulada 'Data Services Products: EMC-PREM Preliminary Reference Earth Model'. Descargamos el fichero PREM500_IDV.csv, o alternativamente, en este enlace
que como apreciamos en la cabecera (primera y segunda fila) contiene los datos del radio en metros (primera columna) y de la densidad en kg/m3 (segunda columna). Las otras columnas de datos no tienen interés. En total 504 filas de pares de datos (radio, densidad). La densidad de la Tierra es esféricamente simétrica, depende solamente de la distancia (radio) del punto considerado al centro de la Tierra
Situamos el fichero, en la carpeta accesible por MATLAB por defecto, e importamos, utilizando en el menú Import data from file, las dos primeras columnas, poniendo en el control titulado Range A3:B506. Pulsamos el botón titulado Import Selection, eligiendo la opción Import Data.
MATLAB genera dos vectores con nombres por defecto que aparecen en la ventana Wokspace.
Si hemos importado todas las columnas, podemos borrar las no deseadas utilizando el botón derecho del ratón y seleccionado en el menú flotante Delete (Suprimir). En esta ocasión, preferimos operar con números más pequeños, creamos dos nuevos vectores
>> densidad=density/1000; >> radio=indexradius/1e6;
y borramos de Workspace los vectores
Tenemos dos vectores de 504 elementos de longitud cada uno. El primer elemento del vector
>> length(radio) ans = 504 >> radio(1) ans = 0 >> radio(end)*1000 ans = 6371
Aceleración de la gravedad
Representamos la densidad en kg/m3 en el eje vertical en función del radio en km en el eje horizontal
plot(radio*1000,densidad*1000) grid on xlabel('r (km)') ylabel('\rho (kg/m^3)') title('Densidad')

Se ha dividido la Tierra en 503 capas: la primera es una esfera de radio r2 de densidad ρ2, la segunda capa esférica de radio interior r2 y radio exterior r3 de densidad ρ3, la tercera capa esférica de radio interior r3 y radio exterior r4 de densidad ρ4...la última capa esférica de radio interior r503 y radio exterior r504 (radio de la Tierra) de densidad ρ504
La masa de la Tierra es
Como el radio se ha dividido entre 106 y la densidad entre 103, el factor de conversión de unidades es 103(106)3=1021, para obtener la masa de la Tierra en kg
masa=0; for i=2:length(radio) masa=masa+densidad(i)*4*pi*(radio(i)^3-radio(i-1)^3)/3; end display(masa*1e21)
5.9687e+24
La masa Mi contenida en la esfera de radio ri es

La aceleración de la gravedad en el interior de una distribución esférica de masa de radio R, es la fuerza con la que la masa M(r) contenida en la esfera de radio r<R, atrae a la unidad de masa situada a una distancia r del centro, tal como se muestra en la figura. La aceleración de la gravedad tiene dirección radial y apunta hacia el centro (no tendremos en cuenta el movimiento de rotación de la Tierra)
Creamos una nueva tabla de datos (radio ri, aceleración de la gravedad gi). Representamos esta tabla del misma forma que se ha hecho con la densidad.
La aceleración de la gravedad en el exterior de la Tierra, disminuye con el cuadrado de la distancia al centro, siendo M(R) la masa de la Tierra
El factor de conversión de unidades es 6.67·10-11·1021/(106)2=6.67·10-2
%aceleración de la gravedad en el interior masa=0; g=zeros(1,length(radio)); for i=2:length(radio) masa=masa+densidad(i)*4*pi*(radio(i)^3-radio(i-1)^3)/3; g(i)=masa*6.67e-2/radio(i)^2; end hold on plot(radio*1000,g) line([0,radio(end)*1000],[g(end), g(end)],'lineStyle','--','color','k') %aceleración de la gravedad en el exterior de la Tierra f=@(x) 6.67e4*masa./(x.^2); %x en km, 10^-11*10^21/(10^3)^2=10^4 fplot(f,[radio(end)*1000,9000]) line([radio(end)*1000,radio(end)*1000],[0,f(radio(end)*1000)], 'lineStyle','--','color','k') hold off grid on xlabel('r (km)') ylabel('g (m/s^2)') title('Aceleración de la gravedad')
La aceleración de la gravedad en la superficie de la Tierra es
>> g(end) ans = 9.8083
Siguiendo el modelo de capas, dado el radio r<R, calculamos la masa M(r) contenida en dicha esfera. Supongamos que ri<r<ri+1
Definimos una función que devuelve la aceleración de la gravedad g(r) para cualquier valor de r<R
function g = gravedad_interior(r, radio,densidad) masa=0; if r>radio(2) indice=find(radio<r); i=indice(end); for j=2:i masa=masa+densidad(j)*4*pi*(radio(j)^3-radio(j-1)^3)/3; end masa=masa+densidad(i+1)*4*pi*(r^3-radio(i)^3)/3; end g=6.67e-2*masa/r^2; end
Utilizamos esta función para representar, alternativamente, la aceleración de la gravedad en el interior de la Tierra
r=linspace(0,radio(length(radio)),300); g=zeros(1,length(r)); for j=1:length(r) g(j)=gravedad_interior(r(j), radio,densidad); end plot(r*1000,g) grid on xlabel('r (km)') ylabel('g (m/s^2)') title('Aceleración de la gravedad')
Túnel a través de un diámetro
La fuerza que ejerce el campo gravitatorio producido por la Tierra sobre una partícula de masa m es el producto de su masa por la aceleración de la gravedad g(r). La aceleración de la gravedad en el interior de la Tierra es M(r)/r2 y en el exterior es M/r2, siendo M la masa de la Tierra.
El nivel cero de energía potencial se encuentra en el infinito, de modo que la energía potencial de una partícula de masa m situada a una distancia r del centro de la Tierra es
Velocidad

Supongamos que excavamos un túnel a través de un diámetro de la Tierra. Soltamos la partícula en el extremo del túnel, en la superficie de la Tierra, su velocidad cuando está a una distancia r del centro se calcula aplicando el principio de conservación de la energía
Calculamos la velocidad v en función en función de la distancia de la partícula al centro de la Tierra r. Utilizamos la función
function suma=simpson(f,x0,xf,n) %n número par de intervalos, n+1 número de puntos en el vector x=linspace(x0,xf,n+1); h=x(2)-x(1); suma=f(x(1))+f(x(n+1)); for i=2:2:n suma=suma+4*f(x(i)); end for i=3:2:n-1 suma=suma+2*f(x(i)); end suma=suma*h/3; end
Como el radio r está dividido entre 106, hay que multiplicar el resultado de la integral por 106, es decir, la velocidad por 1000
r=linspace(0,radio(end),300); f=@(r) gravedad_interior(r,radio, densidad); v=zeros(1,length(r)); for i=2:length(r) v(i)=sqrt(2*simpson(f,r(i),radio(end),100)); end plot(r(2:end)*1000,v(2:end)*1000) grid on xlabel('r (km)') ylabel('v (m/s)') title('Velocidad')
La velocidad cuando llega al centro de la Tierra es
>> v(2)*1000 ans = 9.9127e+03
Ecuación del movimiento
Situamos el origen en el centro de la Tierra y el eje X tal como se muestra en la figura. Cuando la partícula se encuentra a una distancia x del centro, la fuerza que actúa sobre la partícula es mg(x). La ecuación del movimiento es
Se resuelve la ecuación diferencial por el procedimiento
Primero, definimos la función a integrar, que calcule la aceleración de la gravedad para una distancia r al centro de la Tierra, tarea que hemos realizado en la función
function fg=gravedad_variable(~,x, radio,densidad) masa=0; r=x(1); if r>radio(2) indice=find(radio<r); i=indice(end); for j=2:i masa=masa+densidad(j)*4*pi*(radio(j)^3-radio(j-1)^3)/3; end masa=masa+densidad(i+1)*4*pi*(r^3-radio(i)^3)/3; end fg=[x(2);-6.67e-2*masa/x(1)^2]; end
En segundo lugar, definimos una función para que el procedimiento
function [value,isterminal,direction]=stop_gravedad(t,x) value=x(1); isterminal=1; direction=-1; end
Creamos el script para calcular y representar la posición x de la partícula en el túnel en función del tiempo t
ffg=@(t,x) gravedad_variable(t,x,radio, densidad); opts=odeset('events',@stop_gravedad); [t,x, te]=ode45(fg,[0,4],[radio(end),0], opts); plot(t*1000/60,x(:,1)*1000) grid on xlabel('t (min)') ylabel('r (km)'); title('Túnel')
El tiempo que tarda en llegar al centro de la Tierra en minutos y su velocidad son, respectivamente
>> te*1000/60 ans = 19.0930 >> x(end,2)*1000 ans = -9.9515e+03
El perido de una oscilación es cuatro veces este tiempo, 76.4 min
Túnel a través de una cuerda

Cuando la partícula se mueve por un túnel excavado a lo largo de una cuerda distante d del centro, la componente de la fuerza sobre la partícula a lo largo del eje X es
La ecuación del movimiento se escribe
Se resuelve la ecuación diferencial por el procedimiento
Primero, definimos la función a integrar, que calcule la aceleración de la gravedad para una distancia r al centro de la Tierra, tarea que hemos realizado al definir las funciones
function fg=gravedad_variable_1(~,x, radio,densidad, d) r=sqrt(x(1)^2+d^2); masa=0; if r>radio(end) r=radio(end); end if r>radio(2) indice=find(radio<r); i=indice(end); for j=2:i masa=masa+densidad(j)*4*pi*(radio(j)^3-radio(j-1)^3)/3; end masa=masa+densidad(i+1)*4*pi*(r^3-radio(i)^3)/3; end fg=[x(2);-6.67e-2*masa*x(1)/r^3]; end
En segundo lugar, utilizamos la función
Creamos el script para calcular y representar la posición x de la partícula en el túnel en función del tiempo t. El túnel dista R/2 del centro de la Tierra
d=0.5*radio(end); fg=@(t,x) gravedad_variable_1(t,x,radio, densidad, d); opts=odeset('events',@stop_gravedad); [t,x, te]=ode45(fg,[0,4],[sqrt(radio(end)^2-d^2),0], opts); plot(t*1000/60,x(:,1)*1000) grid on xlabel('t (min)') ylabel('r (km)'); title('Túnel')
El tiempo que tarda en llegar al centro del túnel en minutos y su velocidad son, respectivamente
>> te*1000/60, x(end,2)*1000 ans = 19.7047 ans = -8.0039e+03
Modificamos el script para calcular el tiempo de llegada de la partícula al centro del túnel para distintas distancias d del túnel al centro de la Tierra
opts=odeset('events',@stop_gravedad); j=1; distancias=0:0.05:0.95; T=zeros(1,length(distancias)); for d=distancias*radio(end) fg=@(t,x) gravedad_variable_1(t,x,radio, densidad, d); [t,x, te]=ode45(fg,[0,4],[sqrt(radio(end)^2-d^2),0], opts); T(j)=te; j=j+1; end plot(distancias*radio(end)*1000,T*1000/60,'-o','markersize',3,'markeredgecolor', 'r','markerfacecolor','r') grid on xlabel('r (km)') ylabel('T (min)'); title('Túnel')
En el modelo de densidad constante del interior de la Tierra todos los tiempos eran iguales, en el modelo PREM los tiempos son diferentes, siendo mínimo cuando el túnel está excavado a través de un diámetro
Aceleración de la gravedad. Ajuste

La forma de la curva de la acelaración de la gravedad en función de la distancia al centro de la Tierra, nos sugiere la posibilidad de dividirla en dos capas, en la primera la aceleración de la gravedad crece linealmente con el radio y en la segunda, la aceleración de la gravedad disminuye linealmente con el radio. Las ecuaciones de las rectas que se ha dibujado en color rojo en la figura son
Buscamos los valores de los parámetros x1 y ζ1 que hacen que la función ζ(x) se ajuste a los valores calculados de la aceración de la gravedad g(r)/g(R) en función de x=r/R. g(R) es la aceleración de la gravedad en la superficie de la Tierra. Véase la página titulada Ajuste de datos (no lineal), la sección 'Ajuste mediante
masa=0; g=zeros(1,length(radio)); for i=2:length(radio) masa=masa+densidad(i)*4*pi*(radio(i)^3-radio(i-1)^3)/3; g(i)=masa*6.67e-2/radio(i)^2; end hold x=radio'/radio(end); %vector fila y=g/g(end); %vector fila plot(x,y) %ajuste f=@(a,x) a(1)*(x<a(2)).*x/a(2)+(x>=a(2)).*(((1-a(1))*x+a(1)-a(2))/(1-a(2))); error=@(a) sum((y-f(a,x)).^2); a0=[1,0.5]; %valor inicial af=fminsearch(error,a0); g=@(x) f(af,x); fplot(g,[0,1]) line([af(2),af(2)],[0,af(1)],'lineStyle','--', 'color','k') hold off grid on xlabel('r/R') ylabel('g/g(R)') title('Aceleración de la gravedad')
Los valores de los parámetros ζ1 y x1 son
>> af af = 1.0514 0.4869
Túnel a través de una cuerda
Velocidades
Calculamos la velocidad adimensional
Hay dos posibilidades
Para x1≤x≤1
Para 0≤x≤x1
Vamos a buscar, soluciones analíticas para calcular el tiempo que tarda una partícula en recorrer un túnel que atraviesa la Tierra
Primer caso, δ>x1
La partícula se deja caer en el extremo del túnel. El tiempo T que tarda en recorrer el túnel es
Cambiamos la variable y por la variable r,
Introducimos la expresión de la velocidad adimensional , válida en el intervalo x1≤x≤1
En el libro titulado 'Table of Integrals, Series, and Products' buscamos la solución a integrales de este tipo y la encontramos en el apartado 3.148, n° 5 de la página 276
Identificamos las variables
El resultado es
El tiempo T se expresa en términos de las integrales elípticas completas de primera especie y de tercera, F(π/2,q)=K(q),
Segundo caso, δ<x1
La partícula se deja caer en el extremo del túnel. El tiempo T que tarda en recorrer el túnel es
Cambiamos la variable y por la variable r,
Introducimos la expresión de la velocidad adimensional , válida en el intervalo x1≤x≤1 y de , válida en el intervalo 0≤x≤x1
La primera integral, es similar al caso anterior. Denominamos
El resultado es
La segunda integral se resuelve haciendo el cambio u=x2, es similar a la trayectoria de un planeta.
Hacemos un nuevo cambio de variable
Con este cambio la integral es inmediata
Evaluamos el integrando entre x1 y δ. Denominamos
El resultado final es
El segundo término entre corchetes, se puede simplificar aún más, del siguiente modo:
El resultado es
Resultados
Comparamos los tiempos adimensionales en recorrer el túnel mediante cálculo numérico y mediante la aproximación analítica que hemos descrito en este apartado
distancias=0:0.05:0.95; g1=1.0514; x1=0.4869; mu=(g1-1)/(1-x1); T=zeros(1,length(distancias)); j=1; hold on for d=distancias q=sqrt((1-d)*(mu+2+d*mu)/((mu+2-mu*d)*(1+d))); A=1/sqrt((mu+2-mu*d)*(1+d)); if d>x1 T(j)=4*A*((mu+2)*ellipke(q^2)/mu-2*ellipticPi(mu*(1-d)/(mu+2-mu*d),q^2)/mu); else eta=1+g1-x1; k=asin(sqrt((mu+2-mu*d)*(1-x1)/((1-d)*(mu+2-mu*x1)))); T(j)=sqrt(x1/g1)*(pi/2-asin((g1*(2*x1^2-d^2)-eta*x1)/(g1*d^2-eta*x1)))+ 4*A*((mu+2)*ellipticF(k,q^2)/mu-2*ellipticPi(mu*(1-d)/(mu+2-mu*d),k,q^2)/mu); end j=j+1; end plot(distancias,T,'-bo','markersize',3,'markeredgecolor','b','markerfacecolor','b') TT=zeros(1,length(distancias)); opts=odeset('events',@stop_gravedad); j=1; for d=distancias*radio(end) fg=@(t,x) gravedad_variable_1(t,x,radio, densidad, d); [t,x, te]=ode45(fg,[0,4],[sqrt(radio(end)^2-d^2),0], opts); TT(j)=2*te*1000/sqrt(radio(end)*1e6/9.8083); j=j+1; end plot(distancias,TT,'-ro','markersize',3,'markeredgecolor','r','markerfacecolor','r') hold off grid on xlabel('\delta') ylabel('\tau'); legend('analítica','numérica','Location','best') title('Tiempos')
Referencias
Alexander R. Klotz. The gravity tunnel in a non-uniform Earth. Am. J. Phys. 83(3) march 2015, pp. 231-237
Stefan Isermann. Analytical solution of gravity tunnels through an inhomogeneous Earth. Am. J. Phys. 87(1), January 2019, pp. 10-17
I.S. Gradshteyn, I.M. Ryzhik. Table of Integrals, Series, and Products. Seventh Edition. Elsevier, 2007