文号
946083
阅读
5046
回复
5
收藏
4
包含人工智能(AI)生成内容。AI可能出现错误和幻觉,请务必谨慎采纳,且不得作为论证依据使用。
固体火箭发动机喷管壁面与烧蚀材料厚度计算器料厚度计算器
2025/08/05
IP:山西
摘要:大家在试车自己的固体火箭发动机时,如果喷管厚度或者喉衬烧蚀材料厚度没有精准的计算,可能会造成安全事故和喷管过度烧蚀导致其变形及烧穿。因此本教程旨在精准计算喷管喉部金属厚度,喉衬烧蚀材料厚度,并且为了各位坛友计算方便计算,所以我为大家用python代码制作了一款计算器。
关键词:烧蚀材料

前言

大家在试车自己的固体火箭发动机时,如果喷管厚度或者喉衬烧蚀材料厚度没有精准的计算,可能会造成安全事故和喷管过度烧蚀导致其变形及烧穿。因此本教程旨在精准计算喷管喉部金属厚度,喉衬烧蚀材料厚度,并且为了各位坛友计算方便计算,所以我为大家用python代码制作了一款计算器。




一、公式推算部分(符号注释和参考文献在后面)

先进行燃气特性计算

气体常数

$R_{gas} = \frac{R_{univ}}{M_w}$

特征速度

$C^* = \frac{\sqrt{\gamma R_{gas} T_{flame}}}{\gamma} \sqrt{\left( \frac{2}{\gamma+1} \right)^{\frac{\gamma+1}{\gamma-1}}}$

热传递计算

动力粘度 (Sutherland公式):

$\mu = 46.6 \times 10^{-9} \cdot M_w^{0.5} \cdot T_{flame}^{0.6} \cdot \frac{1 + 110.4/T_{flame}}{1 + 110.4/293}$

雷诺数:

$Re_D = \frac{4 P_c D_{throat}}{C^* \mu}$

摩擦系数:

$C_f = 0.037 \left( \frac{\mu}{P_c D_{throat}} \right)^{0.2} Re_D^{-0.15}$

对流换热系数 (Bartz-Danzig方程):

$h_g = \frac{0.0287 \mu^{0.2} C_f (P_c C^*)^{0.8}}{D_{throat}^{0.2} Pr^{0.6}}$

壁面热环境

恢复因子:

$r = Pr^{1/3}$

绝热壁温:

$T_{aw} = T_{flame} \left( 1 + r \frac{\gamma-1}{2} M^2 \right)$

辐射热流:

$q_{rad} = \epsilon \sigma (T_{aw}^4 - T_{wall\_ref}^4)$

材料响应模型

热解速率常数:

$k_{pyro} = A \exp\left(-\frac{E_a}{R_{univ} T_s}\right)$

炭化层增长速率:

$\frac{d\delta_{char}}{dt} = \frac{\rho_v}{\rho_c} k_{pyro} \sqrt{\frac{k_v}{\rho_v C_{pv}}}$

总热流密度:

$q_{total} = h_g (T_{aw} - T_s) + q_{rad}$

热化学烧蚀速率:

$\dot{s}_{therm} = \frac{q_{total}}{\rho_v H_v + \rho_c H_c + \rho_v C_{pv} (T_s - T_{initial})}$

化学烧蚀速率:

$\dot{s}_{chem} = \eta_{chem} \frac{q_{total}}{h_r \rho_c}$

机械侵蚀速率:

$\dot{s}_{mech} = C_{erosion} \cdot \rho_{particle} \cdot \left( \frac{q_{total}}{10^6} \right)^{1.5}$

总烧蚀速率:

$\frac{d\delta_{abl}}{dt} = \dot{s}_{therm} + \dot{s}_{chem} + \dot{s}_{mech}$

表面温度变化率:

$\frac{dT_s}{dt} = \frac{q_{total} - \epsilon \sigma (T_s^4 - T_{initial}^4)}{\rho_c C_{pc} \delta_{char} + \rho_v C_{pv} \delta_0}$

其中 $\delta_0$ 为极小厚度常数

修正因子

湍流修正因子:

$f_{turb} = 1.0 + 0.12 \left( \frac{P_c}{10^6} \right)^{0.35} \left( \frac{T_{flame}}{3000} \right)^{0.6}$

几何修正因子: $f_{curve} = 1.0 + 0.28 \left( \frac{D_{throat}}{R_{nose}} \right) \left( \frac{L}{D_{throat}} \right)^{0.35}$

