Построение модели внутреннего устройства Солнца

28 минут назад

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

Обозначим через p(r) суммарное давление на удалении r от центра, m(r) — массу вещества внутри сферы радиуса r, ρ(r) — локальную плотность, T(r) — температурные показатели и L(r) — интегральную светимость на заданном расстоянии. Сформулируем систему из четырех фундаментальных дифференциальных уравнений, определяющих состояние звездного вещества:

1) Уравнение гидростатического баланса (отражает равновесие между градиентом давления и силами гравитационного притяжения):

$\frac{dp}{dr} = -\frac{Gm(r)\rho(r)}{r^2}$

. Здесь G — гравитационная постоянная.

2) Уравнение неразрывности массы (определяет связь между приращением массы и плотностью среды):

$\frac{dm}{dr} = 4\pi r^2\rho(r)$

3) Уравнение энергетического баланса (описывает изменение мощности излучения по радиусу):

$\frac{dL}{dr} = 4\pi r^2(\varepsilon - \varepsilon')$

, где ε обозначает темп генерации энергии на единицу массы, а ε’ — потери, обусловленные нейтринным излучением.

4) Уравнение переноса энергии: При доминировании лучистого механизма

$\frac{dT}{dr} = -\frac{3K\rho L(r)}{64\pi r^2\sigma T^3}$

, где K — коэффициент непрозрачности вещества, а σ — постоянная Стефана-Больцмана.

В случае конвективного теплопереноса уравнение принимает вид:

$\frac{dT}{dr} = -\nabla_{ad}\frac{T}{P}\frac{dp}{dr}$

(источник: Википедия: Строение звёзд).

Уравнение состояния вещества:

$p = p_{газ} + p_{изл} = \frac{\rho kT}{\mu m_u} + \frac{4\sigma}{3c}T^4$

, где k — постоянная Больцмана, m_u — атомная единица массы, c — скорость света.

Для вырожденных объектов (белых карликов):

$p \sim \rho^{5/3} \quad \text{или} \quad p \sim \rho^{4/3}$

Энерговыделение при протон-протонном цикле:

$\varepsilon \sim \rho T^4$

, для массивных звезд с более высокой температурой:

$\varepsilon \sim \rho T^{19/3}$

, интенсивность нейтринных потерь:

$\varepsilon' \sim \rho T^6$

Коэффициент K остается практически неизменным при экстремальных температурах, однако в умеренном диапазоне он подчиняется закону:

$K \sim \rho T^{-3.5}$

(источник: Лекции по общей астрофизике, К.А. Постнов)

Критерий устойчивости Шварцшильда: если локальный температурный градиент ниже адиабатического, энергия переносится излучением; в противном случае развивается конвекция. Как правило, в горячих внутренних областях преобладает лучистый перенос, а во внешних слоях с более низкой температурой — конвективный.

Граничные условия задачи: m(0)=0, L(0)=0 (в центре); p(R)=0, T(R)=T_eff, m(R)=M, L(R)=L_star (на поверхности радиуса R).

На основе вышеизложенных принципов построим математическую имитацию солнечной структуры. Ниже представлен код на Python для численного моделирования Солнца:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
Фундаментальные константы
G = 6.67430e-11      # м^3 кг^-1 с^-2
sigma = 5.670374e-8  # Вт м^-2 К^-4
k_B = 1.380649e-23   # Дж/К
c = 299792458        # м/с
m_u = 1.660539e-27  # кг
Параметры Солнца
M_sun = 1.989e30    # кг
R_sun = 6.9634e8    # м
L_sun = 3.828e26    # Вт
class RealisticSolarModel:
def init(self):
Состав Солнца (массовые доли)
   self.X = 0.735  # водород
   self.Y = 0.249  # гелий
   self.Z = 0.016  # тяжёлые элементы

   # Границы зон (по современным данным)
   self.r_core = 0.25 * R_sun  # ядро до 0.25 R⊙
   self.r_radiative = 0.71 * R_sun  # лучистая зона до 0.71 R⊙

   # Средняя молекулярная масса
   self.mu = 4 / (5 * self.X + 3 * self.Y + 1)
def equation_of_state(self, rho, T):

"""Упрощённое уравнение состояния идеального газа"""

P_gas = rho  k_B  T / (self.mu  m_u)

P_rad = (4  sigma / (3  c))  T**4

return P_gas + P_rad
def opacity(self, rho, T):

