固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法

摘要

本发明公开了一种固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法，建立了C/C复材喷管的多尺度烧蚀模型，综合考虑了燃气流动、传热和传质，粒子的机械侵蚀、燃气的化学反应等过程，获得了喷管的温度场、流场、压力分布情况，通过粒子轨迹的结果分析判断喷管各部位的机械侵蚀程度，并且利用水平集的方法模拟了C/C复合材料微观表面的退移和形貌，揭示了纤维和基体的热化学烧蚀规律。本发明从宏观和微观两个角度全面分析了C/C复合材料喷管的烧蚀情况，对喷管的热性能研究和精细化设计有重要价值。

权利要求

1.一种固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法,包括以下步骤：步骤1：建立C/C复合材料的喷管几何模型；

采用通道截面形状先收敛后扩张的拉瓦尔喷管，基于喷管及复合材料结构的对称性，建立二维轴对称几何模型，从下到上依次是：收敛段、喉部、扩张段，外端为C/C复合材料；

步骤2：建立气相控制方程；

设燃料气体中的冷凝相颗粒的质量含量为30wt％,质量流率为0.01kg/s；燃烧室压力和温度在喷嘴入口处指定，压力为5.5MPa，温度为3327K，环境温度为300K；考虑高温燃气与喷管内壁的热传导和对流传热，建立能量守恒方程计算温度场：这里，T是温度，Cp是热容，k是热导率；

为了更好地计算跨声速流动及可能存在的旋涡，湍流模型采用realizablek‑ε两方程模型,近壁面处理为标准壁面函数；建立质量守恒和动量守恒方程计算流场：这里，ρ是来自理想气体状态方程的气体混合物的质量密度，u是质量平均速度的矢量，p是压力；I是单位张量；在入口设置压力，出口设置流速边界；其他壁面设为无滑移边界；

利用COMSOL Multiphysics软件进行仿真分析，获得喷管结构的温度场、流场分布数据；计算得到C/C复合材料喷管的温度场分布云图；得到喷管壁与气流交界面温度分布情况；得到喷管燃气的速度、压力轴向分布云图和线图；

步骤3：建立Al2O3粒子在喷管中运动的本构模型；

建立流体流动颗粒追踪模型，计算中不考虑粒子之间的相互作用，粒子在气相中运动时，遵循牛顿第二运动定律，对于一个在流体中下沉的粒子，总力 是重力Fg和曳力FD的总和，ρp是粒子的密度，ρ是周围流体的密度，g是重力引起的加速度(海平面以下约为9.8m/s2)，μ是周围流体的动力粘度，dp是粒径，是周围流体的速度，是粒子的速度；

重力表达式中的(ρp‑ρ)/ρp项代表浮力；曳力表达式来自斯托克斯曳力定律(stokes drag law)；当粒子的相对雷诺数非常小时，此曳力定律较为适用：假设粒子没有改变大小，dp和mp为常数，则球形粒子的质量是，因此，得到粒子运动方程的简化表达式：

这里，引入了常数τp，τp具有时间单位，并且通常被称为拉格朗日时间尺度或者粒子速度响应时间，利用COMSOL Multiphysics软件进行仿真分析，获得Al2O3粒子在喷管中的轨迹分布图和Al2O3粒子在喷管各位置区间的滞留率，并根据滞留率判断Al2O3粒子对喷管各位置的机械侵蚀程度；

步骤4：建立微观尺度C/C复合材料的几何模型；

建立半径为3.5μm的单根纤维、包裹在纤维外围的厚度为0.5μm的基体层，及覆盖在材料表面的高温气流组成的C/C复合材料微观模型；

步骤5：建立相关化学反应、稀物质传递、水平集模型的本构方程；

n

将C与H2O、H、CO2的三个反应输入，并且输入相应的反应频率因子A，反应温度指数T ，反n应活化能指数E，根据三参数Arrhenius公式k＝A*T*exp(‑E/RT)来确定反应速率；

建立稀物质传递模型：采用菲克扩散模型和附加对流传递机理；

这里，ci是物质i相对于固定坐标系的摩尔浓度，Ji为扩散通量；

步骤6：将稀物质传递模型中获得的气体浓度与反应速率耦合起来，并采用水平集方法追踪C/C复合材料界面的退移；

设在高温气流中，φ＝0，在复合材料中，φ＝1；因此，水平集函数可视为复合材料的体积分数；分隔两个相态的流体界面的传递由下式给定ε参数决定界面的厚度；当对水平集方程使用数值稳定方法时，通常界面厚度指定为ε＝hc/2，其中hc是界面经过的区域中的特征网格尺寸；γ参数决定重新初始化的数量，合适的γ值是模型中出现的最大速度大小；

