一种火箭发动机仿真系统深度解读

一、源码整体框架与架构设计1.1 项目结构与模块划分RocketMotor/ ├── src/ │ ├── core/ # 核心计算引擎 │ │ ├── __init__.py │ │ ├── propellant.py # 推进剂类与材料属性 │ │ ├── grain.py # 药柱几何与燃面计算 │ │ ├── chamber.py # 燃烧室模型 │ │ ├── nozzle.py # 喷管模型 │ │ ├── motor.py # 发动机集成类 │ │ ├── solver.py # 数值求解器 │ │ └── physics.py # 物理常数与方程 │ │ │ ├── gui/ # 图形用户界面 │ │ ├── __init__.py │ │ ├── main_window.py # 主窗口控制器 │ │ ├── grain_designer.py # 药柱设计器 │ │ ├── material_editor.py # 材料编辑器 │ │ ├── simulation_panel.py # 仿真控制面板 │ │ ├── results_viewer.py # 结果查看器 │ │ ├── plot_widgets.py # 绘图组件 │ │ └── dialogs.py # 各种对话框 │ │ │ ├── data/ # 数据管理 │ │ ├── __init__.py │ │ ├── database.py # 数据库管理器 │ │ ├── models.py # 数据模型 │ │ ├── exporters.py # 数据导出 │ │ └── importers.py # 数据导入 │ │ │ ├── utils/ # 工具函数 │ │ ├── __init__.py │ │ ├── geometry.py # 几何计算工具 │ │ ├── units.py # 单位转换 │ │ ├── validation.py # 输入验证 │ │ └── logging.py # 日志系统 │ │ │ └── tests/ # 测试套件 │ ├── __init__.py │ ├── test_core.py │ ├── test_gui.py │ └── test_integration.py │ ├── data/ # 静态数据文件 │ ├── propellants.json # 推进剂数据库 │ ├── materials.json # 材料数据库 │ └── examples/ # 示例文件 │ ├── docs/ # 文档 │ ├── api/ # API文档 │ └── user_guide/ # 用户指南 │ ├── examples/ # 示例脚本 ├── requirements.txt # 依赖包 └── setup.py # 安装脚本1.2 架构模式Model-View-Controller (MVC) 变体# 数据流架构 用户输入 → GUI控制器 → 模型更新 → 仿真计算 → 结果返回 → 视图更新二、核心计算模块源码深度分析2.1 推进剂模块 (propellant.py)class Propellant: 推进剂物理属性模型 def __init__(self, name: str, density: float, # 密度 [kg/m³] burn_rate_params: dict, # 燃速参数 thermo_params: dict): # 热力学参数 self.name name self.density density # 燃速模型: r a * P^n self.burn_rate_a burn_rate_params.get(a, 0.003) # 燃速系数 self.burn_rate_n burn_rate_params.get(n, 0.3) # 压力指数 # 热力学属性 self.specific_heat_ratio thermo_params.get(gamma, 1.2) self.molecular_weight thermo_params.get(MW, 0.024) # kg/mol self.chamber_temp thermo_params.get(T0, 2500) # 燃烧温度 K # 计算特征速度 c* sqrt(gamma * R * T0) / (gamma * ((gamma1)/2)^((gamma1)/(2*(gamma-1)))) self._compute_characteristic_velocity() def burn_rate(self, pressure: float) - float: 计算给定压力下的燃速 return self.burn_rate_a * (pressure ** self.burn_rate_n) def _compute_characteristic_velocity(self): 计算特征速度 c* R 8314.462618 / self.molecular_weight # 气体常数 J/(kg·K) gamma self.specific_heat_ratio # 特征速度公式 term gamma * ((gamma 1) / 2) ** ((gamma 1) / (2 * (gamma - 1))) self.c_star math.sqrt(gamma * R * self.chamber_temp) / term2.2 药柱几何模块 (grain.py)class Grain: 药柱几何基类计算燃面推移 def __init__(self, outer_diameter: float, length: float, web_thickness: float, geometry_type: str): self.outer_diameter outer_diameter self.length length self.web web_thickness self.geometry_type geometry_type # 预计算几何参数 self._initial_volume self._compute_volume() self._burn_area_progression self._precompute_burn_area() def _precompute_burn_area(self) - list: 预计算燃面随肉厚的变化 返回: [(肉厚1, 燃面1), (肉厚2, 燃面2), ...] burn_area_data [] n_points 100 for i in range(n_points 1): web i * self.web / n_points area self._compute_burning_area(web) burn_area_data.append((web, area)) return burn_area_data def burning_area(self, web: float) - float: 通过插值获取给定肉厚的燃面 # 查找最近的网格点进行线性插值 idx int(web * len(self._burn_area_progression) / self.web) idx min(idx, len(self._burn_area_progression) - 2) x1, y1 self._burn_area_progression[idx] x2, y2 self._burn_area_progression[idx 1] # 线性插值 if x2 x1: return y1 return y1 (y2 - y1) * (web - x1) / (x2 - x1) class StarGrain(Grain): 星形药柱实现 def __init__(self, outer_diameter: float, length: float, web_thickness: float, num_points: int, point_angle: float, fillet_radius: float): super().__init__(outer_diameter, length, web_thickness, star) self.num_points num_points self.point_angle math.radians(point_angle) self.fillet_radius fillet_radius # 计算星形参数 self._compute_star_parameters() def _compute_burning_area(self, web: float) - float: 计算星形药柱在给定肉厚下的燃面 算法: 基于多边形偏移和布尔运算 # 1. 构建初始星形多边形 points self._generate_star_polygon() # 2. 向内偏移指定肉厚 burned_polygon self._offset_polygon(points, -web) # 3. 计算燃面周长乘以长度 if burned_polygon and len(burned_polygon) 2: perimeter self._polygon_perimeter(burned_polygon) burning_area perimeter * self.length return burning_area return 0.0 def _generate_star_polygon(self) - list: 生成星形多边形顶点 points [] angle_step 2 * math.pi / self.num_points for i in range(self.num_points * 2): angle i * angle_step / 2 # 交替内外半径 if i % 2 0: radius self.outer_diameter / 2 else: # 内半径计算 radius (self.outer_diameter / 2) * math.sin(math.pi / self.num_points) * \ math.tan(self.point_angle / 2) x radius * math.cos(angle) y radius * math.sin(angle) points.append((x, y)) return points2.3 内弹道求解器 (solver.py)class BallisticSolver: 内弹道方程求解器 def __init__(self, method: str RK4, rtol: float 1e-6, atol: float 1e-9): self.method method self.rtol rtol self.atol atol self.results SimulationResults() def solve(self, motor, max_time: float 10.0) - SimulationResults: 求解固体火箭发动机内弹道方程 方程: dP/dt (ṁ_gen - ṁ_nozzle) * R * T0 / V ṁ_gen ρ * A_b * r ṁ_nozzle P * A_t / c* r a * P^n # 初始化状态变量 t 0.0 P motor.initial_pressure web 0.0 mass motor.propellant.density * motor.grain.initial_volume # 存储初始状态 self.results.record(t, P, web, 0.0, mass) # 自适应步长控制 dt 0.001 # 初始步长 min_dt 1e-6 max_dt 0.01 while t max_time and web motor.grain.web: # 计算当前燃面 A_b motor.grain.burning_area(web) if A_b 0: break # 燃烧结束 # 计算燃速 r motor.propellant.burn_rate(P) # 计算质量生成率 m_dot_gen motor.propellant.density * A_b * r # 计算喷管质量流率 m_dot_nozzle P * motor.nozzle.throat_area / motor.propellant.c_star # 计算燃烧室自由容积 V motor.chamber.initial_volume motor.grain.burned_volume(web) # 计算压力导数 dP_dt (m_dot_gen - m_dot_nozzle) * motor.gas_constant * motor.propellant.chamber_temp dP_dt / V # 计算推力 thrust self._compute_thrust(P, motor.nozzle, motor.atmospheric_pressure) # 计算误差并调整步长 if self.method RK4: dt self._adaptive_rk4_step(t, dt, [P, web], motor, dP_dt, r) else: dt self._adaptive_euler_step(t, dt, P, web, motor, dP_dt, r) # 应用约束 dt max(min(dt, max_dt), min_dt) # 更新状态 P dP_dt * dt web r * dt mass - m_dot_gen * dt t dt # 记录结果 self.results.record(t, P, web, thrust, mass) # 检查收敛条件 if P motor.atmospheric_pressure * 0.1: break return self.results def _compute_thrust(self, Pc: float, nozzle, P_amb: float) - float: 计算推力 # 推力系数 Cf self._thrust_coefficient(Pc, nozzle, P_amb) # 推力 推力系数 * 燃烧室压力 * 喉部面积 return Cf * Pc * nozzle.throat_area def _thrust_coefficient(self, Pc: float, nozzle, P_amb: float) - float: 计算推力系数 gamma nozzle.propellant.specific_heat_ratio Pc_over_Pe Pc / nozzle.exit_pressure # 推力系数公式 term1 math.sqrt(2 * gamma**2 / (gamma - 1) * (2 / (gamma 1))**((gamma 1) / (gamma - 1)) * (1 - (nozzle.exit_pressure / Pc)**((gamma - 1) / gamma))) term2 (nozzle.exit_area / nozzle.throat_area) * \ (nozzle.exit_pressure - P_amb) / Pc return term1 term2三、GUI架构与执行流程3.1 主窗口控制器 (main_window.py)class MainWindow(QMainWindow): 主窗口协调所有组件 def __init__(self): super().__init__() # 初始化核心模型 self.motor_model MotorModel() self.simulation_engine BallisticSolver() # 创建GUI组件 self._setup_ui() self._connect_signals() # 加载默认配置 self._load_defaults() def _setup_ui(self): 设置用户界面 # 中心部件使用标签页布局 self.tab_widget QTabWidget() # 创建各个标签页 self.grain_designer_tab GrainDesignerWidget(self.motor_model) self.motor_config_tab MotorConfigWidget(self.motor_model) self.simulation_tab SimulationControlWidget(self.motor_model, self.simulation_engine) self.results_tab ResultsViewerWidget() # 添加标签页 self.tab_widget.addTab(self.grain_designer_tab, 药柱设计) self.tab_widget.addTab(self.motor_config_tab, 发动机配置) self.tab_widget.addTab(self.simulation_tab, 仿真控制) self.tab_widget.addTab(self.results_tab, 结果查看) self.setCentralWidget(self.tab_widget) # 创建菜单栏 self._create_menu() # 状态栏 self.statusBar().showMessage(就绪) def _connect_signals(self): 连接信号与槽 # 仿真完成信号 self.simulation_tab.simulation_completed.connect( self._on_simulation_completed ) # 模型更新信号 self.motor_model.model_updated.connect( self._on_model_updated ) def _on_simulation_completed(self, results): 仿真完成处理 # 更新结果标签页 self.results_tab.update_results(results) # 切换到结果标签页 self.tab_widget.setCurrentWidget(self.results_tab) # 在状态栏显示仿真时间 self.statusBar().showMessage( f仿真完成。总时间: {results.time[-1]:.3f} 秒最大推力: {max(results.thrust):.1f} N )3.2 药柱设计器 (grain_designer.py)class GrainDesignerWidget(QWidget): 药柱图形化设计器 def __init__(self, motor_model): super().__init__() self.motor_model motor_model # 创建绘图区域 self.canvas GrainCanvas() # 创建控制面板 self.control_panel self._create_control_panel() # 布局 layout QHBoxLayout() layout.addWidget(self.canvas, 3) layout.addWidget(self.control_panel, 1) self.setLayout(layout) # 连接信号 self.canvas.grain_modified.connect(self._on_grain_modified) def _create_control_panel(self) - QWidget: 创建控制面板 panel QWidget() layout QVBoxLayout() # 药柱类型选择 self.grain_type_combo QComboBox() self.grain_type_combo.addItems([圆柱形, 星形, 套管形, 车轮形]) self.grain_type_combo.currentTextChanged.connect(self._on_grain_type_changed) # 参数输入 self.param_widgets {} # 通用参数 self.param_widgets[length] self._create_param_input(长度 (mm), 100, 1000, 200) self.param_widgets[diameter] self._create_param_input(直径 (mm), 10, 200, 50) self.param_widgets[web] self._create_param_input(肉厚 (mm), 1, 50, 10) # 添加到布局 layout.addWidget(QLabel(药柱类型:)) layout.addWidget(self.grain_type_combo) for label, widget in self.param_widgets.items(): layout.addWidget(widget) # 添加计算按钮 self.calc_button QPushButton(计算燃面推移) self.calc_button.clicked.connect(self._calculate_burn_area) layout.addWidget(self.calc_button) # 添加燃面推移图 self.burn_area_plot BurnAreaPlot() layout.addWidget(self.burn_area_plot) layout.addStretch() panel.setLayout(layout) return panel四、数据流与执行流程4.1 完整执行流程启动应用 ↓ 加载GUI和默认配置 ↓ 用户设计药柱 → 实时几何计算 → 显示药柱形状 ↓ 用户配置发动机参数 → 验证输入 → 更新模型 ↓ 用户点击仿真按钮 → 触发仿真流程: ├── 1. 收集所有输入参数 ├── 2. 验证参数完整性 ├── 3. 创建仿真任务 ├── 4. 在工作线程中运行仿真 │ ├── 初始化求解器 │ ├── 时间步进循环 │ │ ├── 计算当前燃面 │ │ ├── 计算燃速 │ │ ├── 计算质量流率 │ │ ├── 计算压力导数 │ │ ├── 积分更新状态 │ │ └── 记录数据点 │ └── 返回仿真结果 └── 5. 在主线程中处理结果 ├── 更新结果视图 ├── 绘制曲线 └── 计算性能指标 ↓ 用户分析结果 → 参数调整 → 重新仿真4.2 关键数据流# 数据传递路径 class DataFlow: GUI输入 → 控制器 → 模型更新 → 仿真计算 → 结果处理 → 可视化 def run_simulation(self, inputs: dict) - dict: # 1. 输入验证 if not self._validate_inputs(inputs): raise ValueError(输入参数无效) # 2. 创建推进剂对象 propellant Propellant( nameinputs[propellant_name], densityinputs[density], burn_rate_params{ a: inputs[burn_rate_a], n: inputs[burn_rate_n] }, thermo_params{ gamma: inputs[gamma], T0: inputs[chamber_temp] } ) # 3. 创建药柱对象 grain StarGrain( outer_diameterinputs[grain_diameter], lengthinputs[grain_length], web_thicknessinputs[web_thickness], num_pointsinputs[star_points], point_angleinputs[point_angle] ) # 4. 创建喷管对象 nozzle Nozzle( throat_diameterinputs[throat_diameter], exit_diameterinputs[exit_diameter], efficiencyinputs[nozzle_efficiency] ) # 5. 创建发动机对象 motor SolidMotor( propellantpropellant, graingrain, nozzlenozzle, chamber_pressureinputs[chamber_pressure] ) # 6. 运行仿真 solver BallisticSolver() results solver.solve(motor, max_timeinputs[max_time]) # 7. 后处理 processed_results self._post_process(results) return processed_results五、核心算法实现细节5.1 几何计算优化class GeometryUtils: 几何计算工具优化性能 staticmethod def offset_polygon(polygon: list, offset: float) - list: 多边形偏移算法用于燃面推移计算 使用直线段偏移和圆弧连接的方法 if not polygon or len(polygon) 3: return [] offset_polygon [] n len(polygon) for i in range(n): # 当前点、前一点、后一点 prev polygon[(i - 1) % n] curr polygon[i] next_pt polygon[(i 1) % n] # 计算边向量 v1 (curr[0] - prev[0], curr[1] - prev[1]) v2 (next_pt[0] - curr[0], next_pt[1] - curr[1]) # 计算单位法向量 norm1 GeometryUtils._normalize_vector((-v1[1], v1[0])) norm2 GeometryUtils._normalize_vector((-v2[1], v2[0])) # 计算偏移点 offset_prev ( prev[0] norm1[0] * offset, prev[1] norm1[1] * offset ) offset_curr1 ( curr[0] norm1[0] * offset, curr[1] norm1[1] * offset ) offset_curr2 ( curr[0] norm2[0] * offset, curr[1] norm2[1] * offset ) # 计算两偏移线的交点 intersection GeometryUtils._line_intersection( offset_prev, offset_curr1, offset_curr2, (next_pt[0] norm2[0] * offset, next_pt[1] norm2[1] * offset) ) if intersection: offset_polygon.append(intersection) else: # 平行线情况取中点 offset_polygon.append(( (offset_curr1[0] offset_curr2[0]) / 2, (offset_curr1[1] offset_curr2[1]) / 2 )) return offset_polygon staticmethod def polygon_area(polygon: list) - float: 计算多边形面积鞋带公式 area 0.0 n len(polygon) for i in range(n): j (i 1) % n area polygon[i][0] * polygon[j][1] area - polygon[j][0] * polygon[i][1] return abs(area) / 2.05.2 自适应步长控制算法class AdaptiveStepSolver: 自适应步长控制求解器 def __init__(self, rtol1e-6, atol1e-9, safety0.9, max_factor5.0, min_factor0.2): self.rtol rtol self.atol atol self.safety safety self.max_factor max_factor self.min_factor min_factor def adaptive_rk4_step(self, t, dt, y, func, *args): 自适应RK4步长控制 y: 当前状态向量 [压力, 肉厚] func: 微分方程函数 # 计算两个半步长 k1 func(t, y, *args) k2 func(t dt/2, [y[i] k1[i]*dt/2 for i in range(len(y))], *args) k3 func(t dt/2, [y[i] k2[i]*dt/2 for i in range(len(y))], *args) k4 func(t dt, [y[i] k3[i]*dt for i in range(len(y))], *args) # 计算RK4解 y_rk4 [y[i] (k1[i] 2*k2[i] 2*k3[i] k4[i]) * dt/6 for i in range(len(y))] # 计算两个半步长的RK4用于误差估计 dt_half dt / 2 y_half y for _ in range(2): k1 func(t, y_half, *args) k2 func(t dt_half/2, [y_half[i] k1[i]*dt_half/2 for i in range(len(y_half))], *args) k3 func(t dt_half/2, [y_half[i] k2[i]*dt_half/2 for i in range(len(y_half))], *args) k4 func(t dt_half, [y_half[i] k3[i]*dt_half for i in range(len(y_half))], *args) y_half [y_half[i] (k1[i] 2*k2[i] 2*k3[i] k4[i]) * dt_half/6 for i in range(len(y_half))] t dt_half # 计算误差 error 0.0 for i in range(len(y)): scale self.atol self.rtol * max(abs(y[i]), abs(y_rk4[i])) error ((y_rk4[i] - y_half[i]) / scale) ** 2 error math.sqrt(error / len(y)) # 计算新步长 if error 0.0: new_dt self.max_factor * dt else: new_dt self.safety * dt * (1.0 / error) ** 0.2 new_dt min(self.max_factor * dt, max(self.min_factor * dt, new_dt)) # 如果误差可接受返回RK4解和新步长 if error 1.0: return y_rk4, new_dt else: # 误差太大减小步长重新计算 return None, new_dt六、源码层面的优缺点分析6.1 架构优点模块化设计良好各模块职责单一耦合度低清晰的抽象层次物理模型、数值方法、GUI分离可扩展的接口基类设计允许用户自定义推进剂、药柱形状数据驱动JSON配置文件便于修改和扩展实时反馈GUI与计算引擎通过信号槽机制实时交互6.2 代码质量缺陷数值稳定性问题# 问题代码示例 def compute_pressure_derivative(self, P, A_b, V): # 当V接近0时会导致除零错误 dP_dt (self.m_dot_gen(P, A_b) - self.m_dot_nozzle(P)) * self.R * self.T0 dP_dt / V # 危险V可能为0或接近0 return dP_dt缺乏边界检查def burning_area(self, web): # 没有检查web是否超出药柱厚度 return self._burn_area_data[int(web * 100)] # 可能索引越界硬编码参数class Solver: def __init__(self): self.max_iterations 1000 # 硬编码应作为参数 self.tolerance 1e-6 # 硬编码性能问题# 每次调用都重新计算多边形 def burning_area(self, web): poly self._offset_polygon(self.base_polygon, -web) # 计算量大 return self._polygon_area(poly) # 重复计算错误处理不足def solve(self, motor): try: # 仿真计算 result self._run_simulation(motor) return result except Exception as e: # 过于宽泛的异常捕获 print(f仿真失败: {e}) return None七、重点改进方向源码级实现7.1 数值算法升级class ImprovedSolver: 改进的求解器解决现有问题 def __init__(self): # 使用Scipy的ODE求解器 from scipy.integrate import solve_ivp self.solver solve_ivp self.method RK45 # 自适应Runge-Kutta方法 def solve(self, motor, t_span(0, 10), **kwargs): 使用稳健的ODE求解器 def ode_system(t, y): 定义ODE系统 P, web y # 边界检查 if web motor.grain.web: return [0.0, 0.0] # 燃烧结束 # 计算燃面 A_b motor.grain.burning_area(web) # 计算容积避免除零 V motor.chamber.volume(web) if V 1e-12: # 微小容积处理 V 1e-12 # 计算燃速 r motor.propellant.burn_rate(P) # 质量流率 m_dot_gen motor.propellant.density * A_b * r m_dot_nozzle P * motor.nozzle.At / motor.propellant.c_star # 压力导数 dP_dt (m_dot_gen - m_dot_nozzle) * motor.gas_constant * motor.T0 dP_dt / V # 肉厚导数 dweb_dt r return [dP_dt, dweb_dt] # 事件当燃烧结束时停止 def burn_end(t, y): return y[1] - motor.grain.web burn_end.terminal True # 初始条件 y0 [motor.P_initial, 0.0] # 求解ODE sol self.solver( ode_system, t_span, y0, methodself.method, events[burn_end], rtol1e-6, atol1e-9, max_step0.01, **kwargs ) return sol7.2 性能优化class OptimizedGrain(Grain): 性能优化的药柱基类 def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) # 预计算燃面表并使用插值 self._init_burn_area_interpolator() # 预计算容积表 self._init_volume_interpolator() def _init_burn_area_interpolator(self): 初始化燃面插值器 # 生成燃面-肉厚表 webs np.linspace(0, self.web, 1000) areas np.array([self._compute_burning_area_raw(w) for w in webs]) # 创建三次样条插值 from scipy.interpolate import CubicSpline self._burn_area_spline CubicSpline(webs, areas, bc_typenatural) # 缓存最近的计算结果 self._burn_area_cache {} self._cache_max_size 1000 def burning_area(self, web: float) - float: 使用插值和缓存优化燃面计算 # 检查缓存 if web in self._burn_area_cache: return self._burn_area_cache[web] # 使用样条插值 area float(self._burn_area_spline(web)) # 更新缓存 if len(self._burn_area_cache) self._cache_max_size: self._burn_area_cache[web] area return area7.3 添加高级物理模型class AdvancedMotorModel(Motor): 包含高级物理效应的发动机模型 def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) # 添加侵蚀燃烧效应 self.erosion_enabled kwargs.get(erosion_enabled, False) self.erosion_coefficient kwargs.get(erosion_coeff, 0.0) # 添加两相流损失 self.two_phase_enabled kwargs.get(two_phase_enabled, False) self.particle_size kwargs.get(particle_size, 1e-6) # 颗粒尺寸 # 添加热损失 self.heat_loss_enabled kwargs.get(heat_loss_enabled, False) self.wall_temp kwargs.get(wall_temp, 300.0) # 壁面温度 def burn_rate(self, P, port_velocity0.0): 考虑侵蚀燃烧的燃速模型 r0 super().burn_rate(P) # 基础燃速 if self.erosion_enabled and port_velocity 0: # 侵蚀燃烧模型: r r0 * (1 k * v) r_erosion r0 * (1 self.erosion_coefficient * port_velocity) return r_erosion return r0 def thrust_coefficient(self, Pc, Pe, Pa, gamma): 考虑两相流损失的推力系数 Cf_ideal super().thrust_coefficient(Pc, Pe, Pa, gamma) if self.two_phase_enabled: # 两相流损失估算 loss_factor 0.95 # 假设5%损失 Cf_real Cf_ideal * loss_factor return Cf_real return Cf_ideal7.4 完整的单元测试import unittest import numpy as np class TestSolidMotor(unittest.TestCase): 固体发动机模型测试 def setUp(self): 测试准备 self.propellant Propellant( nameAPCP, density1800.0, burn_rate_params{a: 0.003, n: 0.3}, thermo_params{gamma: 1.2, T0: 2500, MW: 0.024} ) self.grain CylindricalGrain( outer_diameter0.1, length0.3, web_thickness0.02 ) self.motor SolidMotor( propellantself.propellant, grainself.grain, chamber_pressure3e6 ) def test_propellant_burn_rate(self): 测试推进剂燃速计算 P 3e6 # 3 MPa expected_r 0.003 * (3e6) ** 0.3 actual_r self.propellant.burn_rate(P) self.assertAlmostEqual(actual_r, expected_r, delta1e-6) def test_grain_geometry(self): 测试药柱几何计算 # 初始燃面 initial_area self.grain.burning_area(0) self.assertGreater(initial_area, 0) # 最终燃面 final_area self.grain.burning_area(self.grain.web) self.assertGreater(final_area, 0) # 面积应随肉厚增加而减少 area_mid self.grain.burning_area(self.grain.web / 2) self.assertLess(area_mid, initial_area) self.assertGreater(area_mid, final_area) def test_motor_simulation(self): 测试发动机仿真 solver BallisticSolver() results solver.solve(self.motor, max_time5.0) # 检查结果完整性 self.assertIsNotNone(results) self.assertGreater(len(results.time), 10) self.assertGreater(len(results.pressure), 10) self.assertGreater(len(results.thrust), 10) # 检查物理合理性 self.assertTrue(all(p 0 for p in results.pressure)) self.assertTrue(all(t 0 for t in results.thrust)) # 检查总冲为正 total_impulse np.trapz(results.thrust, results.time) self.assertGreater(total_impulse, 0) def test_conservation_approximate(self): 近似检查质量守恒 # 计算生成的总质量 m_propellant self.propellant.density * self.grain.initial_volume # 仿真 solver BallisticSolver() results solver.solve(self.motor, max_time10.0) # 计算通过喷管排出的质量 m_dot_nozzle [P * self.motor.nozzle.At / self.propellant.c_star for P in results.pressure] m_exhausted np.trapz(m_dot_nozzle, results.time) # 排出质量应接近推进剂质量考虑误差 self.assertAlmostEqual(m_exhausted, m_propellant, deltam_propellant*0.1) if __name__ __main__: unittest.main()八、总结该火箭发动机仿真系统实现展示了教育优先、功能完整的设计理念。其代码结构清晰模块化良好适合学习和教学使用。核心算法虽然采用简化模型但准确实现了固体火箭发动机设计的基本原理。关键亮点完整的 MVC 架构分离了计算逻辑和用户界面可扩展的类设计允许用户自定义推进剂、药柱形状实时的可视化反馈增强了用户体验合理的默认参数降低了使用门槛主要局限数值方法简单缺乏工业级求解器的稳健性物理模型过于简化忽略了许多实际效应性能优化不足计算复杂药柱时可能较慢错误处理和边界条件检查不完善最佳改进策略是渐进式重构首先添加全面的单元测试确保现有功能正确性然后逐步替换数值求解核心使用更稳健的算法接着优化性能关键路径特别是几何计算最后添加可选的高级物理模型保持向后兼容性这样的改进路径可以在不破坏现有用户工作流程的前提下显著提升软件的工程实用价值使其从教学工具升级为实用的初步设计工具。