模拟太阳系。
在本文中，我们将看看如何创建一个完整的太阳系模型，不仅展示了引力相互作用，而且还允许我们研究系统关键参数变化的后果，例如中心恒星的质量。 所提出的实现扩展了天体力学的经典问题，包括热力学方面，这使得研究太阳质量变化对行星轨道动力学和气候条件的复杂影响成为可能。
实施的关键方面:
-
牛顿万有引力定律的实现为一个由9个相互作用的物体（太阳+8个行星）组成的系统，每个物体影响所有其他物体，而不仅仅是中心物体。
-
该模型包括恒星的基本质量光度比（L~M3*⁵·，这使我们能够估计太阳质量的变化如何影响其能量释放。 基于这种关系，行星的平衡温度以摄氏度计算。
-
已经创建了一个单独的计算和可视化系统，显示行星的温度如何随太阳质量的变化而变化，范围从1到10太阳质量。 结果以带有温度标记的图形和表格的形式呈现（温度高于1000°C，温度低于0°c）。
-
使用欧拉方法求解具有自适应间距的运动方程，在长达10年的时间尺度上建模时提供精度和性能之间的平衡。
这个项目是学习Engee的一个很好的起点，展示了编程语言如何处理需要计算效率、清晰的代码和可视化的任务。 在以下部分中，我们将详细研究代码的架构，物理基础，温度计算的实现和该模型的教育方面，以及考虑如何在教育过程和科普项目中使用此类模拟。
现在让我们继续分析代码块，并从工作中使用的常量开始。
·引力常数（G=6.67430×10-11m3·s⁻2kg⁻1）**是一个基本的物理常数，它决定了有质量的物体之间的引力。 在经典牛顿力学中，它包含在万有引力定律中，根据该定律，两个物体之间的吸引力与它们的质量乘积成正比，与它们中心之间距离的平方成反比。 正是这个常数将时空的几何形状与物质的分布联系起来，决定了非相对论极限中引力相互作用的强度。 在我们的模拟中，引力常数允许我们计算太阳吸引行星的加速度，反之亦然，形成逼真的轨道动力学。
**天文单位（AU=1.496×1011米）**是天文学中历史上建立的距离测量单位，大约等于从地球到太阳的平均距离。 将此常数引入模型有两个重要目的。 首先，它提供了一个自然尺度来表示行星的轨道，因为它们的轨道的半长轴传统上以天文单位表示。 其次，使用AU简化了模型参数的感知和调整，因为研究人员习惯于在天文背景下使用这些值，而不是以绝对公制单位进行操作。
**斯特凡-玻尔兹曼常数（σ=5.670374419×10-∞Wm-2K⁻∞）**是描述黑体辐射的斯特凡-玻尔兹曼定律的基本物理常数。 根据这一规律，单位时间内黑体表面的一个单位辐射的总能量与其热力学温度的四次方成正比。 在我们的模型背景下，这个常数是计算行星平衡温度的关键元素。 它允许您从太阳的光度（通过质量确定）切换到到达行星的热流，并进一步切换到其表面的计算温度，同时考虑到反照率。 正是由于这个常数，模型获得了一个热力学组件，将引力动力学与行星上的气候条件联系起来。
const G = 6.67430e-11 # 引力常数
const AU = 1.496e11 # 天文单位(m)
const YEAR = 365.25 * 86400.0 # 以秒为单位的年份
const SOLAR_MASS = 1.989e30 # 太阳质量(公斤)
const SOLAR_RADIUS = 6.957e8 # 太阳半径(m)
const STEFAN_BOLTZMANN = 5.670374419e-8
结构 Body 它是模型的一个基本元素，封装了天体的完整物理和视觉轮廓。 由于是可变的，它允许您在模拟的每一步动态更新对象的状态，反映它们随时间的演变。
它包括三类参数:
*动态特性（位置 x, y，速度 vx, vy，加速度 ax, ay)
*物理性质（质量 mass，物理和轨道半径）
*可视化和分析的属性（名称，颜色，温度）。
场尤为重要 temperature，这是根据太阳的光度和到身体的距离来计算的。
mutable struct Body
name::String
mass::Float64 # 重量(公斤)
radius::Float64 # 物理半径(m)
orbit_radius::Float64 # 初始轨道半径(m)
color::Symbol
x::Float64 # 位置X(m)
y::Float64 # Y位置(m)
vx::Float64 # 速度X(m/s)
vy::Float64 # 速度Y(m/s)
ax::Float64 # 加速度X(m/s2)
ay::Float64 # 加速度Y(m/s2)
temperature::Float64 # 平衡温度(°C)
end
阵列 planets_data 它包含了太阳系八大行星的初始天体物理参数和视觉属性，作为创建该类型物体的初始条件 Body 在模拟中。
planets_data = [
("Mercury", 3.285e23, 2.4397e6, 0.39*AU, :orange),
("Venus", 4.867e24, 6.0518e6, 0.72*AU, :green),
("Earth", 5.972e24, 6.371e6, 1.0*AU, :blue),
("Mars", 6.39e23, 3.3895e6, 1.52*AU, :red),
("Jupiter", 1.898e27, 6.9911e7, 5.2*AU, :brown),
("Saturn", 5.683e26, 5.8232e7, 9.58*AU, :gold),
("Uranus", 8.681e25, 2.5362e7, 19.22*AU, :cyan),
("Neptune", 1.024e26, 2.4622e7, 30.05*AU, :darkblue)
]
功能 solar_luminosity 实现了天体物理学定律"质量-光度"（L〇M3·〇），将太阳质量的相对变化转换为其能量释放的适当乘数。 功能 equilibrium_temperature_c 基于这种光度和到行星的距离，它使用能量吸收和再发射模型计算其平衡温度（以摄氏度为单位），同时考虑反照率和斯特凡-玻尔兹曼定律。
function solar_luminosity(mass_factor)
return mass_factor^3.5 # L ~ M^3.5
end
function equilibrium_temperature_c(luminosity_factor, distance_au)
# 简化的温度公式（摄氏度）
albedo = 0.3
solar_constant = 1361.0 # 土地用W/m2
# 开尔文温度
temp_k = (solar_constant * luminosity_factor * (1 - albedo) /
(4 * STEFAN_BOLTZMANN))^0.25
# 距离调整和转换为摄氏
return (temp_k / sqrt(distance_au)) - 273.15
end
功能 create_temperature_stats 进行全面的热力学分析，计算各种太阳质量（从1到10太阳质量）的所有行星的平衡温度，创建两种类型的图表—对恒星质量和距离的温度依赖性，显示结果表，并将两个图表保存到文件中以备将来使用，作为可视化气候情景的自主模块。.
function create_temperature_stats()
sun_mass_factors = [1, 10, 100, 1000, 1e4, 1e5, 1e6, 1e7]
planets_info = [
("Mercury", 0.39, :orange),
("Venus", 0.72, :green),
("Earth", 1.0, :blue),
("Mars", 1.52, :red),
("Jupiter", 5.2, :brown),
("Saturn", 9.58, :gold),
("Uranus", 19.22, :cyan),
("Neptune", 30.05, :darkblue)
]
# 温度的计算
temperatures = Dict()
for (name, distance_au, color) in planets_info
temps = Float64[]
for mass_factor in sun_mass_factors
lum_factor = solar_luminosity(mass_factor)
temp_c = equilibrium_temperature_c(lum_factor, distance_au)
push!(temps, temp_c)
end
temperatures[name] = temps
end
p1 = plot(size=(1000, 600),
title="行星的温度对太阳质量的依赖",
xlabel="太阳的质量（以太阳质量计）",
ylabel="温度(°C)",
xscale=:log10,
yscale=:linear,
grid=true, gridstyle=:dash,
legend=:topleft,
framestyle=:box)
for (name, _, color) in planets_info
plot!(p1, sun_mass_factors, temperatures[name],
label=name, color=color, linewidth=2, marker=:circle)
end
hline!(p1, [0], label="水的冰点", color=:blue, linestyle=:dash)
hline!(p1, [100], label="水的沸点", color=:red, linestyle=:dash)
annotate!(p1, 1e5, 1500, text("极端条件", 10, :gray))
p2 = plot(size=(1000, 600),
title="行星的温度取决于距离",
xlabel="离太阳的距离(AU)",
ylabel="温度(°C)",
xscale=:log10,
yscale=:linear,
grid=true,
legend=:topright)
distances = [d for (_, d, _) in planets_info]
planet_names = [n for (n, _, _) in planets_info]
for (i, mass_factor) in enumerate([1, 10, 100, 1000])
temps = [temperatures[name][i] for name in planet_names]
plot!(p2, distances, temps,
label="M☉ = $(Int(mass_factor))×",
linewidth=2,
marker=:circle)
end
display(p1)
display(p2)
println("\n" * "="^70)
println("行星在不同太阳质量下的温度（°C）")
println("="^70)
print("太阳的质量: ")
for mf in sun_mass_factors
print(lpad("$(Int(mf))×", 10))
end
println()
for (name, _, _) in planets_info
print(rpad(name, 12))
for temp in temperatures[name]
if temp > 1000
print(lpad("🔥$(round(Int, temp))", 10))
elseif temp > 100
print(lpad("$(round(Int, temp))", 10))
elseif temp > 0
print(lpad("$(round(Int, temp))", 10))
else
print(lpad("❄$(round(Int, temp))", 10))
end
end
println()
end
println("="^70)
savefig(p1, "temperature_vs_mass.png")
savefig(p2, "temperature_vs_distance.png")
return temperatures
end
功能 calculate_initial_velocity 计算一个天体的初始轨道速度，确保它沿着围绕太阳的圆形轨迹运动，根据给定的身体位置和中心恒星的质量，按照经典的牛顿力学。
function calculate_initial_velocity(sun_mass, body_x, body_y)
r = sqrt(body_x^2 + body_y^2)
v_circ = sqrt(G * sun_mass / r)
return (-body_y / r * v_circ, body_x / r * v_circ)
end
以下是主要功能 simulate_solar_system 它初始化并执行系统的重力动力学和热状态的完整模拟循环。 它从创建一组身体开始 bodies，其中包括由给定系数调整的质量的太阳 sun_mass_factor，以及所有计算出圆形轨道初始速度的行星。 时间步长等积分参数 dt 和总步数 n_steps，将输入的天文值（年，天）转换为计算间隔，为数值解做好准备。
在主循环中，该函数使用牛顿定律重新计算每个步骤中所有物体对的重力加速度，并使用欧拉方法更新它们的速度和坐标。 与此同时，行星的当前平衡温度是根据太阳的亮度和与太阳的距离计算的，这为纯机械模拟增加了热力学背景。 使用指定的频率，系统状态被保存以供后续可视化。
计算完成后，该函数创建一个动画GIF文件，其中帧顺序显示系统的演变：太阳作为中心的简约点，行星根据其温度（从冷蓝色到红热红色）着色，以及几乎不可 动画标题动态地反映了关键参数，包括地球的温度，提供了一个关于改变中心恒星质量的计算实验结果的可视化表示。
function simulate_solar_system(sun_mass_factor=1.0; years=10, step_days=5, filename="solar_system.gif")
sun_mass = SOLAR_MASS * sun_mass_factor
sun_luminosity = solar_luminosity(sun_mass_factor)
bodies = [Body("Sun", sun_mass, SOLAR_RADIUS, 0.0, :yellow,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5778.0 - 273.15)]
for (name, mass, radius, orbit_r, color) in planets_data
vx, vy = calculate_initial_velocity(sun_mass, orbit_r, 0.0)
body = Body(name, mass, radius, orbit_r, color,
orbit_r, 0.0, vx, vy, 0.0, 0.0, 0.0)
push!(bodies, body)
end
dt = step_days * 86400.0
n_steps = Int(round(years * YEAR / dt))
frames_per_gif = 150
step_skip = max(1, n_steps ÷ frames_per_gif)
positions_history = []
function compute_accelerations!(bodies)
for body in bodies
body.ax = 0.0
body.ay = 0.0
end
n = length(bodies)
for i in 1:n
for j in i+1:n
dx = bodies[j].x - bodies[i].x
dy = bodies[j].y - bodies[i].y
r_sq = dx^2 + dy^2
r = sqrt(r_sq)
if r > 1e-6
force = G * bodies[i].mass * bodies[j].mass / r_sq
ax_i = force * dx / (r * bodies[i].mass)
ay_i = force * dy / (r * bodies[i].mass)
ax_j = -force * dx / (r * bodies[j].mass)
ay_j = -force * dy / (r * bodies[j].mass)
bodies[i].ax += ax_i
bodies[i].ay += ay_i
bodies[j].ax += ax_j
bodies[j].ay += ay_j
end
end
end
end
println("模拟：M÷=÷（sun_mass_factor）×，÷（年）年，步÷（step_days）天")
for step in 0:n_steps
if step % step_skip == 0
pos = []
for body in bodies
if body.mass > 0.0
if body.name != "Sun"
dist_au = sqrt(body.x^2 + body.y^2) / AU
body.temperature = equilibrium_temperature_c(sun_luminosity, dist_au)
end
push!(pos, (body.x, body.y, body.color, body.name, body.temperature))
end
end
push!(positions_history, pos)
end
step == n_steps && break
compute_accelerations!(bodies)
for body in bodies
body.vx += body.ax * dt
body.vy += body.ay * dt
body.x += body.vx * dt
body.y += body.vy * dt
end
end
max_dist = 1.2 * maximum([b.orbit_radius for b in bodies if b.name != "Sun"])
anim = @animate for (frame, positions) in enumerate(positions_history)
time_years = frame * step_skip * dt / YEAR
earth_temp = 0.0
for (x, y, color, name, temp) in positions
if name == "Earth"
earth_temp = temp
break
end
end
p = plot(size=(1000, 800), aspect_ratio=:equal,
xlim=(-max_dist, max_dist), ylim=(-max_dist, max_dist),
title="M☉=$(sun_mass_factor)×|L☉=$(round(sun_luminosity，数字=1))×|$(round(time_years，数字=1))年|T♁=$(round(earth_temp，数字=1))°C",
legend=:none, framestyle=:box,
background_color=:black,
titlefontsize=12)
scatter!(p, [0.0], [0.0],
markersize=2, color=:yellow, marker=:circle,
markerstrokewidth=0, markercolor=:yellow, alpha=0.8)
theta = range(0, 2π, length=100)
zone_radius = 0.05 * AU
circle_x = zone_radius .* cos.(theta)
circle_y = zone_radius .* sin.(theta)
plot!(p, circle_x, circle_y, color=:lightgray, linewidth=0.5, alpha=0.15, label="")
for (x, y, color, name, temp) in positions
if name != "Sun"
body = first([b for b in bodies if b.name == name])
size_factor = 5 + 2 * log10(body.mass / 5.972e24)
marker_color = if temp > 1000
:red
elseif temp > 100
:orange
elseif temp > 0
:green
else
:lightblue
end
scatter!(p, [x], [y],
markersize=max(3, size_factor),
color=marker_color, marker=:circle,
markerstrokewidth=1)
end
end
plot!(p, [max_dist*0.7, max_dist*0.7 + AU], [max_dist*0.9, max_dist*0.9],
color=:white, linewidth=2)
annotate!(p, max_dist*0.7 + AU/2, max_dist*0.88,
text("1 AU", 10, :white))
end
gif(anim, filename, fps=20)
println("动画保存：$文件名")
return bodies
end
让我们运行测试并查看结果。
create_temperature_stats()
simulate_solar_system(1.0, years=10, step_days=3, filename="normal_sun.gif");
.gif)
simulate_solar_system(10.0, years=5, step_days=2, filename="massive_sun_10x.gif");
.gif)
simulate_solar_system(10000.0, years=2, step_days=1, filename="extreme_sun_10000x.gif");
.gif)
结论
所进行的测试证明了太阳质量变化的两个关键物理后果。
温度分析（create_temperature_stats()）清楚地表明，光度对质量的依赖性导致行星温度的极端升高。 即使是太阳质量的10倍，内行星的温度也超过了水的沸点，当质量增加到10-10‰太阳时，所有行星都升温到1000℃以上，这说明了恒星结构与系统气候条件之间的根本关系。
轨道动力学的模拟证实，太阳质量的增加加速了行星的轨道运动（根据开普勒定律），并由于中心体重力的支配而改变了它们轨迹的稳定性。 动画可视化这些效果：在正常质量（normal_sun.gif 行星表现出经典的凯普勒轨道，随着它的生长（massive_sun_10x.gif, extreme_sun_10000x.gif 它们的运动变得更快，颜色从冷色调变为红热。