利用COMSOL Multiphysics软件进行瞬态仿真分析，获取C/C复合材料的微观形貌演变过程及烧蚀率数据；

‑5

在微观模型中，将烧蚀的温度设定为3000K，此时计算出气体的扩散系数为7.8×10 ，‑3基体的反应速率常数Km＝6.02×10 m/s，在该条件下先运用稀物质传递计算出几何模型各区域反应物的浓度，再将浓度与反应速率、材料摩尔体积耦合去计算材料的烧蚀形貌；

特定时间的选择基于以往文献提出的参考时间公式；

τo＝rf/(C0υskf)

这里，rf为纤维的半径，C0为气体的摩尔浓度，vs为材料的摩尔体积，kf为反应速率；

步骤7：建立气流剪切纤维的模型几何；

利用COMSOL Multiphysics软件进行仿真分析，获得气流剪切纤维尖端模型的流速分布和纤维应力分布结果，对气流与粒子对纤维尖端进行受力分析比较，计算粒子对材料造成的机械剥蚀；

得出粒子和气流对纤维的法向作用力Fn，如下所示：

θ为撞击角度，Δt为撞击作用时间，m，v，ρ，V代表质量、速率、密度和体积；ρp，vp分别为粒子的密度和速率，ρg，vg分别为气流的密度和速率；

步骤8：将微观模型中的化学反应速率转化为宏观模型整体的平均反应速率；

利用COMSOL Multiphysics软件进行仿真分析，获得各气体组分的浓度分布；在宏观模型的基础上，考虑氧化组分和C/C复合材料表面发生异质化学反应，加入稀物质传递模型来表征氧化反应中氧化物气体组分和生成物气体组分的浓度变化，计算中忽略氧化铝颗粒的机械剥蚀，在此模型中不考虑喷管几何形貌的变化。

说明书

固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法

技术领域

[0001] 本发明涉及固体火箭发动机技术领域，特别涉及一种宏观或微观尺度的固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法。

背景技术

[0002] C/C复合材料是一种碳纤维增强碳基体的复合材料，具有高强度(尤其是高温强度稳定)、抗热冲击性能好、耐烧蚀性好、耐含固体微粒燃气的冲刷、热膨胀系数小、导热率较低等一系列的优异性能，是目前唯一可用于3000℃以上的高温热防护材料，在火箭发动机特别是固体火箭发动机喷管喉衬已获得广泛应用，但较弱的抗氧化性限制了其高温使用寿命。

[0003] 近年来，国内外采用Hf、Ta、Si等难熔金属化合物对C/C复合材料进行基体改性，以获得性能更优异的热防护材料。但由于对C/C复合材料烧蚀行为、结构演化和烧蚀机理缺乏系统深入的研究，基体改性效果往往不甚理想。因此，研究C/C复合材料在真实火箭发动机燃气环境下的烧蚀行为和烧蚀机理，对正确选材和材料改性都具有重要指导意义。

[0004] 本发明选取C/C复合材料，研究固体火箭发动机在典型的工作状态中喷管喉衬材料的烧蚀行为和机理。碳基材料的烧蚀主要由炭相的热化学烧蚀和机械侵蚀两部分组成。热化学烧蚀指炭相与燃气中的H2O、CO2等组分反应引起的氧化烧蚀，机械侵蚀指在气流压力、剪切力、粒子等机械载荷作用下因材料本身密度不均匀，造成烧蚀差异而引起的颗粒状剥落或因热应力破坏引起的片状剥落。

[0005] 目前，在固体火箭发动机工作过程中，燃烧和高温、高压燃气的流动都是十分恶劣的工作环境，受限于材料和实验，多次实验获取充足的数据需付出高昂的经济和时间成本。而且，现有的绝大多数流场测试实验获取的数据十分有限，不足以充分分析烧蚀行为和机理。与之相比，基于合理假设的数值仿真可以提供大量的工作过程信息，为烧蚀分析提供佐证。但是，由于机理认识不足和准确模型数据的缺乏，模拟结果与真实情况仍有相当的差距。

[0006] 与本发明相关的现有技术一；

