Pole grawitacyjne Ziemi i symulacja spadającej książki

    W poprzednim poście o spadającej do środka Ziemi książce pozostawiłem bez wyjaśnienia przebieg siły grawitacji w zależności od odległości od środka Ziemi. Tutaj pokażę skąd to wziąłem. Ponadto pokażę, co będzie gdy uwzględnimy zmienną gęstość Ziemi i ile czasu wtedy będzie spadać nasza książka. \(\newcommand\ddfrac[2]{\frac{\displaystyle #1}{\displaystyle #2}}\)

    Zacznijmy od umieszczenia Ziemi w układzie współrzędnych. Posłużymy się poniższym rysunkiem:

    Ziemię - zewnętrzna powierzchnia - umieszczamy tak, że jej środek pokrywa się ze środkiem układu współrzędnych. Wewnętrzna powierzchnia reprezentuje tę część Ziemi, która zawiera się w kuli o promieniu $r$. Aby wyznaczyć pole grawitacyjne $\mathbf{g}$ skorzystamy z prawa Gaussa. Prawo Gaussa dla pola grawitacyjnego można wyprowadzić z prawa powszechnego ciążenia Newtona. Najpierw chciałem to tutaj zrobić ale stwierdziłem, że ciekawsze będzie pokazanie symulacji numerycznej wrzucania książki do Ziemi. Zatem, pomijając wyprowadzenie, prawo Gaussa wygląda tak:

\begin{equation}
\oint_S\ \mathbf{g}\mathbf{\hat{n}}dS=-4\pi GM
\label{gauss}
\end{equation}

    Lewa strona równania to strumień pola grawitacyjnego przechodzącego przez dowolną zamkniętą powierzchnę $S$. Po prawej stronie $G$ oznacza stałą grawitacji, a $M$ masę zawartą wewnątrz rozpatrywanej powierzchni $S$. Skoro $S$ może być dowolne wybierzmy sferę o środku w środku Ziemi i promieniu $r$, czyli wewnętrzną sferę z obrazka. Zakładamy, że część Ziemi zawarta wewnątrz $S$ jest sferycznie symetrycza, czyli przechodzi na samą siebie przy obrocie o dowolny kąt wokół dowolnej osi przechodzącej przez środek sfery i przy odbiciu przez dowolną płaszczyznę zawierającą środek sfery. Skoro tak, to pole grawitacyjne generewane przez tę część Ziemi też musi być sferycznie symetryczne. Oznacza to, że na powierzchni $S$ przyjmuje postać

\begin{equation}
\mathbf{g}(\mathbf{r})=g(r)\mathbf{\hat{n}}
\end{equation}

    Jest to pole prostopadłe do powierzchni $S$ w każdym jej punkcie - $\mathbf{\hat{n}}$ to wektor jednostkowy prostopadły do sfery - i mające wartość $g(r)$. Podstawmy tę postać do równania \eqref{gauss}.

\begin{equation}
\oint_S g(r)\mathbf{\hat{n}}\mathbf{\hat{n}}dS=g(r)\oint_S dS=4\pi r^2g(r)
\end{equation}\begin{equation}4\pi r^2g(r)=-4\pi GM\end{equation}\begin{equation}g(r)=-\frac{GM}{r^2}\label{g_field}
\end{equation}

     Pozostaje nam tylko poznać masę $M$ zawartą w kuli o promieniu $r$. W pierwszym przybliżeniu założyliśmy, że gęstość Ziemi jest jednorodna i wynosi $\rho_\oplus$. Wówczas mamy:

\begin{equation}
M(r)=\frac{4}{3}\pi \min(r, R_\oplus)^3 \rho_\oplus
\end{equation}

    Podstawiając tę postać masy do \eqref{g_field} odtworzymy wzór z poprzedniego artykułu. Mamy bowiem dla $r<R_\oplus$:

\begin{equation}
g(r)=-\frac{G}{r^2}\frac{4}{3}\pi r^3 \rho_\oplus=-\frac{G}{R_\oplus^3}(\frac{4}{3}\pi R_\oplus^3 \rho_\oplus) r=-\frac{GM_\oplus}{R_\oplus^3}r
\end{equation}
oraz dla $r\geq R_\oplus$:
\begin{equation}
g(r)=-\frac{G}{r^2}\frac{4}{3}\pi R_\oplus^3 \rho_\oplus=-\frac{GM_\oplus}{r^2}
\end{equation}

    Ciekawie się robi, gdy gęstość Ziemi nie jest jednorodna i jest funkcją położenia. Oznaczmy gęstość Ziemi $\rho(\mathbf{r})=\rho (r)$. Przebieg gęstości weźmiemy z badań sejsmologicznych [1]. Wtedy, aby poznać masę Ziemi w zależności od odległości od środka, a co za tym idzie i wartość pola grawitacyjnego musimy odcałkować tę gęstość.

\begin{equation}
M(r)=\int_V \rho(\mathbf{r'})d^3\mathbf{r'}=\int_0^{2\pi}d\phi \int_0^{\pi}\sin{\theta}d\theta \int_0^r r'^2\rho(r')dr'=4\pi\int_0^r r'^2\rho(r')dr'
\label{mass_int}
\end{equation}

    Na poniższych rysunkach przedstawiony jest przebieg gęstości Ziemi oraz wartość pola grawitacyjnego odcałkowany numerycznie z wykorzystaniem wzorów \eqref{g_field} i \eqref{mass_int}.

Gęstość Ziemi na podstawie pracy [1] - model PREM; kliknij, aby powiększyć.

Wartość pola grawitacyjnego Ziemi w funkcji odległości od jej środka; kliknij, aby powiększyć.

    Jak możemy zauważyć przebieg wartości pola grawitacyjnego wewnątrz Ziemi jest bardzo ciekawy. Idąc od powierzchni przez blisko 2000 km jest prawie stałe - lekko rośnie a następnie lekko maleje. Dla $r=3480 \ \mathrm{km}$ osiąga maksimum równe $g_{max}=10,689\ \ddfrac{\mathrm{m}}{\mathrm{s}^2}$, a następnie monotonicznie maleje w kierunku środka Ziemi. Aby obliczyć tor ruchu książki wrzuconej do Ziemi musimy rozwiązać następujące równanie różniczkowe, wynikające z drugiej zasady dynamiki Newtona:

\begin{equation}
\frac{d^2r}{dt^2}=g(r)
\label{newton}
\end{equation}

    W modelu PREM [1] gęstość Ziemi jest wyrażona przez zestaw wielomianów stopnia co najwyżej trzeciego. Oznacza to, że pole grawitacyjne Ziemi będzie wyrażone przez zestaw wielomianów stopnia co najwyżej czwartego. Rozwiązanie \eqref{newton} wymaga zatem rozwiązania wysoce nieliniowych równań różniczkowych. Nawet nie chcę sprawdzać czy to jest wykonalne, więc posłużę się metodą numeryczną zwaną metodą różnic skończonych. Polega ona na przybliżeniu pochodnej jej definicją tylko bez granicy. Najbardziej stabilne rozwiązania otrzymamy przy zastosowaniu symetrycznej definicji, tzn.:

\begin{equation}
\frac{df}{dx}\approx \frac{f(x+\ddfrac{1}{2}h)-f(x-\ddfrac{1}{2}h)}{h}\end{equation}\begin{equation}\frac{d^2f}{dx^2}\approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}\label{fin_diff}
\end{equation}
    Stosując \eqref{fin_diff} do \eqref{newton} otrzymujemy:
\begin{equation}
\frac{r(t+\delta t)-2r(t)+r(t-\delta t)}{\delta t^2}=g(r(t))\end{equation}\begin{equation}r(t+\delta t)=2r(t)-r(t-\delta t)+g(r(t))\delta t^2 \label{rule}
\end{equation}
Potrzebujemy jeszcze warunków początkowych. Dla naszego problemu mamy:
\begin{equation}
r(0)=R_\oplus, r(-\delta t)=R_\oplus
\end{equation}

    Mamy już wszystko czego potrzeba, żeby przeprowadzić symulację. Korzystamy z \eqref{rule} i warunków początkowych. Na wykresie zamieszczam wykres toru ruchu książki spadającej do Ziemi. Różowa krzywa to symulacja, niebieska to rozwiązanie ścisłe z poprzedniego artykułu.

Kliknij, aby powiększyć

    Widzimy zatem, że książka dotrze do środka Ziemi jeszcze szybciej, bo już w 19 minut i 5 sekund. Następnie będzie lecieć dalej, tym razem zwalniając aż dotrze na drugi koniec Ziemi, a wtedy zawróci. Ten ruch okresowy zaobserwujemy też, gdy przeprowadzimy dalej symulację. Oto rysunek:

Kliknij, aby powiększyć

    Policzyliśmy. Teraz wiemy ile czasu zajmie książce dotarcie do środka Ziemi. Inna sprawa, że pewnie nigdy nie uda nam się wydrążyć takiego tunelu :). Dla zainteresowanych zamieszczam kod w R użyty do przeprowadzenia symulacji i stworzenia rysunków: kod

Bibliografia:
[1] Dziewonski, Adam M., and Don L. Anderson. "Preliminary reference Earth model." Physics of the earth and planetary interiors 25.4 (1981): 297-356. Page 308