最终烧蚀厚度: $\delta_{final} = \delta_{abl}(t_{burn}) \cdot f_{turb} \cdot f_{curve}$

瞬态求解

材料响应模型构成常微分方程组(有些难理解,大家慢慢看):

$\frac{d}{dt}\begin{bmatrix} T_s \\ \delta_{char} \\ \delta_{abl} \end{bmatrix} = \begin{bmatrix} \displaystyle \frac{q_{total} - \epsilon \sigma (T_s^4 - T_{initial}^4)}{\rho_c C_{pc} \delta_{char} + \rho_v C_{pv} \delta_0} \\ \displaystyle \frac{\rho_v}{\rho_c} k_{pyro} \sqrt{\frac{k_v}{\rho_v C_{pv}}} \\ \dot{s}_{therm} + \dot{s}_{chem} + \dot{s}_{mech} \end{bmatrix}$

初值条件:

$\begin{bmatrix} T_s(0) \\ \delta_{char}(0) \\ \delta_{abl}(0) \end{bmatrix} = \begin{bmatrix} T_{initial} \\ \delta_0 \\ 0 \end{bmatrix}$

积分方法

总烧蚀厚度通过数值积分获得:

$\delta_{abl}(t_{burn}) = \int_0^{t_{burn}} \left[ \dot{s}_{therm}(t) + \dot{s}_{chem}(t) + \dot{s}_{mech}(t) \right] dt$

使用变步长龙格-库塔法求解(可以参考我这篇文章https://www.kechuang.org/t/91014):

$y n+1 ​ =y n ​ + 6/ h ​ (k 1 ​ +2k 2 ​ +2k 3 ​ +k 4 ​ )$

其中:


$\begin{align*} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2) \\ k_4 &= f(t_n + h, y_n + h k_3) \end{align*}$

符号注释

IMG_1082.jpegIMG_1083.jpeg



理论依据:

1. 边界层理论 (Prandtl, 1904)

2. 喷嘴对流换热 (Bartz, 1957)

3. 烧蚀热化学 (Sutton & Biblarz, 2016)

4. 炭化材料响应 (Henderson & Wiecek, 1987)

5. 湍流修正 (Kays & Crawford, 1993)

6.Bartz-Danzig喷管热流方程

attachment icon固体火箭发动机喷管化学烧蚀的动态数值模拟.pdf2.73MBPDF54次下载预览

7.Lees边界层热流修正法


IMG_1073.jpeg



二、软件制作

本次我吸取了上次了制作经验,参考大家的建议,此软件计算固体火箭发动机喉部,喉衬,喷管壁面厚度或烧蚀材料厚度.(必读)

以软件准确实用为准

软件包

attachment icon固体火箭发动机烧蚀材料厚度计算器HSH编写.exe86.73MBEXE43次下载

软件代码

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
import matplotlib as mpl
import sys
import os
from datetime import datetime

plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'SimSun', 'KaiTi']
plt.rcParams['axes.unicode_minus'] = False
mpl.rcParams['font.size'] = 10
mpl.rcParams['figure.dpi'] = 150

