Gibt es eine mathematische Möglichkeit, die Dichte, das Volumen usw. eines Planeten zu berechnen, indem man einfach die Masse als Daten verwendet?

Gliederung des Algorithmus

Wie L. Dutch sagte, können Sie dies tun, wenn Sie das Volumen oder die durchschnittliche Dichte des Planeten annehmen. Aber das Volumen und die durchschnittliche Dichte sind keine Grundwerte – sie selbst hängen auf eine Weise vom Material und der Zusammensetzung ab, die ich gleich beschreiben werde. Als ein Wort der Warnung – ich gehe davon aus, dass Sie Kalkül dafür verstehen – wenn Sie dies nicht tun, ist die Antwort von L. Dutch so gut, wie Sie es tun können.

Die grundlegende Schwierigkeit dieses Problems besteht also darin, dass Sie die Dichteverteilung kennen müssen, um die Masse des Planeten zu kennen. Die andere Antwort geht davon aus, dass es konstant ist, aber das ist keine fundierte Annahme - auf das Innere des Planeten drückt viel mehr Druck, was zu höheren Dichten führt. Wie gehen wir also mit dieser Schwierigkeit um? Natürlich mit Differentialgleichungen! Abgesehen von irgendwelchen lustigen Dingen wie einem schnell rotierenden Planeten oder einem, der sich noch ansammelt, ist die Hauptgleichung, die Sie verwenden würden, die hydrostatische Gleichgewichtsgleichung, die im Grunde besagt, dass der Druck, der gegen eine dünne Materialhülle drückt, genau die Schwerkraft aufheben muss, die sie nach unten zieht :