[0007] 电弧风洞是固体火箭发动机用C/C复合材料烧蚀实验较常用的方法之一，主要是利用大功率的电弧加热器注入能量，快速加热流动的空气，提高来流总温，从而进行气动热烧蚀试验。目前已有研究人员利用电弧风洞对石墨材料在空气流中的烧蚀行为进行了研究，对石墨在不同压力(0.3～4.4atm)和不同温度(2570～4030K)下的线烧蚀率和质量烧蚀率进行了测量，并将实验值和利用几种平衡热化学烧蚀理论的计算值进行了比较。研究结果表明：(1)在所有的平衡热化学烧蚀理论中，除在物质的热力学性质JANAF数据表基础上建立的热化学烧蚀理论外，其余热化学烧蚀理论因过高估计了材料的质量烧蚀率而失效；(2)建立在JANAF数据表基础上的热化学烧蚀理论，在氧化扩散控制温度范围内，材料质量烧蚀率的理论计算值与实验值相吻合，但在更高的温度范围内，理论计算值低于实验值。

[0008] 现有技术一的缺点

[0009] 地面实验研制耗费巨大，周期很长。尤其是大型固体火箭发动机设计，每失败一次就会造成巨大的经济损失。且与传统的地面实验相比，电弧风洞实验需要优先调节恢复焓，导致冷壁热流存在最低下限，特别是对于“高焓值，低热流”的情况难以模拟。再者燃烧和高温、高压燃气的流动对于目前现有的绝大多数流场测试手段而言都是十分恶劣的工作环境，因此从一次实验中所能测量的实验数据是极为有限的，不足以用来分析实验现象，而多次实验成本十分高昂。

[0010] 与本发明相关的现有技术二；

[0011] 随着计算机技术的迅速发展，通过模拟了解发动机内部燃气流场的流动本质和传递机理已经成为可能。近年来，有研究人员采用FLUENT公司成熟的前处理软件GAMBIT和通用计算流体动力学(CFD)软件FLUENT按照稳态轴对称问题求解了单纯气相流动问题，计算了实验过程中火箭发动机燃气由燃烧室流出、经喷管排出到外界大气环境整个流场参数的分布情况。

[0012] 现有技术二的缺点

[0013] 该技术没有考虑到多相流，单纯将实验简化为纯气相流动，这与实际并不相符，因为实际情况中氧化铝颗粒对于材料表面的剥蚀影响非常大。另外该技术的流场计算没有考虑喷管材料烧蚀，但实际火箭工作过程中，烧蚀引起的喷管壁面退移会改变喷管的气动几何和表面状态，进而影响到近壁面流场参数和热量、质量的传递过程。随之而来，燃气参数的变化又会强烈地反作用于喷管材料，对烧蚀产生影响。因此，该技术的计算结果虽然为烧蚀机理分析提供了一定基础，但是并不能代表实际情况。

[0014] 缩略语和关键术语定义

[0015] SRM：固体火箭发动机(Solid Rocket Motor)；

[0016] AP/Al/HTPB：复合推进剂(Ammonium Perchlorate/Aluminum powder/Hydroxyl Terminated Polybutadiene)；

[0017] C/C复合材料：C/C复合材料是一种碳纤维增强碳基体的复合材料；

[0018] 数值模拟：数值模拟也叫计算机模拟。依靠电子计算机，结合有限元或有限容积的概念，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类问题研究的目的；

[0019] 水平集：用于跟踪流体流动模型中的移动界面，求解水平集函数。

发明内容

[0020] 本发明针对现有技术的缺陷，提供了一种固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法。

[0021] 为了实现以上发明目的，本发明采取的技术方案如下：

[0022] 一种固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法,包括以下步骤：

[0023] 步骤1：建立C/C复合材料的喷管几何模型；

[0024] 采用通道截面形状先收敛后扩张的拉瓦尔喷管，基于喷管及复合材料结构的对称性，建立二维轴对称几何模型，从下到上依次是：收敛段、喉部、扩张段，外端为C/C复合材料。

[0025] 步骤2：建立气相控制方程；

[0026] 设燃料气体中的冷凝相颗粒的质量含量为30wt％,质量流率为0.01kg/s。燃烧室压力和温度在喷嘴入口处指定，压力为5.5MPa，温度为3327K，环境温度为300K。考虑高温燃气与喷管内壁的热传导和对流传热，建立能量守恒方程计算温度场：

[0027]

[0028] 这里，T是温度，Cp是热容，k是热导率。

[0029] 为了更好地计算跨声速流动及可能存在的旋涡，湍流模型采用realizableκ‑ε两方程模型,近壁面处理为标准壁面函数。建立质量守恒和动量守恒方程计算流场：

[0030]

[0031]

[0032] 这里，ρ是来自理想气体状态方程的气体混合物的质量密度，u是质量平均速度的矢量，p是压力。I是单位张量。在入口设置压力，出口设置流速边界。其他壁面设为无滑移边界。