"""Упрощённая модель непрозрачности"""
Электронное рассеяние
   kappa_es = 0.2 * (1 + self.X)

   # Формула Крамерса
   kappa_ff = 4.3e22 * self.Z * (1 + self.X) * rho / T**3.5

   # H⁻ поглощение (упрощённо)
   kappa_Hminus = 2.5e-32 * self.Z * rho**0.5 * T**9

   # Гармоническое среднее
   if kappa_es > 0 and kappa_ff > 0 and kappa_Hminus > 0:
       kappa = 1 / (1/kappa_es + 1/kappa_ff + 1/kappa_Hminus)
   else:
       kappa = max(kappa_es, kappa_ff, kappa_Hminus)

   return max(kappa, 1e-6)  # ограничение снизу
def energy_generation(self, rho, T):

"""Энерговыделение с учётом pp-цепочки"""

eps_pp = 1.07e-7 * self.X2  (rho / 1e5)  (T / 1e6)4
CNO цикл (упрощённо)
   eps_CNO = 8.68e-12 * self.X * self.Z * (rho / 1e5) * (T / 1e6)**(12)
   # потери на нейтрино
   eps_neyt=1.0e-13*self.X**2 *(rho / 1e5) * (T / 1e6)**6
   return eps_pp + eps_CNO-eps_neyt
def temperature_gradient(self, r, M_r, L_r, P, T, rho):

"""Градиент температуры с критерием Шварцшильда"""
Лучистый градиент
   try:
       nabla_rad = 3 * self.opacity(rho, T) * rho * L_r * P / (64 * np.pi * sigma * G * M_r * T**4)
   except:
       nabla_rad = 0

   # Адиабатический градиент
   nabla_ad = 0.4

   # Критерий Шварцшильда
   if nabla_rad > nabla_ad and r > self.r_core:
       # Конвекция
       dP_dr = -G * M_r * rho / r**2
       return nabla_ad * (T / P) * dP_dr
   else:
       # Лучистый перенос
       return -3 * self.opacity(rho, T) * rho * L_r / (16 * np.pi * sigma * T**3 * r**2)
def equations(self, r, y):

"""Система дифференциальных уравнений"""

M_r, P, L_r, T = y
   # Плотность из уравнения состояния
   try:
       rho = (P - (4 * sigma / (3 * c)) * T**4) * self.mu * m_u / (k_B * T)
       rho = max(rho, 1e-8)  # защита от отрицательных значений
   except:
       rho = 1e-8

   # Градиенты
   dM_dr = 4 * np.pi * r**2 * rho
   dP_dr = -G * M_r * rho / r**2

   epsilon = self.energy_generation(rho, T)
   dL_dr = 4 * np.pi * r**2 * epsilon

   dT_dr = self.temperature_gradient(r, M_r, L_r, P, T, rho)

   return [dM_dr, dP_dr, dL_dr, dT_dr]
def solve_structure(self, r_min=1e6, r_max=R_sun, n_points=500):

"""Решение системы уравнений"""
Начальные условия в центре Солнца
   T_center = 1.57e7  # К
   rho_center = 1.52e5  # кг/м³
   P_center = self.equation_of_state(rho_center, T_center)
   L_center = 0  # светимость в центре

   y0 = [0, P_center, L_center, T_center]

   # Радиальная сетка
   r_points = np.logspace(np.log10(r_min), np.log10(r_max), n_points)

   # Решение системы ОДУ
   solution = solve_ivp(
       self.equations,
       [r_points[0], r_points[-1]],
       y0,
       t_eval=r_points,
       method='BDF',
       rtol=1e-5,
       atol=1e-7
   )

   if not solution.success:
       raise RuntimeError(f"Решение не сошлось: {solution.message}")

   return solution.t, solution.y
Запуск моделирования
print("Начинаем реалистичное моделирование структуры Солнца...")

model = RealisticSolarModel()
try:

r, (M_r, P, L_r, T) = model.solve_structure()

print("Моделирование завершено успешно!")

except Exception as e:

print(f"Ошибка при моделировании: {e}")

raise
Расчёт плотности
rho = np.zeros_like(r)

for i, (ri, Pi, Ti) in enumerate(zip(r, P, T)):

try:

rho[i] = (Pi - (4  sigma / (3  c)) * Ti*4)  model.mu  m_u / (k_B  Ti)

rho[i] = max(rho[i], 1e-8)

except:

rho[i] = 1e-8
Визуализация результатов
fig, axs = plt.subplots(2, 2, figsize=(14, 10))
График 1: распределение массы
axs[0, 0].plot(r / R_sun, M_r / M_sun, 'b-', linewidth=2, label="Модель")

axs[0, 0].set_xlabel('Радиус (R⊙)')

axs[0, 0].set_ylabel('Масса (M⊙)')

axs[0, 0].set_title('Распределение массы внутри Солнца')

axs[0, 0].grid(True, alpha=0.3)

axs[0, 0].axvline(x=0.25, color="r", linestyle="--", alpha=0.7, label="Граница ядра")

axs[0, 0].axvline(x=0.71, color="g", linestyle="--", alpha=0.7, label="Граница лучистой зоны")