Wo$r$ist der radiale Abstand vom Mittelpunkt des Planeten,$P(r)$ist der Druck in dieser Tiefe,$\rho(r)$die Dichte und$G$die Gravitationskonstante.$m(r)$ist eine Funktion, die Ihnen die Masse innerhalb der Radiuskugel gibt$r$, und als solche$m(r) = \int_0^r 4 \pi r^2 \rho(r') dr'$. Aber wir wollen ein System von Differentialgleichungen lösen, also differenzieren wir dieses um zu bekommen

Jetzt wollen wir nur dieses Gleichungssystem unter Einbeziehung Ihrer Randbedingung integrieren. Moment mal, wir sind noch nicht ganz fertig – es gibt drei Funktionen, die wir bestimmen wollen:$\rho$,$P$, Und$m$, aber nur zwei Gleichungen. Wir brauchen eine andere Gleichung, um diese Sache zu lösen!

Diese Gleichung hat die Form der sogenannten Zustandsgleichung, wo die Zusammensetzung des Planeten ins Spiel kommt. Es stellt sich heraus, dass es für ein bestimmtes Material eine Funktion gibt, die die Dichte und den Druck in Beziehung setzt, die Zustandsgleichung genannt wird (die wir werde durch bezeichnen$g$):

Ich sollte zunächst erwähnen, dass das EOS auch eine andere Variable mit diesen beiden in Beziehung setzt – normalerweise Temperatur, Entropie oder spezifische innere Energie. Für die Zwecke der Bestimmung von Dichteprofilen sind die Effekte jedoch ziemlich gering (~ wenige %), sodass Sie die Temperatur normalerweise einfach auf einen vernünftigen Wert wie 5.000 K einstellen und g entlang dieser Isotherme ziehen können.

Jetzt kann es hier etwas knifflig werden – die Zustandsgleichung führt nicht garantiert zu einer Eins-zu-eins-Beziehung zwischen Dichte und Druck. Mit anderen Worten, wir können nicht unbedingt schreiben$P = f(\rho)$Und$\rho = f^{-1}(P)$. Wenn dies der Fall ist, müssen wir einen Löser für eine sogenannte Differential-Algebra-Gleichung verwenden . Sie können dieses Verhalten mit einem normalen ODE-Paket nachahmen, das in den meisten Sprachen verfügbar ist, aber es könnte ziemlich langsam sein. Eine weitere Option ist die Sprache Julia, die über ein großartiges Differentialgleichungspaket verfügt, das algebraische Differentialgleichungen lösen kann . Normalerweise können wir für dieses Regime jedoch zumindest schreiben$P = f(\rho)$obwohl wir es nicht umkehren können, was für dieses spezielle Gleichungssystem gut genug ist, um es umzuschreiben als:

Wie finden wir nun einen EOS für ein bestimmtes Material? Nun, es gibt mehrere Möglichkeiten:

1. Am gebräuchlichsten auf dem Gebiet der Planetenwissenschaft ist die Verwendung einer vorgefertigten Tabelle, die von jemandem erstellt wurde, der mehr über EOS weiß als Sie. SESAME und FPEOS sind zwei, die ich kenne - letzteres ist leicht zugänglich, enthält aber leider nicht viele wichtige Elemente (vor allem Eisen, das für felsige Planeten wichtig ist). Ersteres ist vollständiger, aber Sie müssen den Zugriff über den von mir angegebenen Link anfordern. In jedem Fall müssen Sie bei dieser Option einige Zeit damit verbringen, andere Skripte zu lernen oder Ihre eigenen Skripte zu schreiben, um die Tabellen zu analysieren.
2. Sie können in der Literatur herumstöbern, um eine Arbeit mit einer Formel oder einem Plot zu finden, die Sie digitalisieren können , wie diese hier .
3. Sie können einfach eine Art empirische Formel angeben, die Ihrer Meinung nach vernünftig genug erscheint. Machtrechtliche Beziehungen wie$P = C\rho^{\gamma}$sind dafür oft beliebt.

Puh, wir haben endlich unsere Gleichungen aufgestellt! Jetzt müssen wir nur noch unsere Randbedingungen verwenden: if$R$der Radius eures Planeten ist, diese Bedingung ist das$m(R) = M$, Wo$M$ist die Masse, die der Planet haben soll. Aber leider, da wir nicht wissen, was$R$Sollte dies der Fall sein, müssen wir einen iterativen Ansatz wählen, der durch den folgenden Kontrollfluss beschrieben wird:

1. Satz$m(r = \epsilon) = 0$Und$\rho(r = \epsilon) = \rho_{guess}$. Hier$\epsilon$ist eine kleine positive Zahl, aber nicht ganz$0$weil unsere Gleichungen dort singulär sind. Kann wahrscheinlich eine Variablenänderung vornehmen, um dies zu beheben, aber es ist physikalisch weniger klar, was los ist, und ich bin auch faul.
2. Lösen Sie die beiden obigen gekoppelten Gleichungen, bis Sie den Radius erhalten$R$so dass$P(R) = f(\rho(R)) = 0$(da an der Oberfläche des Planeten kein Druck herrscht). Möglicherweise muss aus numerischen Gründen ein kleiner positiver Druck als Null verwendet werden.
3. Überprüfen$M(R)$.

• Wenn es nahe genug an der gewünschten Masse liegt, sind Sie fertig! Extrahieren Sie die Dichte- und Druckprofile aus der ODE-Lösung.
• Wenn es zu groß ist, machen$\rho_{guess}$kleiner und wiederholen Sie Schritt 2
• Wenn es zu klein ist, machen$\rho_{guess}$größer und wiederholen Sie Schritt 2

Und das ist alles, was dazu gehört! Wenn ich später Zeit habe, werde ich versuchen, ein einfaches Beispielskript hinzuzufügen, aber hoffentlich reicht das aus, um Ihnen den Einstieg zu erleichtern.

EDIT: Bilder und ein einfacher Code

Also habe ich einen einfachen Python-Code erstellt, um Planetenmassen zu berechnen, das Ergebnis ist unten. Es sollte ziemlich einfach zu bedienen sein, alles, was Sie tun müssen, ist, die entsprechenden Pakete installiert zu haben und Mdesired und die Zustandsgleichungsfunktion nach Ihren Wünschen zu ändern:

import scipy.interpolate as interp
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import scipy.integrate as sciint
from math import pi

G = 6.67408e-11   # gravitational constant in [m^3][kg^-1][s^-2]
rho_guess = 1e4   # initial guess for center of planet density in [kg][m^-3]
Mdesired = 6e24          # Desired planet mass in [kg]

# Now we normalize units so that G, M, rho_guess=1.
# This will help numerical stability.
rhostar = rho_guess
Mstar = Mdesired
Tstar = (G * rhostar)**(-1/2)
Lstar = (Mstar / rhostar)**(1/3)
Pstar = Mstar / (Lstar * Tstar**2)
Mdesired = Mdesired / Mstar
rho_guess = rho_guess / rhostar
M0 = 0.0
eps = 1e-4  # start ODE at this value of r to avoid singularity

def dudr(r, u):
    # u[0] = rho, u[1] = m
    drhodr = -(u[1]/(r**2 * dPdrho(u[0])) * u[0])
    dmdr = 4*pi*r**2*u[0]
    return [drhodr, dmdr]

def P(rho):
    """EOS of epsilon iron phase from OSTI 6345571 
    """
    rho0 = 8.430*1e3 / rhostar
    beta0 = 182*1e9 / Pstar
    betaprime0 = 5.0
    eta = rho/rho0

    P = 1.5*beta0*((eta**(7/3) - eta**(5/3))
                   *(1 + (3/4)*(eta**(2/3) - 1)*(betaprime0 - 4)))
    return P

rho_EOS_arr = np.linspace(0, 50, 1000)
dPdrho_arr = np.gradient(P(rho_EOS_arr), rho_EOS_arr)
dPdrho = interp.interp1d(rho_EOS_arr, dPdrho_arr, bounds_error=False, 
                         fill_value=(0.0, 0.0))

def found_surface(r, u):
    """Event that terminates ODE when surface is reached
    """
    Psurf = 1e-7 
    return P(u[0]) - Psurf
found_surface.terminal = True

def M(rho_guess, plot=False, obj = False):
    """Find mass of planet given guess for density at center.
    Can also plot density profiles and give ODE solution object
    """
    sol = sciint.solve_ivp(dudr, (eps, 100), (rho_guess, M0),
                           events = found_surface, 
                           t_eval = np.linspace(eps, 50, 10000))

    if plot:
        plt.figure()
        plt.plot(sol.t*Lstar/1e3, sol.y[0,:]*rhostar/1e3)
        plt.title(r"$\rho$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$\rho$ $(g/cm^3)$")
        
        plt.figure()
        plt.plot(sol.t*Lstar/1e3, sol.y[1,:]*Mstar/1e3)
        plt.title(r"$M$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$M$ $(kg)$")
    
    if obj:    
        return sol
    else:
        return sol.y[1,-1]

rhoc = opt.newton(lambda rho: M(rho) - Mdesired, rho_guess, tol = 1e-3, maxiter=200)
                      
sol = M(rhoc, plot=True, obj=True)
print(f"radius is {sol.t[-1]*Lstar/1e3} km")
print(f"mass is {sol.y[1,-1]*Mstar:e} kg")

plt.show()

Im Moment habe ich es so eingerichtet, dass es Mdesireddie Masse der Erde ist, um es zu messen. Hier ist, was es für das Dichteprofil der Erde vorhersagt . Es ist nicht sehr falsch – es sagt voraus, dass der Erdradius ~5000 km statt ~6400 km beträgt. Es hat auch ein etwas ähnliches Dichteprofil wie der Kern.

Der Hauptgrund für die Diskrepanz ist, dass diese Berechnung davon ausgeht, dass die gesamte Erde aus Eisen besteht$\epsilon$-- dies gilt wirklich nur für den Kern und beschreibt daher die Kruste und den Mantel als zu dicht und das Endergebnis ist, dass der Radius der Erde unterschätzt wird. Um dies zu beheben, müssten Sie die Zustandsgleichung so ändern, dass sich ihr Wert mit ändert$m(r)$. Dies ist eine gute Übung und wenn Sie das überhaupt verstanden haben, empfehle ich, es auszuprobieren!

EDIT 2: Ein genauerer Code

Also bin ich wirklich hingegangen und habe mich auf diesem Nerd selbst geschnippelt . Unten ist ein von mir entwickelter Code, von dem ich denke, dass er Ihnen im Prinzip die genauen Dichteprofile geben kann, die Sie in der Literatur sehen (ohne Temperaturkorrekturen), solange Sie über geeignete Zustandsgleichungen verfügen. Mir ist aufgefallen, dass ich in meiner letzten Bearbeitung einen Fehler gemacht habe, als ich den Weg nach vorne für mehrere EOSs beschrieben habe - wenn Sie die EOS nur abrupt variieren lassen$m(r)$an materiellen Grenzen wird es zu unphysikalischen Ergebnissen kommen. Der Grund dafür ist, dass der Code implizit ein kontinuierliches Dichteprofil annimmt, während in Wirklichkeit das Druckprofil kontinuierlich ist.

Es gibt ein paar Möglichkeiten, damit umzugehen, aber was mein Code tut, ist, dass er Material für Material für das Dichteprofil auflöst und an den Schnittstellen mit den richtigen Randbedingungen zusammenfügt. Als Eingaben benötigen Sie:

1. Eine gewünschte Masse für den Planeten in$kg$
2. Eine Schätzung der Dichte im Zentrum des Planeten in$kg/m^3$. Dies muss nicht sehr genau sein - ~ 10000 wird wahrscheinlich funktionieren, aber wenn Sie mit sehr hohen oder niedrigen Massen arbeiten, müssen Sie möglicherweise Anpassungen vornehmen, damit die Solver konvergieren
3. Eine Sammlung von Materialien, die die Zusammensetzung Ihres Planeten definieren. Es gibt zwei Listen, die diesem Punkt entsprechen: Eine gibt den Anteil der gesamten Planetenmasse für jedes Material an, und eine ist eine Liste von Zustandsgleichungen für jedes Material. Ich habe solche für Eisen und Silikat-Perowskit beigefügt, was eine gute Annäherung an die Erde ist, aber wenn Sie andere Materialien wollen, lesen Sie einfach die Literatur und finden Sie weitere Formeln.

Und das ist es! Alles, was Sie ändern müssen, sind die Variablen im Eingabebereich, alles andere sollte die Schwerarbeit leisten. Ohne weitere Umschweife hier der Code:

import scipy.interpolate as interp
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import scipy.integrate as sciint
from math import pi

#####################################
###### INPUTS #######################
#####################################

def P_silicates(rho):
    """EOS of cold silicate perovskite from Bina 1995
    """
    Ktau0 = 262*1e9 / Pstar
    Kprimetau0 = 4
    rho0 = 4.1 * 1e3 / rhostar
    f = 0.5*((rho/rho0)**(2/3) - 1)
    xi = -(3/4)*(Kprimetau0 - 4)
    Pc = 3*Ktau0*f*(1+2*f)**(5/2)*(1-xi*f)
   
    return Pc

def P_iron(rho):
    """EOS of epsilon iron phase from OSTI 6345571 
    """
    rho0 = 8.430*1e3 / rhostar
    beta0 = 182*1e9 / Pstar
    betaprime0 = 5.0
    eta = rho/rho0
    
    Pc = 1.5*beta0*((eta**(7/3) - eta**(5/3))
                    *(1 + (3/4)*(eta**(2/3) - 1)*(betaprime0 - 4)))
    return Pc

G = 6.67408e-11   # gravitational constant in [m^3][kg^-1][s^-2]
rho0_guess = 1e4   # initial guess for center of planet density in [kg][m^-3]
Mdesired = 6e24          # Desired planet mass in [kg]
mat_Mfracs = [0.3, 0.7]  # Mass fraction for each material type, starting from
                         # core outward
mat_EOSs = [P_iron, P_silicates]  # EOS for each material type


##############################################
######## SOLVER ##############################
##############################################

# Now we normalize units so that G, M, rho_guess=1.
# This will help numerical stability.
msum = sum(mat_Mfracs)
mat_Mfracs = [val/msum for val in mat_Mfracs]
rhostar = rho0_guess
Mstar = Mdesired
Tstar = (G * rhostar)**(-1/2)
Lstar = (Mstar / rhostar)**(1/3)
Pstar = Mstar / (Lstar * Tstar**2)
Mdesired = Mdesired / Mstar
rho0_guess = rho0_guess / rhostar
M0 = 0.0

def dudr(r, u, Mdesired, mat_props):
    # u[0] = rho, u[1] = m
    drhodr = -(u[1]/(r**2 * dPdrho(*u, Mdesired, mat_props)) * u[0])
    dmdr = 4*pi*r**2*u[0]
    return [drhodr, dmdr]

def dPdrho(rho, m, Mdesired, mat_props):
    Pfunc = mat_props["Pfunc"]
    drho = 1e-9
    P1 = Pfunc(rho+drho)
    Pn1 = Pfunc(rho-drho)
    return (P1 - Pn1) / (2*drho)

def end_of_material(r, u, Mdesired, mat_props):
    """Event that terminates ODE when the end of a shell for a given material 
    is reached.
    """
    cum_mass_frac_end = mat_props["cum_mass_frac_end"]
    return u[1]/Mdesired - cum_mass_frac_end
end_of_material.terminal = True    

def found_surface(r, u, Mdesired, mat_props):
    """Event that terminates ODE when surface is reached
    """
    Psurf = 1e-7 
    Pfunc = mat_props["Pfunc"]
    return Pfunc(u[0]) - Psurf
found_surface.terminal = True

def solve_one_mat(rho0, m0, Mdesired, mat_props):
    r0 = mat_props["r0"]
    last_mat = mat_props["last_mat"]
    if last_mat:
        sol = sciint.solve_ivp(
            dudr, (r0, 100), (rho0, m0),
            args = (Mdesired, mat_props),
            events = found_surface, 
            t_eval = np.linspace(r0, 100, 100000),
            tol = 1e-8
            )
    else:
        sol = sciint.solve_ivp(
            dudr, (r0, 100), (rho0, m0),
            args = (Mdesired, mat_props),
            events = (found_surface, end_of_material), 
            t_eval = np.linspace(r0, 100, 100000),
            tol = 1e-8
            )
    return sol

def M(rho0, Mdesired, mat_Mfracs, mat_EOSs, plot=False):
    """Find mass of planet given guess for density at center.
    Can also plot density profiles and give ODE solution object
    """
    cumsum_Mfracs = np.cumsum(mat_Mfracs)
    rho0, m0, r0 = [rho0, 0.0, 1e-7]  # solver starts at 1e-7 radius to avoid
                                      # singularity
    P0 = mat_EOSs[0](rho0)
    r_arr = np.array([])
    rho_arr = np.array([])
    m_arr = np.array([])
    
    for i in range(len(mat_Mfracs)):
        # Build the mat_props dict for a given material, this tells the
        # solver what material we're using
        mat_props = {}
        mat_props["Pfunc"] = mat_EOSs[i]
        mat_props["cum_mass_frac_end"] = cumsum_Mfracs[i]
        mat_props["r0"] = r0
        mat_props["last_mat"] = (i == len(mat_EOSs)-1)
        rho0 = opt.newton(lambda rho: mat_EOSs[i](rho) - P0, rho0, tol = 1e-8)
        
        sol = solve_one_mat(rho0, m0, Mdesired, mat_props)        
        rho0, m0, r0 = [sol.y[0,-1], sol.y[1,-1], sol.t[-1]]
        P0 = mat_EOSs[i](rho0)
        r_arr = np.concatenate((r_arr, sol.t))
        rho_arr = np.concatenate((rho_arr, sol.y[0,:]))
        m_arr = np.concatenate((m_arr, sol.y[1,:]))
        
        # If surface is reached before we can get through all materials,
        # we must end the loop
        if len(sol.t_events[0]) != 0 and not(mat_props["last_mat"]):
            missed_layers = True
            break
        else: 
            missed_layers = False
    
    if plot:
        plt.figure()
        plt.plot(r_arr*Lstar/1e3, rho_arr*rhostar/1e3)
        plt.title(r"$\rho$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$\rho$ $(g/cm^3)$")
        
        plt.figure()
        plt.plot(r_arr*Lstar/1e3, m_arr*Mstar)
        plt.title(r"$M$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$M$ $(kg)$")
    
    return (m_arr[-1], r_arr, rho_arr, m_arr, missed_layers)

# Solve for density at planet center that gives the desired planetary mass
rhoc = opt.newton(
    lambda rho0: M(rho0, Mdesired, mat_Mfracs, mat_EOSs)[0] - Mdesired, 
    rho0_guess, tol = 1e-4, maxiter = 200
    )
# Solve for and plot final values                          
R, r_arr, rho_arr, m_arr, missed_layers = M(
    rhoc, Mdesired, mat_Mfracs, mat_EOSs, plot=True
    )
print(f"radius is {r_arr[-1]*Lstar/1e3} km")
print(f"mass is {m_arr[-1]*Mstar:e} kg")

Im Moment habe ich einige Eingabeparameter, die für die Erde ungefähr korrekt sind, nämlich$M_{desired} = 6\times 10^{24} kg$, mit einer Zusammensetzung von 30 % Eisen und 70 % Silikaten. Wie funktioniert es also? Eigentlich ganz gut! Hier ist das Dichteprofil, das es ausgibt:

Das kommt den Profiltypen ziemlich nahe, die Sie sehen, wenn Sie eine schnelle Google-Suche durchführen! Es schätzt auch den Erdradius auf 6250 km, verdammt nahe am tatsächlichen Wert von 6370 km. Hoffentlich war dies nützlich für Sie!

Ein paar Vereinfachungen, die ich vorgenommen habe, wenn Sie mehr vertiefen möchten

1. Ich ignoriere die Temperaturabhängigkeit. In Wirklichkeit wird das EOS 3 Variablen in Beziehung setzen und Sie benötigen eine dritte Gleichung, die Ihnen sagt, wie sich die Temperatur mit der Tiefe ändert. Dies ist eine kniffligere Gleichung, da sie von der Zeit von Akkretion und radioaktivem Zerfall und einer ganzen Reihe anderer Dinge abhängt, die Sie angeben müssen.
2. In Wirklichkeit wird die Form des EOS nicht das EOS für ein einzelnes Material sein – wenn Sie tiefer gehen, ändert sich die Zusammensetzung des Planeten. Ich bin mir nicht sicher, wie Sie dies nach den ersten Prinzipien berechnen - es gibt wahrscheinlich eine Art thermodynamisches / Auftriebsargument, aber ich bin mir nicht sicher. BEARBEITEN - mein zweiter Code befasst sich damit, wenn es klare Grenzen zwischen Materialtypen gibt. In Wirklichkeit wird es eine Übergangszone von einiger Dicke geben, aber es scheint, dass die Literatur dies sowieso ignoriert.
3. Wenn es auf Ihrem Planeten andere Kräfte wie Zentripetalkräfte durch Drehung gibt, kann dies die Druckgleichgewichtsgleichung ändern.
4. Bei superschweren Objekten wie Neutronensternen müssen Sie GR über die Tolman-Oppenheimer-Volkoff-Gleichung berücksichtigen

Um die Dichte zu berechnen$\rho$oder die Lautstärke$V$eines Körpers können Sie die Beziehung verwenden$\rho=M/V$.

Da es sich um zwei unbekannte Werte und eine einzige Gleichung handelt, können Sie sie nicht lösen, es sei denn, Sie haben eine andere Einschränkung, die Sie einbringen können. Wenn Sie beispielsweise den Radius angeben können, können Sie das Volumen berechnen als$V=4/3 \pi r^3$und daraus die Dichte.

Eine andere mögliche Extrapolation ist, dass Sie, wenn Sie davon ausgehen, dass die durchschnittliche Dichte der Erde oder einem anderen Planeten entspricht, den Sie in Betracht ziehen möchten, aus der Dichte das Volumen berechnen können.