Movimiento bajo una fuerza central y una perturbación (II)

Supongamos un planeta de masa M que está rodeado de una capa de polvo, de densidad uniforme ρ. Un satélite artificial de masa m se mueve alrededor del planeta. Vamos a estudiar su movimiento bajo la fuerza de atracción del planeta y de la capa de polvo. No consideraremos el posible rozamiento del satélite al chocar con las partículas de polvo
La fuerza de atracción del planeta es inversamente proporcional al cuadrado de la distancia
Demostramos en la página titulada La aceleración de la gravedad en el interior y en el exterior de una distribución esférica y uniforme de masa que la fuerza de atracción de un distribución esférica y uniforme de masa ejerce sobre una partícula situada a una distancia r del centro, es la que ejercería sobre dicha partícula, la masa contenida en una esfera concéntrica de radio r.
El satélite se mueve bajo la acción de dos fuerzas centrales, una inversamente proporcional al cuadrado de la distancia r, y otra más pequeña, proporcional a la distancia r al centro fijo de fuerzas. Las ecuaciones del movimiento en coordendas polares son
La segunda ecuación nos indica que el momento angular L, expresado en coordenadas polares, es constante
Teniendo en cuenta la constancia del momento angular L, escribimos la ecuación diferencial que describe el movimiento en la dirección radial
Escribimos esta ecuación en términos de la densidad relativa del polvo k=ρ/ρT respecto de la densidad del planeta y del radio del planeta rT. Recordando que la masa del planeta M es
El sistema de ecuaciones diferenciales que tenemos que resolver es
Ecuaciones de Lagrange
La fuerza de atracción -GMm/r2 es conservativa y su energía potencial es -GMm/r. La fuerza atractiva -4πGρr/3 es similar a la fuerza que ejerce una muelle elástico -kx sobre una partícula, la energía potencial de esta fuerza es kx2/2. Por tanto, la fórmula de la energía potencial para esa fuerza es 2πGρr2/3
La energía cinética en coordenadas polares y la energía potencial del satélite de masa m son
La energía total es la suma de las dos cantidades. La lagrangiana L (no confundir con el momento angular) es la diferencia
Obtenemos las mismas ecuaciones
Resolvemos el sistema de dos ecuaciones diferenciales, con las siguientes condiciones iniciales: en el instante t=0, r=r0, dr/dt=0.