[0033] 利用COMSOL Multiphysics软件进行仿真分析，获得喷管结构的温度场、流场分布数据。计算得到C/C复合材料喷管的温度场分布云图。得到喷管壁与气流交界面温度分布情况。得到喷管燃气的速度、压力轴向分布云图和线图。

[0034] 步骤3：建立Al2O3粒子在喷管中运动的本构模型；

[0035] 建立流体流动颗粒追踪模型，计算中不考虑粒子之间的相互作用，粒子在气相中运动时，遵循牛顿第二运动定律，

[0036]

[0037] 对于一个在流体中下沉的粒子，总力 是重力Fg和曳力FD的总和，

[0038]

[0039]

[0040] ρp是粒子的密度,ρ是周围流体的密度,g是重力引起的加速度(海平面以下约为2
9.8m/s),μ是周围流体的动力粘度,dp是粒径, 是周围流体的速度, 是粒子的速度。

[0041] 重力表达式中的(ρp‑ρ)/ρp项代表浮力。曳力表达式来自斯托克斯曳力定律(stokes drag law)。当粒子的相对雷诺数非常小时，此曳力定律较为适用：

[0042]

[0043] 假设粒子没有改变大小，dp和mp为常数，则球形粒子的质量是，

[0044]

[0045] 因此，得到粒子运动方程的简化表达式：

[0046]

[0047] 这里，引入了常数τp，τp具有时间单位，并且通常被称为拉格朗日时间尺度或者粒子速度响应时间，

[0048]

[0049] 利用COMSOL Multiphysics软件进行仿真分析，获得Al2O3粒子在喷管中的轨迹分布图和Al2O3粒子在喷管各位置区间的滞留率，并根据滞留率判断Al2O3粒子对喷管各位置的机械侵蚀程度。

[0050] 步骤4：建立微观尺度C/C复合材料的几何模型；

[0051] 建立半径为3.5μm的单根纤维、包裹在纤维外围的厚度为0.5μm的基体层，及覆盖在材料表面的高温气流组成的C/C复合材料微观模型。

[0052] 步骤5：建立相关化学反应、稀物质传递、水平集模型的本构方程；

[0053] 将C与H2O、H、CO2的三个反应输入，并且输入相应的反应频率因子A，反应温度指数n nT，反应活化能指数E，根据三参数Arrhenius公式k＝A*T*exp(‑E/RT)来确定反应速率。

[0054] 建立稀物质传递模型：采用菲克扩散模型和附加对流传递机理。

[0055]

[0056]

[0057] 这里,ci是物质i相对于固定坐标系的摩尔浓度，Ji为扩散通量。

[0058] 步骤6：将稀物质传递模型中获得的气体浓度与反应速率耦合起来，并采用水平集方法追踪C/C复合材料界面的退移。

[0059] 设在高温气流中，φ＝0，在复合材料中，φ＝1。因此，水平集函数可视为复合材料的体积分数。分隔两个相态的流体界面的传递由下式给定

[0060]

[0061] ε参数决定界面的厚度。当对水平集方程使用数值稳定方法时，通常界面厚度指定为ε＝hc/2，其中hc是界面经过的区域中的特征网格尺寸。γ参数决定重新初始化的数量，合适的γ值是模型中出现的最大速度大小。

[0062] 利用COMSOL Multiphysics软件进行瞬态仿真分析，获取C/C复合材料的微观形貌演变过程及烧蚀率数据。

[0063] 在微观模型中，将烧蚀的温度设定为3000K，此时计算出气体的扩散系数为7.8×‑5 ‑310 ，基体的反应速率常数Km＝6.02×10 m/s，在该条件下先运用稀物质传递计算出几何模型各区域反应物的浓度，再将浓度与反应速率、材料摩尔体积耦合去计算材料的烧蚀形貌。

[0064] 特定时间的选择基于以往文献提出的参考时间公式；

[0065] τ0＝rf/(C0vskf)

[0066] 这里，rf为纤维的半径，C0为气体的摩尔浓度，vs为材料的摩尔体积，kf为反应速率。

[0067] 步骤7：建立气流剪切纤维的模型几何；

[0068] 利用COMSOL Multiphysics软件进行仿真分析，获得气流剪切纤维尖端模型的流速分布和纤维应力分布结果，对气流与粒子对纤维尖端进行受力分析比较，计算粒子对材料造成的机械剥蚀。

[0069] 得出粒子和气流对纤维的法向作用力Fn，如下所示：

