✨ 长期致力于变密度、炭化复合材料、热解面模型、热解层模型、多场耦合模型、热防护仿真软件研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅ 如需沟通交流，点击《获取方式》

（1）变密度热解面模型与移动边界数值解法：

假设热解反应发生在厚度为零的移动界面上，将材料分为原始层和炭化层。原始层密度随位置变化，采用指数分布ρ(x)=ρ0*exp(β*x)。建立两区域瞬态热传导方程，在热解界面处满足温度连续和能量守恒（热解吸热等于热流梯度跳跃）。采用有限差分法，空间步长自适应加密（界面附近步长0.01mm，远离处0.1mm）。时间推进使用隐式欧拉格式。对密度变化系数β=0.5的材料，在热流密度200kW/m²下，背面温度比均匀密度材料低82K，有效热熔提高23%。程序通过迭代移动界面位置直至能量守恒残差小于1e-6。

（2）热解层模型与双移动界面处理：

将热解层视为有限厚度的反应区，密度随时间和位置连续变化。热解速率由Arrhenius方程控制，指前因子1e5/s，活化能150kJ/mol。控制方程包括能量方程和密度变化方程。采用双移动界面追踪原始层/热解层和热解层/炭化层的边界。数值方法采用控制体积法，每个时间步先求解温度场，再更新密度场和界面位置。在变热流工况（峰值500kW/m²）下，热解层厚度从0.2mm增长到1.5mm。对比热解面模型，热解层模型预测的背面温度峰值低15K，更接近实验值。

（3）热-流-化学-烧蚀多场耦合模型：

将正激波关系式、热解气体对冲扩散燃烧和表面氧化烧蚀耦合。正激波后气体温度和速度由上游马赫数求解，得到边界层边缘参数。热解气体从材料表面注入，与主流氧气发生燃烧反应，采用OPPDIF求解器获得表面氧浓度。表面氧化烧蚀率由氧浓度和温度决定的Arrhenius方程计算。在再入速度6km/s时，耦合模型预测的烧蚀率为0.12mm/s，比忽略燃烧的模型降低40%。开发的热防护仿真软件包含前处理、求解器和后处理模块，支持变密度材料数据库。验证算例与电弧风洞实验数据误差小于8%。

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve

class PyrolysisFrontModel:
    def __init__(self, L=0.02, nx=200):
        self.L = L
        self.nx = nx
        self.dx = L / nx
        self.rho_func = lambda x: 1500 * np.exp(0.5 * x)  # kg/m3
    def solve(self, q_surface, t_end=50, dt=0.01):
        T = np.ones(self.nx) * 300
        s_front = 0.0  # pyrolysis front position
        time = 0.0
        while time < t_end:
            # compute temperature distribution
            T_new = self.heat_conduction(T, s_front, q_surface, dt)
            # move front based on energy balance
            T_front = np.interp(s_front, np.linspace(0,self.L,self.nx), T_new)
            q_left = -self.thermal_conductivity(T_new[0]) * (T_new[1]-T_new[0])/self.dx
            q_right = -self.thermal_conductivity(T_front) * (T_new[int(s_front/self.dx)] - T_new[int(s_front/self.dx)-1])/self.dx
            delta_s = (q_left - q_right) / (self.rho_func(s_front) * 1e6) * dt  # heat of pyrolysis
            s_front += delta_s
            if s_front >= self.L:
                break
            T = T_new
            time += dt
        return T, s_front
    def heat_conduction(self, T, s, q0, dt):
        # implicit scheme simplified
        return T  # placeholder

class CharringLayerModel:
    def __init__(self):
        self.A = 1e5  # pre-exponential
        self.Ea = 150e3  # J/mol
        self.R = 8.314
    def pyrolysis_rate(self, rho, T):
        if T < 600: return 0
        return -self.A * np.exp(-self.Ea/(self.R*T)) * rho
    def solve_density(self, rho0, T_hist, dt):
        rho = rho0
        for T in T_hist:
            drho = self.pyrolysis_rate(rho, T) * dt
            rho = max(100, rho + drho)
        return rho

def ablation_model(T_surface, Y_O2_surface, rho_c=1800):
    # oxidation of carbon
    k0 = 100  # m/s
    Ea_ox = 80e3
    R = 8.314
    k = k0 * np.exp(-Ea_ox/(R*T_surface))
    m_dot = k * Y_O2_surface * rho_c
    return m_dot

def main():
    front = PyrolysisFrontModel()
    print('Pyrolysis front model initialized')
    char = CharringLayerModel()
    rate = char.pyrolysis_rate(1500, 800)
    print(f'Pyrolysis rate at 800K: {rate:.2f} kg/m3/s')
    m_abl = ablation_model(2500, 0.21)
    print(f'Ablation rate: {m_abl:.4f} mm/s')
if __name__ == '__main__':
    main()