Nos queda determinar el valor del momento angular L que es una constante del movimiento. Supondremos que el satélite describe una órbita circular de radio r0 alrededor del planeta de masa M, no cubierto de la nube de polvo.
Su velocidad v0, se obtiene de la ecuación de la dinámica del movimiento circular
El momento angular L será
De esto modo, comprobaremos cómo se modifica la órbita circular del satélite artificial por el efecto de la nube de polvo.
Escalas
Consideraremos que el planeta es la Tierra, cuya masa M=5.98·1024 kg y cuyo radio es rT=6.37·106 m. La constante G=6.67·10-11 Nm2/kg2. Para efectuar los cálculos, establecemos un sistema de unidades de modo que el tiempo τ se mide en horas, t=3600·τ, y la distancia R en unidades de 1000 km, r=106R. Las ecuaciones del movimiento se escriben
Resolvemos el sistema de ecuaciones diferenciales, con las siguientes condiciones iniciales: en el instante τ=0, R=R0 y dR/dτ=0
RT=6.37; %radio de la Tierra R0=8; %Posición inicial del satálite k=0.05; %densidad relativa del polvo x0=[R0,0,0]; %condiciones iniciales %x(1) es r, x(2) es dr/dt, x(3) es theta fg=@(t,x) [x(2);5169.303*(R0/x(1)^3-1/x(1)^2-k*x(1)/RT^3); 71.89787*sqrt(R0)/x(1)^2]; [t,x]=ode45(fg,[0,4],x0); plot(x(:,1).*cos(x(:,3)), x(:,1).*sin(x(:,3))) grid on axis equal xlabel('x') ylabel('y') title('Movimiento en una nube') figure subplot(2,1,1) plot(t,x(:,1)) grid on xlabel('\tau') ylabel('R') subplot(2,1,2) plot(t,x(:,3)) for i=1:floor(x(end,3)/(2*pi)) line([0,t(end)],[i*2*pi,i*2*pi],'lineStyle','--','color','k') end grid on xlabel('\tau') ylabel('\theta')
Se representa la trayectoria hasta el instante τ=4
Si no hubiese polvo de densidad relativa k, el satélite describiría una órbita circular de radio R0=8, con un periodo
Con r0=106R0. En las unidades que hemos establecido, Pθ=1.9774 horas
Comprobamos que la energía del satélite se mantiene constante, véase el apartado más abajo
>> E=x(:,2).^2/2+5169.3034*(R0./(x(:,1).^2*2)-1./x(:,1)+x(:,1).^2*k/(2*RT^3)) E = -291.0826 -291.0826 -291.0826 .....
Representamos en la parte superior de la figura, la distancia radial R en función del tiempo τ. El periodo radial Pr, es el intervalo de tiempo entre dos máximos o dos mínimos
Representamos en la parte inferior de la figura, la posición angular θ en función del tiempo τ. El periodo Pθ, es el tiempo que tarda la posición angular en incrementarse 2π radianes. La primera línea a trazos marca la posición 2π, la segunda 4π, el periodo es el intervalo de tiempo entre dos intersecciones consecutivas
Energía
La suma de la energía cinética y potencial es
Sabiendo que el momento angular , la energía E/m se escribe
En el sistema de unidades que hemos establecido,
Energía potencial efectiva
Cuando se estudia el movimiento en la dirección radial, la energía potencial efectiva es
Con las condiciones iniciales especificadas, t=0, r=r0 y dr/dt=0, , la energía potencial efectiva es
En la segunda fila expresamos la energía potencial efectiva en el sistema de unidades que hemos establecido
La energía es constante e igual a su valor en el instante inicial t=0, r=r0 y dr/dt=0
Representamos, la energía potencial efectiva eef, y superponemos la energía total e mediante una línea horizontal. Los puntos de intersección señalan las distancias R máxima y mínima del satélite al centro de fuerzas. Son las raíces reales de un polinomio de cuarto orden. Utilizamos la función roots de MATLAB para calcularlas y representarlas mediante líneas a trazos.
La energía potencial efectiva tiene un mínimo, que se obtiene derivando eef respecto de R e igualando a cero
Utilizamos la función roots de MATLAB para calcular la raíz real positiva.
RT=6.37; %radio de la Tierra R0=8; %Posición inicial del satálite k=0.05; f=@(R) R0./(2*R.^2)-1./R+k*R.^2/(2*RT^3); fplot(f,[RT,R0+1]) line([RT,R0+1],[f(R0),f(R0)],'color','k') raices=roots([k/(2*RT^3),0,-f(R0),-1,R0/2]) raiz=zeros(2,1); j=1; for i=1:length(raices) if isreal(raices(i)) raiz(j)=raices(i); j=j+1; end end %mínimo raices=roots([k/RT^3,0,0,1,-R0]); for i=1:length(raices) if isreal(raices(i)) & raices(i)>0 Rmin=raices(i); end end for j=1:2 line([raiz(j),raiz(j)],[f(Rmin), f(R0)],'lineStyle','--') end grid on xlabel('R') ylabel('e_{ef}') title('Energía potencial efectiva')
Aproximación para órbitas casi circulares.
Partimos de la ecuación diferencial en la dirección radial. Hacemos que r=r0+ξ,
Suponiendo que ξ/r0<<1, hacemos la aproximación
La ecuación diferencial se escribe
Si r0 es el radio de la órbita circular, d2r/dt2=0
obtenemos una ecuación diferencial sencilla
Que es la ecuación diferencial de un Movimiento Armónico Simple. La frecuencia angular de la oscilación radial es,
y el periodo Pr=2π/ωr
Para , densidad relativa k=0.05, r0=8·106 m y rT=6.37·106 m, obtenemos Pr=1.7362 h. La distancia entre dos mímimos de R en la gráfica es aproximadamente 1.6 h
El tiempo Pθ que tarda en dar una vuelta completa,
Pθ=2π/ωθ
Con los datos anteriores, Pθ=1.9774 h
Ambos periodos Pr y Pθ son distintos, lo que da lugar a la precesión de la órbita
Referencias
Lim Yung-kuo. Problems and Solutions on Mechanics. World Scientific (1994). Problem 1047, pp. 66-69