import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
class Params:
def __init__(self):
self.g = 9.8
self.R = 287
self.m_initial = 3400
self.m_final = 3000
self.h_mass_change = 6000
self.parachutes = {
"引导伞": {"A": 24, "Cd": 0.8, "收口比": 1.0, "收口时间": 0},
"减速伞": {"A": 200, "Cd": 1.0, "收口比": 0.35, "收口时间": 8},
"主伞": {"A": 1200, "Cd": 1.2, "收口比": 0.15, "收口时间": 8}
}
self.seasonal_offsets = {"summer": 0.02}
params = Params()
def get_actual_pressure(h: float, season: str) -> float:
if h < 11000:
P_standard = 101325 * (1 - h / 44331) ** 5.256
else:
P_standard = 22632 * np.exp((11000 - h) / 6340)
random_factor = 1 + np.random.uniform(-0.05, 0.05)
return P_standard * random_factor * (1 + params.seasonal_offsets.get(season, 0.0))
def calculate_temperature(h: float) -> float:
T_standard = 288.15 - 0.0065 * h
random_perturb = np.random.normal(0, 2.0)
return T_standard + random_perturb
def dynamics(t: float, y: np.ndarray, season: str) -> np.ndarray:
v, h = y
if h >= 8000:
phase = "引导伞"
elif 7000 <= h < 8000:
phase = "减速伞"
else:
phase = "主伞"
para = params.parachutes[phase]
t_since_deploy = t if phase == "引导伞" else t - para["收口时间"]
if t_since_deploy < para["收口时间"]:
A = para["A"] * para["收口比"]
Cd = para["Cd"] * para["收口比"]
else:
A = para["A"]
Cd = para["Cd"]
m = params.m_initial if h >= params.h_mass_change else params.m_final
P = get_actual_pressure(h, season)
T = calculate_temperature(h)
rho = P / (params.R * T)
dvdt = params.g - (0.5 * rho * Cd * A * v**2) / (2 * m)
dhdt = -v
return np.array([dvdt, dhdt])
def simulate_single(season: str):
# 定义关键事件检测函数
events = []
def event_8000(t, y): return y[1] - 8000 # 引导伞切换减速伞
event_8000.terminal = False
event_8000.direction = -1
def event_7000(t, y): return y[1] - 7000 # 减速伞切换主伞
event_7000.terminal = False
event_7000.direction = -1
def event_6000(t, y): return y[1] - params.h_mass_change # 质量变化
event_6000.terminal = True
event_6000.direction = -1
# 第一阶段积分（质量变化前）
sol1 = solve_ivp(
lambda t, y: dynamics(t, y, season),
[0, 200],
[200, 10000],
events=[event_8000, event_7000, event_6000],
method='RK45',
t_eval=np.linspace(0, 200, 1000)
)
# 记录关键事件
event_data = []
for i, event_list in enumerate(sol1.t_events):
if event_list.size > 0:
t = event_list[0]
y = sol1.sol(t)
if i == 0:
event_data.append(('引导伞切换', t, y[0], y[1]))
elif i == 1:
event_data.append(('减速伞切换', t, y[0], y[1]))
elif i == 2:
event_data.append(('质量变化', t, y[0], y[1]))
# 第二阶段积分（质量变化后）
if sol1.t_events[2].size > 0:
y_event = sol1.y_events[2][0]
sol2 = solve_ivp(
lambda t, y: dynamics(t, y, season),
[sol1.t_events[2][0], 400],
y_event,
method='RK45',
t_eval=np.linspace(sol1.t_events[2][0], 400, 1000))
# 合并数据
t_total = np.concatenate([sol1.t, sol2.t[1:]])
v_total = np.concatenate([sol1.y[0], sol2.y[0][1:]])
h_total = np.concatenate([sol1.y[1], sol2.y[1][1:]])
else:
t_total = sol1.t
v_total = sol1.y[0]
h_total = sol1.y[1]
# 记录收口结束时间
for phase in ["减速伞", "主伞"]:
if phase == "减速伞":
start_event = next((e for e in event_data if e[0] == '减速伞切换'), None)
else:
start_event = next((e for e in event_data if e[0] == '主伞切换'), None)
if start_event:
start_time = start_event[1]
para = params.parachutes[phase]
end_time = start_time + para["收口时间"]
if end_time <= t_total[-1]:
# 使用插值获取精确值
v = np.interp(end_time, t_total, v_total)
h = np.interp(end_time, t_total, h_total)
event_data.append((f"{phase}收口结束", end_time, v, h))
return {'time': t_total, 'velocity': v_total, 'altitude': h_total, 'events': event_data}
# 获取单次模拟数据
np.random.seed(42) # 固定随机种子
data = simulate_single("summer")
# 绘制双轴图表
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
# 速度-时间曲线
ax1.plot(data['time'], data['velocity'], color='royalblue', lw=2)
ax1.set_xlabel('时间 (s)', fontsize=12)
ax1.set_ylabel('速度 (m/s)', fontsize=12)
ax1.set_title('返回舱速度随时间变化', fontsize=14)
# 速度-高度曲线
ax2.plot(data['velocity'], data['altitude'], color='darkorange', lw=2)
ax2.set_xlabel('速度 (m/s)', fontsize=12)
ax2.set_ylabel('高度 (m)', fontsize=12)
ax2.set_title('返回舱速度随高度变化', fontsize=14)
# 标注关键事件
annotation_params = {
'xytext': (15, -15),
'textcoords': 'offset points',
'arrowprops': dict(arrowstyle="->", color='gray'),
'bbox': dict(boxstyle="round", alpha=0.9, fc="w"),
'fontsize': 9
}
for ax, x_key, y_key in zip([ax1, ax2], ['time', 'velocity'], ['velocity', 'altitude']):
for event in data['events']:
name, t, v, h = event
x_val = t if x_key == 'time' else v
y_val = v if y_key == 'velocity' else h
ax.scatter(x_val, y_val, color='crimson', zorder=10)
ax.annotate(
f"{name}\n时间: {t:.1f}s\n速度: {v:.1f}m/s\n高度: {h:.0f}m",
(x_val, y_val),
**annotation_params
)
# 设置公共属性
for ax in [ax1, ax2]:
ax.grid(True, alpha=0.3)
ax.set_axisbelow(True)
plt.tight_layout()
plt.show()更正上述代码，写出完整正确代码
最新发布
本文讨论了不同编程语言中参数传递的方式，尤其是按值传递和按引用传递的区别。通过实验发现某些类型如array和record似乎是按引用传递的，但在特定方法调用后行为发生了变化。
扫一扫
1732
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