axs[0, 0].legend()
График 2: распределение давления
axs[0, 1].plot(r / R_sun, P, 'r-', linewidth=2)

axs[0, 1].set_xlabel('Радиус (R⊙)')

axs[0, 1].set_ylabel('Давление (Па)')

axs[0, 1].set_title('Распределение давления')

axs[0, 1].set_yscale('log')

axs[0, 1].grid(True, alpha=0.3)

axs[0, 1].axvline(x=0.25, color="r", linestyle="--", alpha=0.7)

axs[0, 1].axvline(x=0.71, color="g", linestyle="--", alpha=0.7)
График 3: распределение светимости
axs[1, 0].plot(r / R_sun, L_r / L_sun, 'g-', linewidth=2)

axs[1, 0].set_xlabel('Радиус (R⊙)')

axs[1, 0].set_ylabel('Светимость (L⊙)')

axs[1, 0].set_title('Распределение светимости')

axs[1, 0].grid(True, alpha=0.3)

axs[1, 0].axvline(x=0.25, color="r", linestyle="--", alpha=0.7)

axs[1, 0].axvline(x=0.71, color="g", linestyle="--", alpha=0.7)
График 4: распределение температуры
axs[1, 1].plot(r / R_sun, T, 'm-', linewidth=2)

axs[1, 1].set_xlabel('Радиус (R⊙)')

axs[1, 1].set_ylabel('Температура (К)')

axs[1, 1].set_title('Распределение температуры')

axs[1, 1].set_yscale('log')

axs[1, 1].grid(True, alpha=0.3)

axs[1, 1].axvline(x=0.25, color="r", linestyle="--", alpha=0.7)

axs[1, 1].axvline(x=0.71, color="g", linestyle="--", alpha=0.7)
plt.tight_layout()

plt.show()
Распределение плотности
plt.figure(figsize=(10, 6))

plt.plot(r / R_sun, rho, 'orange', linewidth=2)

plt.xlabel('Радиус (R⊙)')

plt.ylabel('Плотность (кг/м³)')

plt.title('Распределение плотности внутри Солнца')

plt.grid(True, alpha=0.3)

plt.yscale('log')

plt.axvline(x=0.25, color="r", linestyle="--", alpha=0.7, label="Ядро")

plt.axvline(x=0.71, color="g", linestyle="--", alpha=0.7, label="Лучистая зона")

plt.legend()

plt.show()
Аналитика
print("\n" + "="60)

print("АНАЛИЗ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ")

print("="60)
print(f"\nКЛЮЧЕВЫЕ ПАРАМЕТРЫ В ЦЕНТРЕ:")

print(f"Температура: {T[0]:.2e} К (стандарт: 1.57×10⁷ К)")

print(f"Плотность: {rho[0]:.2e} кг/м³ (стандарт: 1.52×10⁵ кг/м³)")

print(f"Давление: {P[0]:.2e} Па")
print(f"\nПАРАМЕТРЫ НА ПОВЕРХНОСТИ (R⊙):")

print(f"Температура поверхности: {T[-1]:.0f} К (наблюдаемая: 5778 К)")

print(f"Светимость на поверхности: {L_r[-1]/L_sun:.4f} L⊙")
Определение конвективной границы
nabla_T = np.abs(np.diff(T) / np.diff(r))

conv_boundary_idx = np.argmax(nabla_T[len(r)//3:]) + len(r)//3

conv_boundary = r[conv_boundary_idx] / R_sun
print(f"\nГРАНИЦЫ ЗОН:")

print(f"Ядро: до {model.r_core/R_sun:.2f} R⊙")

print(f"Лучистая зона: {model.r_core/R_sun:.2f}–{conv_boundary:.2f} R⊙")

print(f"Конвективная зона: от {conv_boundary:.2f} R⊙")


   Визуализированные результаты симуляции приведены далее:
Динамика изменения массы, давления, светимости и температуры в зависимости от радиуса
Изменение плотности вещества от центра к периферии
        Графики наглядно демонстрируют, что около половины солнечной массы и практически весь объем генерируемой энергии сосредоточены непосредственно в ядре. Температурная кривая характеризуется постепенным снижением с последующим резким падением до фотосферных значений (порядка 6000 К). Аналогичным образом ведут себя показатели давления и плотности, стабилизируясь у поверхности.
    В ходе исследования на основе численного моделирования была воссоздана структура Солнца и получены профили распределения ключевых астрофизических параметров по его радиусу.
Анализ полученных данных подтверждает адекватность модели.
Представленный математический аппарат и программный код являются универсальными и могут быть применены для изучения широкого класса звездных объектов.
Реализованное решение затрагивает фундаментальные аспекты современной астрофизики и открывает возможности для дальнейших вычислений.




 
Источник