[0070]

[0071]

[0072] θ为撞击角度，Δt为撞击作用时间，m，v，ρ，V代表质量、速率、密度和体积。ρp，vp分别为粒子的密度和速率，ρg，vg分别为气流的密度和速率。

[0073] 步骤8：将微观模型中的化学反应速率转化为宏观模型整体的平均反应速率；

[0074] 利用COMSOL Multiphysics软件进行仿真分析，获得各气体组分的浓度分布。在宏观模型的基础上，考虑氧化组分和C/C复合材料表面发生异质化学反应，加入稀物质传递模型来表征氧化反应中氧化物气体组分和生成物气体组分的浓度变化，计算中忽略氧化铝颗粒的机械剥蚀，在此模型中不考虑喷管几何形貌的变化。

[0075] 与现有技术相比，本发明的优点在于：

[0076] 更加详细具体的将湍流，流体传热，粒子追踪，稀物质传递和水平集进行多场耦合，综合且全面地考虑了C/C复合材料喷管在宏观、微观尺度的烧蚀行为和机理，并将两种尺度建立了一定联系，获得了一系列与实际相接近的结果，为固体火箭发动机在典型的工作状态中喷管喉衬材料的烧蚀分析提供佐证，对喷管的热性能研究和精细化设计有重要参考价值。

附图说明

[0077] 图1是本发明实施例固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法流程图；

[0078] 图2是本发明实施例C/C喷管示意图；

[0079] 图3是本发明实施例喷管中温度场分布图；

[0080] 图4是本发明实施例喷管壁与气流交界面温度分布图；

[0081] 图5是本发明实施例喷管燃气的速度、压力轴向分布云图；

[0082] 图6是本发明实施例喷管燃气的速度、压力轴向分布线图；

[0083] 图7是本发明实施例喷管中Al2O3粒子的轨迹分布图；

[0084] 图8是本发明实施例喷管中Al2O3粒子的平均滞留率分布图；

[0085] 图9是本发明实施例C/C复合材料喷管实验后宏观照片图；

[0086] 图10是本发明实施例微观尺度C/C复合材料烧蚀模型图；

[0087] 图11是本发明实施例a2、a3、a4与a1反应速率在不同温度下的比值曲线图；

[0088] 图12是本发明实施例不同时间材料形貌演变过程图；

[0089] 图13是本发明实施例C/C复合材料喷管表面烧蚀形貌图；

[0090] 图14是本发明实施例C/C复合材料在不同时间节点的线烧蚀率图；

[0091] 图15是本发明实施例气流剪切纤维模型几何图；

[0092] 图16是本发明实施例气流的流速分布和纤维的应力分布图；

[0093] 图17是本发明实施例粒子和气流对复合材料的受力分析图；

[0094] 图18是本发明实施例不同时刻宏观模型中的氧化反应速率拟合曲线图；

[0095] 图19是本发明实施例各气体组分的浓度分布图，其中(a)为H2O、(b)为CO2、(c)为CO、(d)为H2。

具体实施方式

[0096] 为使本发明的目的、技术方案及优点更加清楚明白，以下根据附图并列举实施例，对本发明做进一步详细说明。

[0097] 如图1所示，一种固体火箭发动机C/C复合材料喷管烧蚀行为建模仿真方法,包括以下步骤：

[0098] 步骤1：建立C/C复合材料的喷管几何模型；

[0099] 喷管是火箭发动机的一个重要的能量转换部件，它对发动机的性能有着重要的影响。为了使燃气的流动速度能够从亚声速加速到超声速，采用通道截面形状先收敛后扩张的拉瓦尔喷管。为减少计算时间及计算内存，基于喷管及复合材料结构的对称性，建立二维轴对称几何模型，如图2所示。从下到上依次是：收敛段、喉部、扩张段，外端为C/C复合材料。

[0100] 步骤2：建立气相控制方程；

[0101] 发动机中使用的推进剂是AP/Al/HTPB，其中铝粉的质量含量为17wt％。假设铝粉在进入喷嘴部时已经完全反应形成氧化铝，则燃料气体中的冷凝相颗粒的质量含量为30wt％,质量流率为0.01kg/s。燃烧室压力和温度在喷嘴入口处指定，压力为5.5MPa，温度为3327K，环境温度为300K。主要考虑高温燃气与喷管内壁的热传导和对流传热，建立能量守恒方程计算温度场：

[0102]

[0103] 这里，T是温度，Cp是热容量，k是热导率。

