Tyngdekraft

Man kan simulere tyngde-interaktionen mellem massive objekter, ved at bestemme kraftvektoren mellem hvert objekt og bestemme den resulterende effekt på objekts position med numerisk integration.

To-legeme system

Tyngdekraften mellem to objekter er beskrevet af Newtons ligning: \begin{equation} \vec{F}_g = G \frac{M \ m}{ \left| \vec{r} \right|^2} \hat{r} \end{equation} Hvor \(M\) og \(m\) er massen af de respektive objekter, \(G\) er gravitationskonstanten (\( G = 6.67 \times 10^{-11} N \ m^2 \ kg^{-2} \)), \(\vec{r}\) er afstandsvektoren mellem de to objekter og \(\hat{r}\) er retningsvektoren mellem objekterne.
Hvis vi antager at det ene objekt er meget større end det andet (f.eks. Jordens masse \(m_\oplus\) i forhold til en satellits masse \( m_{sat} \)), og derfor ikke ændrer sin position under simuleringen, kan vi nøjes med at følge det lille objekts bane. For hvert tidsskridt beregner vi følgende udtryk:

Positionsopdatering af satellit i kredsløb om Jorden

Newtonisk tyngdekraft

For hvert tidsskridt \(\Delta t\) beregn, \begin{align} 1) \ & \vec{r} = \vec{r}_{sat_n} - \vec{r}_{\oplus_n} & \text{Afstandsvektor mellem de to objekter} \\ 2) \ & \hat{r} = \frac{\vec{r}}{\left| \vec{r} \right|} \ & \text{Retningsvektor mellem de to objekter} \\ 3) \ & \vec{F}_g = - G \frac{m_{\oplus} \ m_{sat}}{ \left| \vec{r} \right|^2} \hat{r} & \text{Kraften mellem de to objekter} \\ 4) \ & \vec{p}_{sat_{n+1}} = \vec{p}_{sat_n} + \Delta t \ \vec{F}_g & \text{Impulsvektor opdatering af satellitten} \\ 5) \ & \vec{r}_{sat_{n+1}} = \vec{r}_{sat_n} + \Delta t \ \frac{\vec{p}_{sat_{n+1}}}{m_{sat}} & \text{Positionsopdatering af det satellitten} \end{align}

\begin{align} 6^*) \ & \vec{p}_{\oplus_{n+1}} = \vec{p}_{\oplus_n} + \Delta t \ \vec{F}_g & \text{Impulsvektor opdatering af Jorden} \\ 7^*) \ & \vec{r}_{\oplus_{n+1}} = \vec{r}_{\oplus_n} + \Delta t \ \frac{\vec{p}_{\oplus_{n+1}}}{m_{\oplus}} & \text{Positionsopdatering af Jorden} \end{align}

* Det primære objekt (Jorden i dette eksempel) kan ligeledes opdateres. Da kraften er den samme (Newtons 3. lov) kan vi nøjes med at opdatere objektets impuls og positionsvektorer

r_jorden = [0, 0, 0]; % (m) Position af Jorden
r_sat = [9e7, 0, 0]; % (m) Position af satellit

m_sat = 1500; % (kg) Satellittens masse
m_jorden = 5.97e24; % (kg) Jordens masse

v_sat = [0, 2400, 600]; % (m/s) begyndelseshastighed for satellitten
p_sat = m_sat .* v_sat; % (kg*m/s) satellittens impuls 

G = 6.67e-11; % (N m^2 kg^-2) Gravitationskonstanten

dt = 1000; % (s) tidsskridt størrelse
t = 0; % (s) Begyndelses tid 
tf = 360*24*60; % (s) slut tid

% Plot startbetingelser
plot3(r_jorden(1), r_jorden(2), r_jorden(3),'ok')
hold on
plot3(r_sat(1), r_sat(2), r_sat(3), 'o')
xlabel('x (m)')
ylabel('y (m)')
zlabel('z (m)')

% Løkke (Euler-integration)
while t < tf
    r = r_sat - r_jorden; % (afstanden mellem jorden og objektet)
    rhat = r./norm(r); % retning (enhedsvektor)
    F_g = -G .* m_sat .* m_jorden ./(norm(r).^2) .* rhat; % (N) Kraft
    
    p_sat = p_sat + dt .* F_g ; % Impulsopdatering
    r_sat = r_sat + dt .* p_sat ./ m_sat; % positionsopdatering
    
    plot3(r_sat(1), r_sat(2), r_sat(3),'.') % Plot den opdaterede position
    t = t + dt; % adder tidsskridt
    
    pause(0.1); % vent 100 ms...
end

n-legeme simulering

clf;clear;clc;

% Protostjerne
stjerne = struct();
stjerne.mass = 10e31;
stjerne.pos = [0, 0, 0];
stjerne.vel = [0, 0, 0];
stjerne.p = stjerne.mass .* stjerne.vel;


% Tyngdekraft mellem to masser 
G = 6.67e-11;
F_g = @(s1, s2) - G .* s1.mass .* s2.mass /(norm(s1.pos-s2.pos).^2) .*  (s1.pos-s2.pos)/norm(s1.pos-s2.pos);

N = 3; % Antal stjerner

stjerner(1:N) = struct(stjerne);

% Tilfældige startbetingelser for stjernerne
for i=1:N
    stjerner(i).pos = [rand * 10e10, rand * 10e10, rand* 10e10];
    stjerner(i).vel = [rand * 5000, rand*5000, rand*500];
    stjerner(i).p = stjerne.mass .* stjerne.vel;
    
    plot3(stjerner(i).pos(1), stjerner(i).pos(2), stjerner(i).pos(3))
    hold on
end


dt = 500; % (s) tidsskridt
t = 0; % (s) tid
t_end = 10000000;
% Kør simulering

while t < t_end
    
    F = zeros(N,3);% (N) den resulterende kraft
    for i=1:N % Løb over alle stjerner
        for j=1:N % Løb over alle stjerner 
            if i ~= j % ingen selv-interaktion
                 F(i,:) = F(i) + F_g(stjerner(i), stjerner(j)); 
            end
        end
        
        stjerner(i).p = stjerner(i).p + F(i,:) .* dt;
        stjerner(i).vel = stjerner(i).p ./ stjerner(i).mass;
    end
    
    % opdater position
    for i=1:N
        stjerner(i).pos = stjerner(i).pos + stjerner(i).vel .* dt;
        if mod(t, 10000) == 0
            plot3(stjerner(i).pos(1), stjerner(i).pos(2), stjerner(i).pos(3),'Color',[i/(4+N) i/N i/N],'Marker','.'); 
        end
    end
    pause(0.01)
    
    
    t = t+dt;
end