零维内弹道方程可填写参数生成曲线python代码（软件自备）

import math
import matplotlib.pyplot as plt

# 用户可修改参数区域 ==============================================
engine_params = {
    'throat_area': 0.0001,       # 喷喉面积At [m²]
    'chamber_volume': 0.005,      # 燃烧室容积Vc [m³]
    'prop_density': 1600.0,       # 推进剂密度ρ_p [kg/m³]
    'burn_rate_coeff': 0.0001,    # 燃速系数a [m/(s·Pa^n)]
    'burn_rate_exponent': 0.5,    # 燃速指数n
    'burn_area': 0.1,             # 燃烧面积Ab [m²]
    'initial_pressure': 1e6,      # 初始压力P0 [Pa]
    'propellant_mass': 10.0,      # 推进剂初始质量 [kg]
    'gamma': 1.2,                 # 比热比γ
    'gas_constant': 320.0,        # 气体常数R [J/(kg·K)]
    'flame_temp': 2500.0,         # 火焰温度T0 [K]
    'time_step': 0.001,           # 时间步长dt [s]
    'max_time': 10.0              # 最大模拟时间 [s]
}

# 模拟执行 =======================================================
def simulate_rocket(params):
    # 解包参数
    At = params['throat_area']
    Vc = params['chamber_volume']
    rho_p = params['prop_density']
    a = params['burn_rate_coeff']
    n = params['burn_rate_exponent']
    Ab = params['burn_area']
    P0 = params['initial_pressure']
    m_prop = params['propellant_mass']
    gamma = params['gamma']
    R = params['gas_constant']
    T0 = params['flame_temp']
    dt = params['time_step']
    t_max = params['max_time']

    # 初始化记录列表
    time = []
    pressure = []
    thrust = []
    mass_flow = []
    prop_remain = []

    # 初始状态
    t = 0.0
    P = P0
    m_prop_remaining = m_prop

    while t <= t_max and m_prop_remaining > 0:
        # 计算燃烧速率
        r = a * (P ** n)
        m_dot_prop = rho_p * Ab * r

        # 计算喷管质量流量
        critical_flow_term = math.sqrt(gamma/R) * (2/(gamma+1))**((gamma+1)/(2*(gamma-1)))
        m_dot = P * At / math.sqrt(T0) * critical_flow_term

        # 压力变化计算
        dP_dt = (R*T0/Vc) * (m_dot_prop - m_dot - (gamma*m_dot*R*T0)/((gamma-1)*P))

        # 更新状态
        P += dP_dt * dt
        m_prop_remaining -= m_dot_prop * dt
        current_thrust = m_dot * math.sqrt(2*gamma*R*T0/(gamma-1))  # 简化推力公式

        # 记录数据
        time.append(t)
        pressure.append(P)
        thrust.append(current_thrust)
        mass_flow.append(m_dot)
        prop_remain.append(max(m_prop_remaining, 0))

        t += dt

    return {
        'time': time,
        'pressure': pressure,
        'thrust': thrust,
        'mass_flow': mass_flow,
        'prop_remain': prop_remain
    }

# 结果可视化 =====================================================
def plot_results(results, params):
    fig, axs = plt.subplots(4, 1, figsize=(10, 12))
    
    # 压力曲线
    axs[0].plot(results['time'], [p/1e6 for p in results['pressure']], 'r-')
    axs[0].set_ylabel('Chamber Pressure (MPa)')
    axs[0].grid(True)

    # 推力曲线
    axs[1].plot(results['time'], [f/1e3 for f in results['thrust']], 'b-')
    axs[1].set_ylabel('Thrust (kN)')
    axs[1].grid(True)

    # 质量流量
    axs[2].plot(results['time'], results['mass_flow'], 'g-')
    axs[2].set_ylabel('Mass Flow (kg/s)')
    axs[2].grid(True)

    # 剩余推进剂
    axs[3].plot(results['time'], results['prop_remain'], 'm-')
    axs[3].set_ylabel('Propellant Mass (kg)')
    axs[3].set_xlabel('Time (s)')
    axs[3].grid(True)

    plt.tight_layout()
    plt.suptitle(f"Rocket Engine Simulation (Propellant Mass: {params['propellant_mass']}kg)", y=1.02)
    plt.show()

# 执行模拟并绘图
results = simulate_rocket(engine_params)
plot_results(results, engine_params)

注：执行后会显示4个子图，分别展示：

   - 燃烧室压力变化（MPa）

   - 发动机推力（kN）

   - 喷管质量流量（kg/s）

   - 剩余推进剂质量（kg）