[0104] 为了更好地计算跨声速流动及可能存在的旋涡，湍流模型采用realizableκ‑ε两方程模型,近壁面处理为标准壁面函数。建立质量守恒和动量守恒方程计算流场：

[0105]

[0106]

[0107] 这里，ρ是来自理想气体状态方程的气体混合物的质量密度，u是质量平均速度的矢量，p是压力。I是单位张量。在入口设置压力，出口设置流速边界。其他壁面设为无滑移边界。

[0108] 利用COMSOL Multiphysics软件进行仿真分析，获得喷管结构的温度场、流场分布数据。如图3所示，计算得到C/C复合材料喷管的温度场分布云图。图4则显示了喷管壁与气流交界面温度分布。图5、6分别为喷管燃气的速度、压力轴向分布云图和线图。

[0109] 步骤3：建立Al2O3粒子在喷管中运动的本构模型；

[0110] 建立流体流动颗粒追踪模型，计算中不考虑粒子之间的相互作用，粒子在气相中运动时，遵循牛顿第二运动定律，

[0111]

[0112] 对于一个在流体中下沉的粒子，总力 是重力Fg和曳力FD的总和，

[0113]

[0114]

[0115] ρp是粒子的密度,ρ是周围流体的密度,g是重力引起的加速度(海平面以下约为2
9.8m/s),μ是周围流体的动力粘度,dp是粒径, 是周围流体的速度, 是粒子的速度。

[0116] 重力表达式中的(ρp‑ρ)/ρp项代表浮力。曳力表达式来自斯托克斯曳力定律(stokes drag law)。当粒子的相对雷诺数非常小时，此曳力定律较为适用：

[0117]

[0118] 假设粒子没有改变大小(dp和mp为常数)，则球形粒子的质量是，

[0119]

[0120] 因此，我们得到粒子运动方程的简化表达式：

[0121]

[0122] 这里，引入了常数τp，τp具有时间单位，并且通常被称为拉格朗日时间尺度或者粒子速度响应时间，

[0123]

[0124] 利用COMSOL Multiphysics软件进行仿真分析，获得Al2O3粒子在喷管中的轨迹分布图和Al2O3粒子在喷管各位置区间的滞留率，并根据滞留率判断Al2O3粒子对喷管各位置的机械侵蚀程度。图7显示的是Al2O3粒子在喷管内速度的轨迹分布图。采用不同总量的粒子数来考察粒子在喷管壁面的滞留率，表1是100～500个粒子在壁面9‑19mm位置区间的滞留结果，按照每隔1mm的间距进行统计。图8则为粒子在喷管壁面的平均滞留率曲线，可以看出，Al2O3粒子主要滞留在喷管的收敛段及收敛段与喉部的过渡区域，因此可以判断Al2O3粒子对喷管的机械剥蚀主要发生在收敛段。

[0125] 表1不同Al2O3粒子总数在各位置区间的滞留率

[0126]

[0127] 颗粒的轨迹分布与喷管烧蚀实验结果非常吻合，如图9所示，实验报道烧蚀后喷管在喉部上游沉积了一层Al2O3，在下游却没有观察到，沉积层厚度为0.2～0.4mm。据研究表明，直径在1～10μm的Al2O3粒子占总粒子质量的71％，剩余29％的粒子直径在10～100μm之间。撞击壁面的粒子大部分偏小，因为其粒径小，随流性好，而大粒径粒子由于惯性较大，往往会向喷管的中轴线靠。因此假设50％的小粒径粒子(平均直径5μm)和5％的大粒径粒子(平均直径50μm)撞击在壁面上后迅速凝固在壁面上，不随气流冲刷，那么各个位置区间的沉积厚度如表2所示，平均沉积厚度为0.52mm，而实验中收敛段沉积层厚度为0.2～0.4mm，其中实际实验包含了沉积层消融后被气流冲刷下去的情况，因此计算的平均沉积厚度0.52mm比实际值偏大是合理的。

[0128] 表2 Al2O3粒子在各位置区间的沉积厚度

[0129]

[0130] 步骤4：建立微观尺度C/C复合材料的几何模型；

[0131] 固体火箭发动机喷嘴中的高温气流由AP/Al/HTPB复合推进剂燃烧的热燃烧产物组成。燃烧产物中主要物质是H2O、CO2、CO、OH、H和H2，以及少量的HCl和N2和可忽略不计的O2和O。建立半径为3.5μm的单根纤维、包裹在纤维外围的厚度为0.5μm的基体层，及覆盖在材料表面的高温气流组成的C/C复合材料微观模型，如图10所示。

