绘制完成

更改后的代码

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
d = 88e-3  # 装药内径 (m)
D = 106e-3  # 装药外径 (m)
L = 200e-3  # 装药长度 (m)
rho_p = 1775  # 装药密度 (kg/m^3)
V_c = 0.0079  # 燃烧室初始自由容积 (m^3)
r_gas = 1.2  # 燃气比热比
c_star = 1500  # 燃气特征速度 (m/s)
d_t = 17e-3  # 喷喉直径 (m)
p_c_t = 1.5e6  # 点火压强 (Pa)
p_a = 0.098e6  # 环境压强 (Pa)
T_c = 3260  # 燃烧温度 (K)
burn_rate_coeff = 8.3e-5  # 燃速系数
burn_rate_exp = 0.3  # 燃速指数

# 计算装药燃烧厚度
e = (D - d) / 2

# 燃速公式
def burn_rate(p_c):
    return burn_rate_coeff * p_c ** burn_rate_exp

# 燃烧时间计算
def burn_time(p_c):
    return e / burn_rate(p_c)

# 零维内弹道方程
def dpcdt(p_c, t):
    if t <= burn_time(p_c):  # 燃烧阶段
        m_dot_p = burn_rate(p_c) * rho_p * np.pi * (D**2 - d**2) / 4 * L
    else:  # 燃烧结束，停止产生燃气
        m_dot_p = 0
    A_t = np.pi * (d_t / 2)**2  # 喷喉面积
    m_dot_g = p_c * A_t / (c_star * np.sqrt(T_c))  # 喷喉质量流量
    return (r_gas - 1) / V_c * (m_dot_p - m_dot_g)

# 四阶龙格-库塔法
def runge_kutta_4th_order(p_c0, t, h):
    p_c = np.zeros(len(t))
    p_c[0] = p_c0
    for i in range(1, len(t)):
        k1 = h * dpcdt(p_c[i-1], t[i-1])
        k2 = h * dpcdt(p_c[i-1] + 0.5 * k1, t[i-1] + 0.5 * h)
        k3 = h * dpcdt(p_c[i-1] + 0.5 * k2, t[i-1] + 0.5 * h)
        k4 = h * dpcdt(p_c[i-1] + k3, t[i-1] + h)
        p_c[i] = p_c[i-1] + (k1 + 2*k2 + 2*k3 + k4) / 6
    return p_c

# 时间数组
t_end = 20  # 总时间 (s)，确保覆盖燃烧和排气阶段
h = 0.01  # 时间步长 (s)
t = np.arange(0, t_end, h)

# 计算燃烧室压强随时间变化
p_c = runge_kutta_4th_order(p_c_t, t, h)

# 绘制燃烧室压强随时间变化曲线
plt.plot(t, p_c / 1e6, label="Chamber Pressure (MPa)")
plt.xlabel("Time (s)")
plt.ylabel("Pressure (MPa)")
plt.title("Chamber Pressure vs Time")
plt.grid()
plt.show()
    

效果图