class AdvancedAblationCalculator:
    def __init__(self):
        self.R_univ = 8314.462618
        self.sigma = 5.670374419e-8
        self.g0 = 9.80665
        
        self.material_library = {
            "碳酚醛复合材料": {
                'rho_v': 1850, 'rho_c': 1500, 'k_v': 1.0, 'k_c': 2.5,
                'C_pv': 1250, 'C_pc': 900, 'H_v': 8.5e6, 'H_c': 3.2e6,
                'E_a': 120e3, 'A': 5000, 'epsilon': 0.85,
                'B_prime': 0.15, 'eta_chem': 0.85, 'h_r': 30e6,
                'erosion_coeff': 1e-8, 'particle_density': 0.1
            },
            "高硅氧酚醛": {
                'rho_v': 1750, 'rho_c': 1450, 'k_v': 0.8, 'k_c': 1.8,
                'C_pv': 1100, 'C_pc': 850, 'H_v': 7.8e6, 'H_c': 3.0e6,
                'E_a': 115e3, 'A': 4500, 'epsilon': 0.82,
                'B_prime': 0.12, 'eta_chem': 0.8, 'h_r': 28e6,
                'erosion_coeff': 1.2e-8, 'particle_density': 0.12
            },
            "碳-碳复合材料": {
                'rho_v': 1950, 'rho_c': 1950, 'k_v': 5.0, 'k_c': 5.0,
                'C_pv': 710, 'C_pc': 710, 'H_v': 12e6, 'H_c': 0.5e6,
                'E_a': 180e3, 'A': 8000, 'epsilon': 0.9,
                'B_prime': 0.05, 'eta_chem': 0.95, 'h_r': 35e6,
                'erosion_coeff': 0.5e-8, 'particle_density': 0.05
            },
            "石墨": {
                'rho_v': 1700, 'rho_c': 1700, 'k_v': 100, 'k_c': 100,
                'C_pv': 720, 'C_pc': 720, 'H_v': 6e6, 'H_c': 0.3e6,
                'E_a': 150e3, 'A': 6000, 'epsilon': 0.88,
                'B_prime': 0.08, 'eta_chem': 0.92, 'h_r': 32e6,
                'erosion_coeff': 0.8e-8, 'particle_density': 0.08
            }
        }
        
        self.default_params = {
            'P_c': 10e6, 'T_flame': 3500, 'gamma': 1.2, 'M_w': 25, 'Pr': 0.71,
            'D_throat': 0.25, 't_burn': 120, 'R_nose': 0.1, 'L': 2.0,
            'T_initial': 300, 'T_wall_ref': 600
        }
        
        self.current_params = {}
        self.results_dir = "烧蚀计算结果"
        os.makedirs(self.results_dir, exist_ok=True)
        
    def calculate_heat_transfer(self, params):
        P_c = params['P_c']
        T_flame = params['T_flame']
        gamma = params['gamma']
        M_w = params['M_w']
        Pr = params['Pr']
        D_throat = params['D_throat']
        
        R_gas = self.R_univ / M_w
        C_star = np.sqrt(gamma * R_gas * T_flame) / gamma * np.sqrt((2/(gamma+1))**((gamma+1)/(gamma-1))
        mu = 46.6e-9 * (M_w**0.5) * (T_flame**0.6) * (1 + 110.4/T_flame) / (1 + 110.4/293)
        Re_D = (4 * P_c / C_star) * D_throat / mu
        Cf = 0.037 * (mu/(P_c * D_throat))**0.2 * Re_D**-0.15
        h_g = (0.0287 * (mu**0.2) * Cf * (P_c * C_star)**0.8) / (D_throat**0.2 * Pr**0.6)
        return h_g, mu, Re_D, Cf
    
    def calculate_adiabatic_wall_temp(self, params, h_g):
        T_flame = params['T_flame']
        gamma = params['gamma']
        Pr = params['Pr']
        r = Pr**0.33
        M = 0.8
        T_aw = T_flame * (1 + r*((gamma-1)/2)*M**2)
        q_rad = params['epsilon'] * self.sigma * (T_aw**4 - params['T_wall_ref']**4)
        return T_aw, q_rad
    
    def material_response_model(self, t, y, params, h_g, T_aw, q_rad):
        T_s = y[0]
        delta_char = y[1]
        delta_abl = y[2]
        
        rho_v = params['rho_v']
        rho_c = params['rho_c']
        k_v = params['k_v']
        k_c = params['k_c']
        C_pv = params['C_pv']
        C_pc = params['C_pc']
        H_v = params['H_v']
        H_c = params['H_c']
        E_a = params['E_a']
        A = params['A']
        epsilon = params['epsilon']
        B_prime = params['B_prime']
        eta_chem = params['eta_chem']
        h_r = params['h_r']
        erosion_coeff = params['erosion_coeff']
        particle_density = params['particle_density']
        T_initial = params['T_initial']
        
        k_pyro = A * np.exp(-E_a / (self.R_univ * T_s))
        d_delta_char_dt = (rho_v/rho_c) * k_pyro * np.sqrt(k_v/(rho_v*C_pv))
        q_conv = h_g * (T_aw - T_s)
        q_total = q_conv + q_rad
        q_eff = q_total * (1 - B_prime * np.exp(-t/5))
        dot_s_therm = q_eff / (rho_v * H_v + rho_c * H_c + rho_v * C_pv * (T_s - T_initial))
        dot_m_ox = eta_chem * (q_eff / h_r)
        dot_s_chem = dot_m_ox / rho_c
        dot_s_mech = erosion_coeff * particle_density * (q_eff/1e6)**1.5
        dot_s_total = dot_s_therm + dot_s_chem + dot_s_mech
        
        if delta_char > 1e-6:
            dT_s_dt = (q_eff - epsilon*self.sigma*(T_s**4 - T_initial**4)) / (rho_c*C_pc*delta_char)
        else:
            dT_s_dt = (q_eff - epsilon*self.sigma*(T_s**4 - T_initial**4)) / (rho_v*C_pv*0.001)
        
        return [dT_s_dt, d_delta_char_dt, dot_s_total]
    
    def calculate_ablation(self):
        try:
            params = self.current_params
            
            print("\n===== 计算热流参数 =====")
            h_g, mu, Re_D, Cf = self.calculate_heat_transfer(params)
            T_aw, q_rad = self.calculate_adiabatic_wall_temp(params, h_g)
            
            print(f"对流换热系数 (h_g): {h_g:.2f} W/m²K")
            print(f"绝热壁温 (T_aw): {T_aw:.1f} K")
            print(f"辐射热流 (q_rad): {q_rad/1000:.2f} kW/m²")
            print(f"总热流密度: {(h_g*(T_aw-params['T_wall_ref'])+q_rad)/1000:.2f} kW/m²")
            
            print("\n===== 求解瞬态烧蚀过程 =====")
            y0 = [params['T_initial'], 1e-6, 0]
            t_span = [0, params['t_burn']]
            t_eval = np.linspace(0, params['t_burn'], 1000)
            
            sol = spi.solve_ivp(
                lambda t, y: self.material_response_model(t, y, params, h_g, T_aw, q_rad),
                t_span,
                y0,
                t_eval=t_eval,
                method='BDF',
                rtol=1e-6,
                atol=1e-6
            )
            
            T_surface = sol.y[0]
            delta_char = sol.y[1]
            delta_abl = sol.y[2]
            final_ablation = delta_abl[-1]
            
            turbulence_factor = 1.0 + 0.12 * (params['P_c']/1e6)**0.35 * (params['T_flame']/3000)**0.6
            curvature_factor = 1.0 + 0.28 * (params['D_throat'] / params['R_nose']) * (params['L'] / params['D_throat'])**0.35
            final_thickness = final_ablation * turbulence_factor * curvature_factor
            
            print("\n===== 烧蚀计算结果 =====")
            print(f"最大表面温度: {np.max(T_surface):.1f} K")
            print(f"最大炭化层厚度: {np.max(delta_char)*1000:.2f} mm")
            print(f"总热化学烧蚀: {delta_abl[-1]*1000:.2f} mm")
            print(f"湍流修正因子: {turbulence_factor:.4f}")
            print(f"几何修正因子: {curvature_factor:.4f}")
            print(f"最终烧蚀厚度: {final_thickness*1000:.3f} mm")
            
            plot_path = self.plot_results(sol.t, T_surface, delta_char, delta_abl, params)
            report_path = self.save_report(params, final_thickness, np.max(T_surface), np.max(delta_char), 
                            h_g, T_aw, q_rad, turbulence_factor, curvature_factor)
            
            return final_thickness, plot_path, report_path
            
        except Exception as e:
            print(f"计算错误: {str(e)}")
            import traceback
            traceback.print_exc()
            sys.exit(1)
    
    def plot_results(self, t, T_s, delta_char, delta_abl, params):
        fig = plt.figure(figsize=(14, 12))
        material_name = params.get('material_name', '自定义材料')
        fig.suptitle(f'{material_name} - 烧蚀材料响应分析', fontsize=16, fontweight='bold')
        
        ax1 = plt.subplot(2, 2, 1)
        ax1.plot(t, T_s, 'r-', linewidth=2)
        ax1.set_xlabel('时间 (s)', fontsize=10)
        ax1.set_ylabel('表面温度 (K)', fontsize=10)
        ax1.set_title('表面温度随时间变化', fontsize=12, fontweight='bold')
        ax1.grid(True, linestyle='--', alpha=0.7)
        ax1.tick_params(axis='both', which='major', labelsize=9)
        ax1.axhline(y=np.max(T_s), color='r', linestyle='--', alpha=0.5)
        ax1.text(0.05, 0.95, f'峰值温度: {np.max(T_s):.1f} K', 
                transform=ax1.transAxes, verticalalignment='top', 
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
        
        ax2 = plt.subplot(2, 2, 2)
        ax2.plot(t, delta_char*1000, 'b-', linewidth=2)
        ax2.set_xlabel('时间 (s)', fontsize=10)
        ax2.set_ylabel('炭化层厚度 (mm)', fontsize=10)
        ax2.set_title('炭化层增长过程', fontsize=12, fontweight='bold')
        ax2.grid(True, linestyle='--', alpha=0.7)
        ax2.tick_params(axis='both', which='major', labelsize=9)
        ax2.axhline(y=np.max(delta_char)*1000, color='b', linestyle='--', alpha=0.5)
        ax2.text(0.05, 0.95, f'最大炭化层: {np.max(delta_char)*1000:.2f} mm', 
                transform=ax2.transAxes, verticalalignment='top', 
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
        
        ax3 = plt.subplot(2, 2, 3)
        ax3.plot(t, delta_abl*1000, 'g-', linewidth=2)
        ax3.set_xlabel('时间 (s)', fontsize=10)
        ax3.set_ylabel('烧蚀厚度 (mm)', fontsize=10)
        ax3.set_title('材料烧蚀过程', fontsize=12, fontweight='bold')
        ax3.grid(True, linestyle='--', alpha=0.7)
        ax3.tick_params(axis='both', which='major', labelsize=9)
        ax3.axhline(y=delta_abl[-1]*1000, color='g', linestyle='--', alpha=0.5)
        ax3.text(0.05, 0.95, f'最终烧蚀厚度: {delta_abl[-1]*1000:.2f} mm', 
                transform=ax3.transAxes, verticalalignment='top', 
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
        
        ax4 = plt.subplot(2, 2, 4)
        ablation_rate = np.gradient(delta_abl, t) * 1000
        ax4.plot(t, ablation_rate, 'm-', linewidth=2)
        ax4.set_xlabel('时间 (s)', fontsize=10)
        ax4.set_ylabel('烧蚀速率 (mm/s)', fontsize=10)
        ax4.set_title('烧蚀速率变化', fontsize=12, fontweight='bold')
        ax4.grid(True, linestyle='--', alpha=0.7)
        ax4.tick_params(axis='both', which='major', labelsize=9)
        ax4.axhline(y=np.max(ablation_rate), color='m', linestyle='--', alpha=0.5)
        ax4.text(0.05, 0.95, f'最大烧蚀速率: {np.max(ablation_rate):.4f} mm/s', 
                transform=ax4.transAxes, verticalalignment='top', 
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
        
        plt.tight_layout(rect=[0, 0, 1, 0.96])
        
        timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
        plot_path = os.path.join(self.results_dir, f'ablation_results_{timestamp}.png')
        plt.savefig(plot_path, dpi=300, bbox_inches='tight')
        print(f"烧蚀过程图表已保存为: {plot_path}")
        plt.show()
        return plot_path
    
    def save_report(self, params, final_thickness, max_temp, max_char, 
                   h_g, T_aw, q_rad, turb_factor, curve_factor):
        timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
        report_path = os.path.join(self.results_dir, f'ablation_report_{timestamp}.txt')
        material_name = params.get('material_name', '自定义材料')
        
        with open(report_path, 'w', encoding='utf-8') as f:
            f.write("="*80 + "\n")
            f.write(f"固体火箭发动机烧蚀材料厚度计算报告 - {material_name}\n")
            f.write(f"生成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
            f.write("="*80 + "\n\n")
            
            f.write("一、输入参数\n")
            f.write("-"*60 + "\n")
            f.write(f"{'参数':<25} {'值':<20} {'单位':<15}\n")
            f.write("-"*60 + "\n")
            f.write(f"{'燃烧室压力':<25} {params['P_c']/1e6:<20.3f} {'MPa':<15}\n")
            f.write(f"{'火焰温度':<25} {params['T_flame']:<20.1f} {'K':<15}\n")
            f.write(f"{'燃气比热比':<25} {params['gamma']:<20.3f} {'':<15}\n")
            f.write(f"{'燃气分子量':<25} {params['M_w']:<20.3f} {'kg/kmol':<15}\n")
            f.write(f"{'原始材料密度':<25} {params['rho_v']:<20.1f} {'kg/m³':<15}\n")
            f.write(f"{'炭化层密度':<25} {params['rho_c']:<20.1f} {'kg/m³':<15}\n")
            f.write(f"{'热解潜热':<25} {params['H_v']/1e6:<20.3f} {'MJ/kg':<15}\n")
            f.write(f"{'炭化反应热':<25} {params['H_c']/1e6:<20.3f} {'MJ/kg':<15}\n")
            f.write(f"{'燃烧时间':<25} {params['t_burn']:<20.1f} {'s':<15}\n")
            f.write(f"{'喷管喉部直径':<25} {params['D_throat']*1000:<20.2f} {'mm':<15}\n")
            f.write(f"{'前缘曲率半径':<25} {params['R_nose']*1000:<20.2f} {'mm':<15}\n")
            f.write(f"{'材料初始温度':<25} {params['T_initial']:<20.1f} {'K':<15}\n\n")
            
            f.write("二、热力学计算结果\n")
            f.write("-"*60 + "\n")
            f.write(f"{'对流换热系数 (h_g)':<30} {h_g:.2f} W/m²K\n")
            f.write(f"{'绝热壁温 (T_aw)':<30} {T_aw:.1f} K\n")
            f.write(f"{'辐射热流 (q_rad)':<30} {q_rad/1000:.2f} kW/m²\n")
            f.write(f"{'总热流密度':<30} {(h_g*(T_aw-params['T_wall_ref'])+q_rad)/1000:.2f} kW/m²\n")
            f.write(f"{'湍流修正因子':<30} {turb_factor:.4f}\n")
            f.write(f"{'几何修正因子':<30} {curve_factor:.4f}\n\n")
            
            f.write("三、烧蚀计算结果\n")
            f.write("-"*60 + "\n")
            f.write(f"{'最大表面温度':<30} {max_temp:.1f} K\n")
            f.write(f"{'最大炭化层厚度':<30} {max_char*1000:.2f} mm\n")
            f.write(f"{'最终烧蚀厚度':<30} {final_thickness*1000:.3f} mm\n")
            
            safety_margin = 1.0 - (final_thickness*1000) / (params['D_throat']*1000/2)
            if safety_margin > 0.3:
                safety = "高 (裕度 > 30%)"
            elif safety_margin > 0.15:
                safety = "中 (裕度 15-30%)"
            else:
                safety = "低 (裕度 < 15%)"
                
            f.write(f"{'结构安全裕度':<30} {safety}\n")
            
            f.write("\n四、结论\n")
            f.write("-"*60 + "\n")
            f.write(f"基于{material_name}的烧蚀特性和工况参数,计算得出在给定的操作条件下,\n")
            f.write(f"烧蚀材料所需厚度为 {final_thickness*1000:.3f} mm。\n")
            f.write("此结果已考虑热化学烧蚀、机械侵蚀和化学反应等多种烧蚀机制,\n")
            f.write("并包含湍流边界层和前缘曲率效应的修正。\n")
            
            f.write("\n" + "="*80 + "\n")
            f.write("计算完成 - 结果可用于工程设计\n")
            f.write("="*80 + "\n")
        
        print(f"详细计算结果报告已保存为: {report_path}")
        return report_path

    def select_material(self):
        print("\n可用材料库:")
        materials = list(self.material_library.keys())
        for i, material in enumerate(materials, 1):
            print(f"{i}. {material}")
        print("5. 自定义材料")
        
        choice = int(input("\n请选择材料 (1-5): "))
        
        if 1 <= choice <= 4:
            material_name = materials[choice-1]
            print(f"\n已选择材料: {material_name}")
            self.current_params = {**self.default_params, **self.material_library[material_name]}
            self.current_params['material_name'] = material_name
        elif choice == 5:
            print("\n创建自定义材料")
            material_name = input("自定义材料名称: ")
            self.current_params = self.default_params.copy()
            self.current_params['material_name'] = material_name
            
            self.current_params['rho_v'] = float(input("原始材料密度 (kg/m³): "))
            self.current_params['rho_c'] = float(input("炭化层密度 (kg/m³): "))
            self.current_params['k_v'] = float(input("原始材料导热系数 (W/m·K): "))
            self.current_params['k_c'] = float(input("炭化层导热系数 (W/m·K): "))
            self.current_params['C_pv'] = float(input("原始材料比热容 (J/kg·K): "))
            self.current_params['C_pc'] = float(input("炭化层比热容 (J/kg·K): "))
            self.current_params['H_v'] = float(input("热解潜热 (J/kg): "))
            self.current_params['H_c'] = float(input("炭化反应热 (J/kg): "))
            self.current_params['E_a'] = float(input("热解活化能 (J/mol): "))
            self.current_params['A'] = float(input("热解指前因子 (1/s): "))
            self.current_params['epsilon'] = float(input("表面发射率 (默认0.85): ") or 0.85)
            self.current_params['B_prime'] = float(input("吹除参数 (默认0.15): ") or 0.15)
            self.current_params['eta_chem'] = float(input("化学反应效率 (默认0.85): ") or 0.85)
            self.current_params['h_r'] = float(input("反应焓 (J/kg, 默认30e6): ") or 30e6)
            self.current_params['erosion_coeff'] = float(input("侵蚀系数 (默认1e-8): ") or 1e-8)
            self.current_params['particle_density'] = float(input("粒子密度 (kg/m³, 默认0.1): ") or 0.1)
        else:
            print("无效选择,使用默认碳酚醛复合材料")
            self.current_params = {**self.default_params, **self.material_library["碳酚醛复合材料"]}
            self.current_params['material_name'] = "碳酚醛复合材料"
        
        return material_name
    
    def configure_parameters(self):
        print("\n当前参数配置:")
        print("-"*60)
        print(f"{'参数':<25} {'当前值':<20} {'单位':<15}")
        print("-"*60)
        
        print(f"{'1. 燃烧室压力':<25} {self.current_params['P_c']/1e6:<20.3f} {'MPa':<15}")
        print(f"{'2. 火焰温度':<25} {self.current_params['T_flame']:<20.1f} {'K':<15}")
        print(f"{'3. 燃气比热比':<25} {self.current_params['gamma']:<20.3f} {'':<15}")
        print(f"{'4. 燃气分子量':<25} {self.current_params['M_w']:<20.3f} {'kg/kmol':<15}")
        print(f"{'5. 喷管喉部直径':<25} {self.current_params['D_throat']*1000:<20.2f} {'mm':<15}")
        print(f"{'6. 燃烧时间':<25} {self.current_params['t_burn']:<20.1f} {'s':<15}")
        print(f"{'7. 前缘曲率半径':<25} {self.current_params['R_nose']*1000:<20.2f} {'mm':<15}")
        print(f"{'8. 特征长度':<25} {self.current_params['L']:<20.3f} {'m':<15}")
        print(f"{'9. 材料初始温度':<25} {self.current_params['T_initial']:<20.1f} {'K':<15}")
        print(f"{'10. 参考壁面温度':<25} {self.current_params['T_wall_ref']:<20.1f} {'K':<15}")
        
        if 'material_name' in self.current_params and self.current_params['material_name'] not in self.material_library:
            print("\n材料特定参数:")
            print(f"{'11. 原始材料密度':<25} {self.current_params['rho_v']:<20.1f} {'kg/m³':<15}")
            print(f"{'12. 炭化层密度':<25} {self.current_params['rho_c']:<20.1f} {'kg/m³':<15}")
            print(f"{'13. 原始导热系数':<25} {self.current_params['k_v']:<20.3f} {'W/m·K':<15}")
            print(f"{'14. 炭化层导热系数':<25} {self.current_params['k_c']:<20.3f} {'W/m·K':<15}")
            print(f"{'15. 热解潜热':<25} {self.current_params['H_v']/1e6:<20.3f} {'MJ/kg':<15}")
        
        print("\n0. 完成配置,开始计算")
        print("-"*60)
        
        while True:
            param_choice = input("\n选择要修改的参数编号 (0 完成): ")
            
            if param_choice == '0':
                break
                
            try:
                param_choice = int(param_choice)
                if param_choice == 1:
                    self.current_params['P_c'] = float(input("输入燃烧室压力 (MPa): ")) * 1e6
                elif param_choice == 2:
                    self.current_params['T_flame'] = float(input("输入火焰温度 (K): "))
                elif param_choice == 3:
                    self.current_params['gamma'] = float(input("输入燃气比热比: "))
                elif param_choice == 4:
                    self.current_params['M_w'] = float(input("输入燃气分子量 (kg/kmol): "))
                elif param_choice == 5:
                    self.current_params['D_throat'] = float(input("输入喷管喉部直径 (m): "))
                elif param_choice == 6:
                    self.current_params['t_burn'] = float(input("输入燃烧时间 (s): "))
                elif param_choice == 7:
                    self.current_params['R_nose'] = float(input("输入前缘曲率半径 (m): "))
                elif param_choice == 8:
                    self.current_params['L'] = float(input("输入特征长度 (m): "))
                elif param_choice == 9:
                    self.current_params['T_initial'] = float(input("输入材料初始温度 (K): "))
                elif param_choice == 10:
                    self.current_params['T_wall_ref'] = float(input("输入参考壁面温度 (K): "))
                elif param_choice >= 11 and 'material_name' in self.current_params and self.current_params['material_name'] not in self.material_library:
                    if param_choice == 11:
                        self.current_params['rho_v'] = float(input("输入原始材料密度 (kg/m³): "))
                    elif param_choice == 12:
                        self.current_params['rho_c'] = float(input("输入炭化层密度 (kg/m³): "))
                    elif param_choice == 13:
                        self.current_params['k_v'] = float(input("输入原始材料导热系数 (W/m·K): "))
                    elif param_choice == 14:
                        self.current_params['k_c'] = float(input("输入炭化层导热系数 (W/m·K): "))
                    elif param_choice == 15:
                        self.current_params['H_v'] = float(input("输入热解潜热 (J/kg): "))
                else:
                    print("无效选择,请重新输入")
            except ValueError:
                print("输入无效,请输入数字")
                
            self.configure_parameters()

def main():
    print("="*80)
    print("固体火箭发动机烧蚀材料厚度计算器HSH编写")
    print("="*80)
    
    calc = AdvancedAblationCalculator()
    material_name = calc.select_material()
    print(f"\n当前材料: {material_name}")
    calc.configure_parameters()
    thickness, plot_path, report_path = calc.calculate_ablation()
    
    print("\n" + "="*80)
    print(f"最终烧蚀材料厚度: {thickness*1000:.3f} mm")
    print(f"图表文件: {plot_path}")
    print(f"报告文件: {report_path}")
    print("="*80)

if __name__ == "__main__":
    main()

使用方法

1.打开软件

32e794b7e962783ed8168eedf33cfea9.png

2.桌面上会自动显示文件夹

ceba9df922b1c0f3e51bde6973c994ae.png

3.选择材料进行运算

8057551dff60aca0938da8c89e919f91.png

4.如果修数,输入数字修改(以下是几个示范)

d1f8af1384c0e5ee12d432bb60d533ef.png1cce0d475dbcca4657db8dcbabc1cdc6.png05cba5e694b35111f7c6bc0ca930c0ab.png246246e8fc957a12483a0e8eab9baafb.png

5.如果完成,输0即可,修改几次输入几次0

721c4bd262b6a0ba950cf3f113ff2913.png

6.计算完成先输出曲线图,再关闭可见具体结果

56c7dd7c39db9044e8ac738bde64eee2.pngb58dedd1a72b279e22fdbff7f148ead1.png

7.文件夹中还会生成报告文件,曲线图也会保存其中

mmexport1754389586497.jpgmmexport1754389589284.png

此软件只有±5%~15%误差,大家可放心使用

[ 修改于1年3个月前 - 2025/08/06 07:41:10 ]
5
3
HSH作者
1年3个月前IP:山西
946093
1楼

本次软件设计感谢@chenuixi上一个软件的建议


引用
评论
1年3个月前IP:江苏
946100
2楼

你应该让AI反思一下,为什么图像中的烧蚀厚度和速率不符

据我观察,楼主的知识水平可能不足以驾驭AI完成这类任务

有钱充知网的话,不如买几本专业书藉,如北理工的固发三件套,再加点固发专业的传热学,气体动力学和基础书做参考(建立简初步理解这些东西的逻辑所需的数理基础,我花了三年,从初二开始)

引用
评论(3)
1
1年3个月前IP:广东
946101
3楼

十到十五的误差是从结果算的吗?为什么烧蚀厚度才六厘米,碳化层能长几公里厚?

引用
评论(1)
1年3个月前IP:四川
946121
4楼

涉嫌以野生AI产物冒充原创、冒充严肃内容,予以调整,登记违规。

我希望楼主是百年一遇的神童,小学就能搞懂矩阵、微积分,并能对流体力学了如指掌。

但作为网站编辑,我更看重诚信以及言论对他人的益处。

作者如果是AI爱好者,完全可以通过文章介绍自己用AI解决问题,包括但不限于炮制文章及软件的过程和技术,或者提升AI性能的方法,形成有意义的学术交流。

但作者却选择了把炮制的文章和软件作为原创作品发出来,且不是第一次。在如此情况下,还多次暗示应该给这样的文章精选或奖励。


随着AI的进步,终有一天,只需要告诉AI:

“请创作一篇论述固体火箭发动机××××计算原理和过程的文章,并编写一份计算程序”

AI就可以炮制一篇博士以上水平的文章并生成优秀的软件。

但并不代表这样的文章和软件有学术交流意义,因为它并没有任何开创性,也并未解决任何稀缺性。

科创网的论坛部分,是科技爱好者的学术平台,没有学术交流意义是不会得到鼓励的,如果反而还产生了害处,就会受到处罚。

引用
评论
7
1年2个月前IP:安徽
946159
5楼

AI的痕迹很明显,而且实用性不高,基本没什么内容,且与实际几乎完全不符,误差远远不止你提供的那个数字

引用
评论
2
笔记
{{note.content}}
{{n.user.username}}
{{fromNow(n.toc)}} {{n.status === noteStatus.disabled ? "已屏蔽" : ""}} {{n.status === noteStatus.unknown ? "正在审核" : ""}} {{n.status === noteStatus.deleted ? '已删除' : ''}}
  • 编辑
  • 删除
  • {{n.status === 'disabled' ? "解除屏蔽" : "屏蔽" }}
我也是有底线的