[0132] 步骤5：建立相关化学反应、稀物质传递、水平集模型的本构方程；

[0133] 在高温高压环境下，喷嘴表面容易受到推进剂燃烧产物的化学侵蚀。建立化学反应，查阅资料，发现有两组化学反应动力学数据差异较大，如下表3a、b，将本发明气体反应n输入，并且输入相应的反应频率因子(A)，反应温度指数(T)，反应活化能指数(E)根据三参n
数Arrhenius公式k＝A*T *exp(‑E/RT)来确定反应速率。我们将两组数据分别应用到模型计算中，最后将算出的烧蚀率进行了对比。

[0134] 表3气相反应模型列表

[0135]

[0136] a反应速率以阿仑尼乌斯形式表示：k＝A*Tn*e‑E/RT，表中的反应速率k是在3000K情况下所得。

[0137] 图11为a2、a3、a4与a1反应速率在不同温度下的比值曲线，从图中的比值可以看出，a4反应速率太小，因此在正式计算中忽略不计，只考虑C与H2O、H、CO2的烧蚀反应。

[0138] 建立稀物质传递模型：采用菲克扩散模型和附加对流传递机理。

[0139]

[0140]

[0141] 这里,ci是物质i相对于固定坐标系的摩尔浓度，Ji为扩散通量。

[0142] 步骤6：将稀物质传递模型中获得的气体浓度与反应速率耦合起来，并采用水平集方法追踪C/C复合材料界面的退移。

[0143] 在高温气流中，φ＝0，在复合材料中，φ＝1。因此，水平集函数可视为复合材料的体积分数。分隔两个相态的流体界面的传递由下式给定

[0144]

[0145] ε参数决定界面的厚度。当对水平集方程使用数值稳定方法时，通常界面厚度可指定为ε＝hc/2，其中hc是界面经过的区域中的特征网格尺寸。γ参数决定重新初始化的数量，合适的γ值是模型中出现的最大速度大小。

[0146] 利用COMSOL Multiphysics软件进行瞬态仿真分析，获取C/C复合材料的微观形貌演变过程及烧蚀率数据。

[0147] 在微观模型中，将烧蚀的温度设定为3000K，此时计算出气体的扩散系数为7.8×‑5 ‑310 ，基体的反应速率常数Km＝6.02×10 m/s，在该条件下先运用稀物质传递计算出几何模型各区域反应物的浓度，再将浓度与反应速率、材料摩尔体积耦合去计算材料的烧蚀形貌。
‑5 ‑3 ‑4
图12为扩散系数D＝7.8×10 ，Km＝6.02×10 m/s，Kf＝7.08×10 m/s，舍伍德数Sh＝2.7×‑4
10 时计算出来的C/C复合材料的微观形貌演变过程。由于材料的形貌演变过程十分缓慢，只保存形貌演变过程中特定时间的数据。特定时间的选择基于以往文献提出的参考时间公式；

[0148] τ0＝rf/(C0vskf)

[0149] 这里，rf为纤维的半径，C0为气体的摩尔浓度，vs为材料的摩尔体积，kf为反应速率。

[0150] 从图12中可以看到，C/C复合材料经过烧蚀最终形成了具有尖角的纤维微观烧蚀形貌，这与实际烧蚀实验中观察到的表面形貌一致(图13)。烧蚀过程中基体被大量消耗，这是由于C/C复合材料基体的密度小于纤维的密度，基体的烧蚀速率大于纤维烧蚀速率，于是基体退移速度比纤维束快，随着烧蚀的进行，纤维边缘的缝隙逐渐加深，大部分纤维变细变尖，部分纤维表面存在小的缺口，可能是局部杂质或缺陷产生的优先烧蚀。

[0151] 图14为C/C复合材料在不同时间节点的线烧蚀率，a和b是用表3中两组不同的动力‑3 ‑3学数据计算的结果。a的纤维、基体平均线烧蚀率分别为0.337×10 mm/s、0.912×10 mm/‑3 ‑3
s，b的纤维、基体平均线烧蚀率分别为10.27×10 mm/s、30.17×10 mm/s。实验报道中的线‑3 ‑3
烧蚀率为2.2×10 mm/s～4.6×10 mm/s，该实验数据处在这两组数据的中间，说明计算结果存在一定的误差，但是在合理的范围之内。从数据中也可看出，两组基体的线烧蚀率值都约为纤维的3倍，这进一步验证了微观形貌的变化。

[0152] 步骤7：建立气流剪切纤维的模型几何；

