Моделирование строения Солнца / Хабр

Рассмотрим задачу о строении звёзд: Примем сферически-симметричную квазистатическую модель строения звезды (звезда это огромный шар, все параметры симметричны относительно центра звезды, находятся в равновесии друг с другом), никаких турбулентностей не происходит.

Пусть p(r) - полное давление на расстоянии r от центра, m(r)- масса , заключённая в шаре радиуса r, ρ(r)- плотность, T(r)- температура, L(r)- светимость на расстоянии r от центра. Запишем 4 основных дифференциальных уравнения, описывающих состояние звезды:

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}$

(источник https://ru.wikipedia.org/wiki/Строение_звёзд ).

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

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

, где R- газовая постоянная, m- молярная масса, 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}$

(источник https://dl.libcats.org/genesis/44000/208da625d87ceef3ed66f01bcbae383a/_as/[Postnov_K.A.]_Lekcii_po_obshei_astrofizike_dlya_f(libcats.org).pdf )

Критерий Шварцшильда: Если фактический температурный градиент ниже адиабатического, то перенос лучистый, иначе конвективный ( источник https://kpfu.ru/portal/docs/F_1009337479/Belyaeva.E.E..Fizika.zvyozd.pdf ). В целом, при высоких температурах перенос лучистый, а при

низких- конвективный.

Граничные условия: m(0)=0, L(0)=0, p(R) =0,T( R)=Teff, m( R)=M,

L( R)=Lз, где 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.X**2 * (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⊙ (точно: 1.000 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⊙ (стандарт: 0.25–0.71 R⊙)")
print(f"Конвективная зона: от {conv_boundary:.2f} R⊙ до поверхности")

Результат моделирования представлен ниже в виде графиков:

Распределение массы, плотности, светимости и температуры по радиусу

Из графиков видно, что почти половина массы Солнца и почти вся его светимость находятся в ядре, температура сначала плавно, а потом резко падает до температуры поверхности(6000 К), оставаясь постоянной в фото сфере. Давление и плотность также уменьшается, затем стабилизируясь на некотором постоянном значении.

Таким образом, в данной статье при помощи моделирования исследования структура Солнца, построены распределения по радиус основных его параметров.

Проведён неглубокий анализ полученных графиков.

Уравнения и структура кода примененимы к любой звезде, не обязательно к Солнцу.

Решена очень интересная и важная задача в области астрофизики.