Есть ли математический способ рассчитать плотность, объем и т. д. планеты, просто используя массу в качестве данных?

Схема алогритма

Итак, как сказал Л. Датч, если вы предполагаете объем или среднюю плотность планеты, вы можете это сделать. Но объем и средняя плотность не являются величинами из первых принципов — они сами зависят от материала и состава, что я вскоре опишу. В качестве предупреждения - я собираюсь предположить, что вы понимаете исчисление для этого - если вы этого не сделаете, ответ Л. Датча будет настолько хорош, насколько вы можете это сделать.

Итак, принципиальная трудность этой задачи состоит в том, что для того, чтобы узнать массу планеты, нужно знать распределение плотности. Другой ответ предполагает, что она постоянна, но это не совсем обоснованное предположение — внутри планеты сжимается гораздо большее давление, что приводит к более высокой плотности. Так как же нам справиться с этой трудностью? Конечно же, с дифференциальными уравнениями! За исключением каких-либо забавных вещей, таких как быстро вращающаяся планета или планета, которая все еще срастается, основное уравнение, которое вы должны использовать, — это уравнение гидростатического равновесия, которое в основном говорит, что давление, воздействующее на тонкую оболочку материала, должно точно компенсировать силу гравитации, тянущую ее вниз. :

где$r$- радиальное расстояние от центра планеты,$P(r)$давление на этой глубине,$\rho(r)$плотность и$G$гравитационная постоянная.$m(r)$это функция, которая дает вам массу внутри сферы радиуса$r$, и как таковой$m(r) = \int_0^r 4 \pi r^2 \rho(r') dr'$. Но мы хотим, чтобы система дифференциальных уравнений решалась, поэтому мы дифференцируем ее, чтобы получить

Теперь мы просто хотим интегрировать эту систему уравнений, включающую ваше граничное условие. Подождите, мы еще не совсем готовы - есть три функции, которые мы хотим определить:$\rho$,$P$, и$m$, а только два уравнения. Нам нужно другое уравнение, чтобы решить эту штуку!

Это уравнение имеет форму так называемого уравнения состояния, в которое входит состав планеты. Оказывается, для данного материала существует функция, которая связывает плотность и давление, называемая уравнением состояния (которое мы буду обозначать через$g$):

Прежде всего, я должен упомянуть, что УС связывает с этими двумя переменными еще одну переменную — обычно температуру, энтропию или удельную внутреннюю энергию. Однако для целей определения профилей плотности эффекты довольно малы (~ несколько%), поэтому обычно вы можете просто установить температуру на какое-то разумное значение, например 5000 К, и вытянуть g вдоль этой изотермы.

Здесь все может стать немного сложнее — уравнение состояния не обязательно приводит к прямой связи между плотностью и давлением. Другими словами, мы не можем обязательно написать$P = f(\rho)$и$\rho = f^{-1}(P)$. Если это так, мы должны использовать решатель для так называемого дифференциально-алгебраического уравнения . Вы можете имитировать это поведение, используя обычный пакет ODE, доступный на большинстве языков, но это может быть довольно медленным. Другим вариантом является язык Julia, который имеет большой пакет дифференциальных уравнений, способный решать алгебро-дифференциальные уравнения . Однако обычно для этого режима мы можем хотя бы написать$P = f(\rho)$хотя мы не можем инвертировать его, что для этой конкретной системы уравнений достаточно хорошо, чтобы переписать их как:

Теперь, как нам найти EOS для данного материала? Ну есть несколько способов:

1. Наиболее распространенным в области планетологии является использование заранее подготовленной таблицы, созданной кем-то, кто знает об EOS больше, чем вы. SESAME и FPEOS — два, о которых я знаю — последний легко доступен, но, к сожалению, не содержит много важных элементов (в первую очередь железа, которое важно для каменистых планет). Первый более полный, но вы должны запросить доступ по ссылке, которую я дал. В любом случае с этой опцией вам придется потратить некоторое время на изучение других сценариев или на написание собственных сценариев для анализа таблиц.
2. Вы можете копаться в литературе, чтобы найти документ с формулой или графиком, который вы можете оцифровать , как этот .
3. Вы можете просто дать какую-нибудь эмпирическую формулу, которая, по вашему мнению, кажется достаточно разумной. Отношения по степенному закону, такие как$P = C\rho^{\gamma}$часто популярны для этого.

Уф, наконец-то мы настроили наши уравнения! Теперь нам просто нужно использовать наши граничные условия: если$R$это радиус вашей планеты, это условие заключается в том, что$m(R) = M$, где$M$- это масса, которую вы хотите, чтобы планета имела. Но, к сожалению, поскольку мы не знаем, что$R$должно быть, мы должны использовать итеративный подход, описанный следующим потоком управления:

1. Набор$m(r = \epsilon) = 0$и$\rho(r = \epsilon) = \rho_{guess}$. Здесь$\epsilon$небольшое положительное число, но не совсем$0$потому что наши уравнения там сингулярны. Вероятно, можно изменить переменную, чтобы исправить это, но физически менее понятно, что происходит, и я также ленив.
2. Решайте два связанных уравнения выше, пока не получите радиус$R$такой, что$P(R) = f(\rho(R)) = 0$(поскольку на поверхности планеты нет давления). Возможно, придется использовать некоторое небольшое положительное давление, отличное от нуля, по числовым причинам.
3. Проверять$M(R)$.

• Если она достаточно близка к желаемой массе, все готово! Извлеките профили плотности и давления из решения ОДУ.
• Если он слишком большой, сделайте$\rho_{guess}$меньше и повторите шаг 2
• Если он слишком мал, сделайте$\rho_{guess}$больше и повторите шаг 2

И это все, что нужно! Если позже у меня будет время, я попытаюсь добавить простой пример сценария, но, надеюсь, этого достаточно, чтобы вы начали.

РЕДАКТИРОВАТЬ: картинки и простой код

Поэтому я сделал простой код на Python для расчета массы планет, результат ниже. Он должен быть довольно прост в использовании, все, что вам нужно сделать, это установить соответствующие пакеты и изменить Mdesired и функцию уравнения состояния по своему вкусу:

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()

Прямо сейчас я настроил его так, чтобы это Mdesiredбыла масса Земли, чтобы сравнить его. Вот что он предсказывает для профиля плотности Земли . Это не так уж и неправильно — он предсказывает, что радиус Земли составляет ~ 5000 км вместо ~ 6400 км. Он также имеет несколько схожий профиль плотности с ядром.

Основная причина несоответствия заключается в том, что в этом расчете предполагается, что вся Земля состоит из железа в$\epsilon$фаза - на самом деле это верно только для ядра, поэтому он описывает кору и мантию как слишком плотные, и конечным результатом является то, что он предсказывает радиус Земли. Чтобы исправить это, вам придется изменить уравнение состояния так, чтобы его значение менялось с$m(r)$. Это хорошее упражнение, и если вы вообще это поняли, я рекомендую попробовать!

РЕДАКТИРОВАТЬ 2: более точный код

Так что я действительно пошел и ботаник подстрелил себя на этом. Ниже приведен код, который я разработал, и я думаю, что в принципе он может дать вам точные профили плотности, которые вы видите в литературе (без температурных поправок), если у вас есть соответствующие уравнения состояния. Я понял, что допустил ошибку в моем последнем редактировании при описании пути вперед для нескольких EOS - если вы просто сделаете EOS резко меняющимся на$m(r)$на материальных границах это даст нефизические результаты. Причина в том, что код неявно предполагает непрерывный профиль плотности, тогда как в действительности непрерывным является профиль давления.

Есть несколько способов справиться с этим, но мой код решает для профиля плотности один материал за раз, сшивая интерфейсы с правильными граничными условиями. В качестве входных данных вам понадобятся:

1. Желаемая масса планеты в$kg$
2. Предположение о плотности в центре планеты в$kg/m^3$. Это не должно быть очень точным - ~ 10000, вероятно, будет работать, но если вы работаете с очень большими или низкими массами, вам может потребоваться настроить, чтобы решатели сходились.
3. Коллекция материалов, определяющих состав вашей планеты. Этому пункту соответствуют два списка: один дает долю общей планетарной массы для каждого материала, а другой представляет собой список уравнений состояния для каждого материала. Я включил формулы для железа и силиката перовскита, что является хорошим приближением для Земли, но если вам нужны другие материалы, просто просмотрите литературу и найдите больше формул.

Вот и все! Все, что вам нужно изменить, это переменные в разделе ввода, все остальное должно выполнять тяжелую работу. Без лишних слов, вот код:

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")

Прямо сейчас у меня есть некоторые входные параметры, которые примерно соответствуют Земле, а именно$M_{desired} = 6\times 10^{24} kg$, с составом 30% железа и 70% силикатов. Итак, как это сделать? На самом деле неплохо! Вот профиль плотности, который он выводит:

Что довольно близко к типам профилей, которые вы увидите, если выполните быстрый поиск в Google! Он также оценивает радиус Земли в 6250 км, что чертовски близко к фактическому значению в 6370 км. Надеюсь, это было полезно для вас!

Несколько упрощений, которые я сделал, если вы хотите копнуть глубже

1. Температурную зависимость игнорирую. На самом деле, EOS будет связывать 3 переменные, и вам понадобится третье уравнение, которое покажет вам, как температура зависит от глубины. Это более сложное уравнение, потому что оно зависит от времени аккреции и радиоактивного распада и целого ряда других вещей, которые вам нужно указать.
2. На самом деле форма ЭОС не будет ЭОС для отдельного материала — по мере того, как вы идете глубже, планетарная композиция меняется. Tbh, я не уверен, как вы вычисляете это из первых принципов - вероятно, есть какой-то аргумент термодинамики / плавучести, но я не уверен. РЕДАКТИРОВАТЬ - мой второй код решает эту проблему, когда между типами материалов есть четкие границы. На самом деле будет переходная зона некоторой толщины, но похоже, что большая часть литературы все равно игнорирует это.
3. Если на вашей планете есть другие силы, такие как центростремительные силы от вращения, это может изменить уравнение баланса давления.
4. Для сверхтяжелых объектов, таких как нейтронные звезды, вы должны учитывать ОТО с помощью уравнения Толмана-Оппенгеймера-Волкова.

Для расчета плотности$\rho$или объем$V$тела можно использовать соотношение$\rho=M/V$.

Будучи двумя неизвестными значениями и одним уравнением, вы не можете решить его, если у вас нет другого ограничения. Если вы можете указать, например, радиус, вы можете рассчитать объем как$V=4/3 \pi r^3$оттуда и плотность.

Другая возможная экстраполяция заключается в том, что если вы предполагаете, что средняя плотность такая же, как у Земли или любой другой планеты, которую вы, возможно, захотите рассмотреть, то по плотности вы можете рассчитать объем.