[0153] 由微观模型形貌变化可知，烧蚀过程中基体退移纤维逐渐变尖。变尖的纤维头部抗屈应力下降，Al2O3粒子撞击到材料表面后可能会打断烧尖的纤维，之后呈现出扁平的头部。基于层流和固体力学理论建立流固耦合模型，研究高温气流和Al2O3粒子对于纤维尖端的剪切作用，如图15所示，烧尖的纤维几何从微观烧蚀模型中提取，高温气流域的高度由宏观模型中边界层(收敛段与喉部的过渡区域)的厚度决定，为0.2mm。

[0154] 利用COMSOL Multiphysics软件进行仿真分析，获得气流剪切纤维尖端模型的流速分布和纤维应力分布结果，对气流与粒子对纤维尖端进行受力分析比较，计算粒子对材料造成的机械剥蚀。

[0155] 图16为气流剪切纤维尖端模型的流速分布和纤维应力分布结果，从应力分布云图可看到，纤维剪切应力最大值为5.22MPa，处在纤维底部边缘区域，剪切应力最小为1.45MPa，在纤维中部区域。

[0156] 据文献报道，C/C复合材料的剪切强度约为17.6MPa，也就是说气流的剪切力不足以切断纤维头部。但是根据理论分析，粒子对纤维的作用力远大于气流的作用力，图17粒子和气流对复合材料的受力分析。研究认为撞击力起主要破坏作用(法向作用力)，在此忽略切向作用力和旋转效应，假设粒子和气流接触材料后均不反弹。根据动量定理和牛顿第三定律，可以得出粒子和气流对纤维的法向作用力Fn，如下所示：

[0157]

[0158]

[0159] 其中θ为撞击角度，Δt为撞击作用时间，m，v，ρ，V代表质量、速率、密度和体积。ρp，vp分别为粒子的密度和速率，ρg，vg分别为气流的密度和速率。

[0160] 假设α为单位体积粒子对纤维的作用力Fnp和气流对纤维的作用力Fng之比，则由气3 3
流密度6.3833kg/m、速度900m/s，Al2O3粒子密度3850kg/m、速度332m/s，计算得α＝222.5。
由此可知，粒子对纤维的作用力远大于气流的作用力，足以打断烧尖的纤维。表4为粒子撞击引起的机械剥蚀数据，其中s代表的是粒径1～10μm的小粒径，这些小粒子主要是会对纤维尖端造成影响，b代表的是10～100μm的大粒子，其质量和惯性都比较大，通常会在大的惯性作用下直接撞击到材料表面产生机械剥蚀造成材料动态断裂，炭纤维和基体被压溃，同时被气流剥离，引起较大的质量损失。从两组不同粒径的机械剥蚀数据可以看出，对C/C复合材料的机械剥蚀起主要作用的是粒径较大的粒子。

[0161] 表4粒子撞击引起的机械剥蚀

[0162]

[0163]

[0164] 步骤8：将微观模型中的化学反应速率转化为宏观模型整体的平均反应速率；

[0165] 微观模型中由于氧化烧蚀形貌起伏大，而宏观模型中喷管几何是平滑的壁面，故不可直接将微观模型中的化学反应速率应用在宏观模型中，需要经过一定的计算转化。表5中k＝Sf*kf+Sm*km为换算后不同时间节点的宏观模型中的反应速率，Sf、Sm分别代表的是微观烧蚀中纤维和基体形貌变化的影响因子，是由纤维和基体在不同时刻的烧蚀长度变换得来，kf、km分别代表的是微观烧蚀中纤维和基体的反应速率。图18为不同时刻宏观模型中的氧化反应速率拟合曲线，拟合公式为 从图18可以看出，随着时间的增长，材料的氧化反应速率逐渐增大，最后在慢慢接近最大值时趋向稳定。

[0166] 表5宏观模型中的平均反应速率

[0167]

[0168] 利用COMSOL Multiphysics软件进行仿真分析，获得各气体组分的浓度分布。在宏观模型的基础上，考虑氧化组分和C/C复合材料表面发生异质化学反应，加入稀物质传递模型来表征氧化反应中氧化物气体组分和生成物气体组分的浓度变化，在此计算中忽略氧化铝颗粒的机械剥蚀。由于微观上的烧蚀形貌在宏观上显得微不足道，因此在此模型中并不考虑喷管几何形貌的变化。

[0169] 图19显示了各气体组分的浓度变化，从分布结果来看，随着反应的进行，氧化物质H2O、CO2被消耗，浓度沿着气流与复合材料交界面逐渐减小，CO、H2的浓度则沿交界面逐渐增大。

[0170] 本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。