一种固体推进剂烟雾影响可见光制导信号虚拟试验方法

摘要

本发明公开了一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，包括以下步骤：一、羽流场计算及数据存储；二、气相产物羽流场数据结构网格化处理；三、凝聚相产物羽流场数据结构网格化处理；四、可见光制导信号影响虚拟测试：初始参数设定、制导信号入射方向确定、粒子轨迹途经制导信号影响网格确定、粒子轨迹网格化数据计算、各粒子轨迹衰减系数换算和入射方向上可见光制导信号衰减系数计算。本发明方法步骤简单、设计合理且实现方便、使用效果好，能简便、快速且准确获得固体推进剂烟雾对可见光制导信号的影响结果。

说明书

技术领域

本发明属于固体推进剂烟雾光学制导信号影响研究技术领域，尤其是 涉及一种固体推进剂烟雾影响可见光制导信号虚拟试验方法。

背景技术

随着机械化、信息化技术的发展，现代武器系统越来越多的采用雷达、 光电等自动跟踪、控制和制导手段，武器系统的自动化程度、打击精度等 方面得到显著提高，极大程度上提升了弹箭武器系统的作战效能。但与此 同时，弹箭武器发射和运载过程中的烟焰问题危害性日益显现，对武器系 统信息化、智能化、隐身化性能影响很大，对武器系统快速、精确打击、 隐身等功能造成了严重的不利影响，甚至已成为制约我军弹箭武器发挥正 常效能、新装备订货以及新型武器装备发展的主要技术瓶颈。尤其是近年 来，重型反坦克导弹研制过程中出现推进剂装药烟雾遮挡制导信号的问 题，解决问题耗费巨大。制导信号通常采用红外、紫外、激光、可见光等。 现如今，可见光系统已开始大量应用于制导武器系统，其显著提高了武器 系统的自动化程度和打击精度，极大程度上提升了弹箭武器系统的作战效 能。

固体推进剂是弹箭武器系统的发射能源，其低特征信号性能是保证武 器系统战场使命完成的重要因素。羽流是从火箭发动机喷管喷射出来的羽 毛状的高速高温燃气流，火箭羽流是一种气体分子浓度大、电子密度和电 子碰撞频率都很高的弱等离子体，其与雷达微波之间会发生相互作用，使 微波信号功率大为降低，影响导弹的制导。而羽流场是火箭发动机喷管喷 射出来的羽毛状的高速高温燃气流的流体运动所占据的空间。目前，对推 进剂燃烧和燃烧产物流动的数学模型进行模拟试验时，通常都采用FLUNT 软件对固体推进剂燃烧后的羽流场进行模拟，且采用FLUNT软件完成所试 验固体推进剂的羽流场计算后，便能自动输出*.out文件供微波衰减计算、 红外辐射计算以及光电特性计算使用。火箭排气羽流的特征信号（Exhaust  Plume Signature）是一种包含有系统或火箭发动机排气全部性能或特性 的术语，此性能或特性可被用作探测、识别、或拦截执行任务的发射平台 或导弹。羽流的特征信号特性主要包括烟、辐射能的散发、能见度（视程 /能见距离）和雷达波吸收等4项参数。固体推进剂烟雾是羽流特征信号 的重要表征参数，其从表观特征可以细分为：一次烟雾和二次烟雾，一次 烟雾主要由火箭发动机喷管喷出的凝相颗粒构成，主要存在形式为金属氢 氧化物、氧化物和氯化物等，复合推进剂中的金属粉（铝、镁等）燃烧产 生的金属氧化物（Al₂O₃、MgO等）是一次烟雾的主要来源；二次烟雾由燃 气中可凝气体组成，是燃气与大气相互作用的结果。复合推进剂中的AP 和含氟组分燃烧分解后形成HCl、HF，在特定的温度和湿度条件下，与空 气相结合形成水和HCl、HF的共沸液滴烟云。

但现如今在低特征信号推进剂研制及评价过程中，解决问题需要耗费 大量的人力物力，研制周期大幅度延长，因而带来巨大的经济损失和政治 影响，现如今迫切需要虚拟试验对固体推进剂相关性能进行快速评价。而 虚拟试验是借助于计算机的高速解算，按实际试验要求对基于固体推进剂 燃烧、燃烧产物流动的数学模型以及可见光制导信号穿越固体推进剂羽流 进行计算机建模与求解，不仅可以作为真实试验的前期准备工作或在一定 程度上替代传统的试验（如一些极限工况）；并且能大幅减少真实试验次 数，降低试验费用，缩短试验周期；同时，具有较好的交互性，使得各种 试验信息及时反馈；另外，虚拟试验不受气象条件、场地、时间和次数的 限制，试验过程可方便实现回放、再现和重复。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一 种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其方法步骤简单、 设计合理且实现方便、使用效果好，能简便、快速且准确获得固体推进剂 烟雾对可见光制导信号的影响结果。

为解决上述技术问题，本发明采用的技术方案是：一种固体推进剂烟雾 影响可见光制导信号虚拟试验方法，其特征在于该方法包括以下步骤：

步骤一、羽流场计算及数据存储：采用数据处理器调用FULENT软件 对所试验固体推进剂进行羽流场计算后，将羽流场计算结果自动输出并存 储至与所述数据处理器相接的数据存储单元内；其中，所输出的羽流场计 算结果中包括气相产物羽流场计算结果和凝聚相产物羽流场计算结果；

步骤二、气相产物羽流场数据结构网格化处理，其处理过程如下：

步骤201、喷流区域气相产物流场数据读取：采用所述数据处理器从 步骤一中所输出的气相产物羽流场计算结果中，读取喷流区域内所有非结 构网格节点的气相产物流场数据；所述喷流区域为步骤一中羽流场计算时 所用发动机喷管出口后方的矩形区域；

步骤202、轴向坐标轴上非结构网格点提取：采用所述数据处理器从 步骤201中所述喷流区域内的所有非结构网格节点中，提取位于轴向坐标 轴上的所有非结构网格点，本步骤中所提取的位于轴向坐标轴上的非结构 网格点总数量为N_X；其中，轴向坐标轴为所述发动机喷管的中心轴线所 在的横坐标轴，位于轴向坐标轴上的非结构网格点的径向坐标y_h=0且其轴 向坐标x_h≥0，其中h为正整数，且h=1，2，…，N_X；

步骤203、径向坐标轴上非结构网格点提取：采用所述数据处理器从 步骤201中所述喷流区域内的所有非结构网格节点中，提取位于径向坐标 轴上的所有非结构网格点，且所提取的位于径向坐标轴上的非结构网格点 总数量为N_Y出；其中，径向坐标轴为所述发动机喷管出口处所在的纵坐标 轴且发动机喷管出口处的轴向坐标值为0，位于径向坐标轴上的非结构网 格点的轴向坐标x_k1=0且其径向坐标y_k1≥0，其中k1为正整数，且k1=1， 2，…，N_Y出；

步骤204、构建结构网格图：将N_X条直线x=x_h和N_Y出条直线y=y_k1正 交后，构建出一个包含（N_X－1）×（N_Y出－1）个矩形网格的结构网格图；

步骤205、气相产物羽流场数据结构网格化处理：采用所述数据处理 器对步骤204中所构建结构网格图中各矩形网格的四个顶点的气相产物流 场数据分别进行重新赋值；所有矩形网格中各顶点的气相产物流场数据重 新赋值方法均相同，其中对于所构建结构网格图中任一个矩形网格的任一 顶点的气相产物流场数据进行重新赋值时，所述数据处理器先自步骤201 中所述喷流区域内的所有非结构网格节点中找出与当前被赋值顶点距离 最近的非结构网格节点，并将所找出的非结构网格节点的气相产物流场数 据赋给当前被赋值的顶点；

步骤三、凝聚相产物羽流场数据结构网格化处理：所述数据处理器利 用步骤204中所构建的结构网格图对M个不同粒径粒子的粒子轨迹数据分 别进行结构网格化处理，过程如下：

步骤301、结构网格化处理初始参数设定：采用所述参数输入单元对 M的取值和M个不同粒径粒子的粒径D_nr分别进行设定；其中，r为正整数， 且r=1,2，…，M；

步骤302、粒子轨迹数据读取：采用所述数据处理器从步骤一中所输 出的凝聚相产物羽流场计算结果中读取所设计固体推进剂的所有粒子轨 迹数据；其中，所读取的凝聚相产物羽流场计算结果包括粒子直径随轨迹 变化的文件和粒子轨迹的时间步长文件；

步骤303、发动机喷管入口矩形网格数量获取及各喷管入口矩形网格 的上下边界确定：首先，采用所述数据处理器从步骤一中输出的气相产物 羽流场计算结果中读取发动机喷管区域内所有非结构网格节点的流场数 据；之后，采用所述数据处理器从所读取的发动机喷管区域内的所有非结 构网格节点中提取位于直线x=-Δd上的所有非结构网格节点，所提取的位 于直线x=-Δd上的非结构网格点总数量为N_Y入；其中，位于直线x=-Δd上 的非结构网格节点的轴向坐标x_k2=-Δd且其径向坐标y_k2≥0，其中k2为正 整数，且k2=1，2，…，N_Y入；Δd为所述发动机喷管入口至喷管出口之间 的距离；步骤204中所构建结构网格图中位于所述发动机喷管入口处的矩 形网格为喷管入口矩形网格，所述喷管入口矩形网格数量为（N_Y入-1）个， 各喷管入口矩形网格的上下边界分别为上下相邻的两条直线y=y_k2；

步骤四、可见光制导信号影响虚拟测试，其测试过程如下：

步骤401、初始参数设定：采用所述参数输入单元输入试验用可见光 的波长λ；

步骤402、制导信号入射方向确定：先采用所述参数输入单元对可见 光的入射位置进行设定，再根据所设定的入射位置对可见光的入射方向 x=x_入进行确定，并采用所述数据处理器找出位于直线x=x_入上的所有矩形 网格；其中，x_入为所设定入射位置与所述发动机喷管出口之间的间距；位 于直线x=x_入上的所有矩形网格均为制导信号影响网格；

步骤403、粒子轨迹途经制导信号影响网格确定：采用所述数据处理 器对以步骤303中所述的（N_Y入-1）个喷管入口矩形网格作为起始网格的 所有粒子轨迹途经的制导信号影响网格进行确定；其中，以（N_Y入-1）个 喷管入口矩形网格中任一个喷管入口矩形网格作为起始网格的粒子轨迹 轨迹均包括M个不同粒径粒子轨迹，以（N_Y入-1）个喷管入口矩形网格作 为起始网格的所有粒子轨迹数量为M×（N_Y入-1）个；

步骤404、粒子轨迹网格化数据计算：采用所述数据处理器对M× （N_Y入-1）个粒子轨迹的网格化数据分别进行计算，所有粒子轨迹的网格 化数据计算方法均相同；其中，对M×（N_Y入-1）个粒子轨迹中任一条粒 子轨迹的网格化数据进行计算时，采用所述数据处理器计算得出当前所计 算粒子轨迹在所途经制导信号影响网格内的粒子数密度和粒子平均直径；

步骤405、各粒子轨迹衰减系数换算：采用所述数据处理器对M× （N_Y入-1）个粒子轨迹的衰减系数分别进行换算，所有粒子轨迹的衰减系 数换算方法均相同；其中，对M×（N_Y入-1）个粒子轨迹中任一条粒子轨 迹的衰减系数进行换算时，所述数据处理器根据公式κ₀=πr²×Q×n_p进行换 算，式中κ₀为当前所换算粒子轨迹的衰减系数，d为步骤404中计 算得出的当前所换算粒子轨迹在所途经制导信号影响网格内的粒子平均 直径，n_p为步骤404中计算得出的当前所换算粒子轨迹在所途经制导信号 影响网格内的粒子数密度；Q为衰减效率且式中 a_n和b_n为入射波长为λ时半径为r的球形粒子的Mie散射系数， Re(a_n+b_n)表示取a_n+b_n的实部，n为正整数；

步骤406、入射方向上可见光制导信号衰减系数计算：采用所述数据 处理器将步骤405中计算得出的M×（N_Y入-1）个粒子轨迹的衰减系数叠 加后，获得入射方向上可见光制导信号的衰减系数κ。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤403中进行初始参数设定时，还需采用所述参数输入单元输入试 验用可见光的入射光程长L，其中L=2×y_m，y_m为喷流径向长度且其为发 动机喷管轴线到喷流径向外边界的长度；步骤406中获得入射方向上可见 光制导信号的衰减系数κ后，所述数据处理器还需根据公式T=e^-κ×L计算得 出试验用可见光的光谱透射比T。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤404中对M×（N_Y入-1）个粒子轨迹中任一条粒子轨迹的网格化 数据进行计算时，过程如下：

步骤4041、发动机喷管入口起始网格确定及粒子质量流率计算： 采用所述数据处理器从步骤302中所读取的所有粒子轨迹数据中找出当前 所计算粒子轨迹的粒子轨迹数据，并根据所找出的粒子轨迹数据中的发动 机喷管入口处轨迹点的径向坐标和步骤303中所确定的各喷管入口矩形网 格的上下边界，对当前所计算粒子轨迹的喷管入口起始网格进行确定；

待当前所计算粒子轨迹的喷管入口起始网格确定后，所述数据处理器 找出所确定喷管入口起始网格的顶点一和顶点二的气相产物流场数据，并 从顶点一的气相产物流场数据中找出顶点一处的气相密度ρ_g1和气相轴向 速度u_g1，且从顶点二的气相产物流场数据中找出顶点二处的气相密度ρ_g2和气相轴向速度u_g2，其中顶点一为当前所计算粒子轨迹的喷管入口起始网 格的左上方顶点，且顶点二为当前所计算粒子轨迹的喷管入口起始网格的 左下方顶点；之后，所述数据处理器根据公式计算得出当前 所计算粒子轨迹的喷管入口起始网格内的粒子质量流率式中，M为步 骤301中所设定的需进行处理的不同粒径粒子数量；f_ptog为凝聚相与气相 的流率比且式中n_x为所试验固体推进剂燃烧后所产生的所有 凝聚相产物的平衡组分之和；为当前所处理粒子轨迹的喷管入口起始网 格内的气相流率，且${\stackrel{·}{m}}_g=\frac{{\stackrel{·}{m}}_{g1}+{\stackrel{·}{m}}_{g2}}{2},$式中其中S_网格为当前所处理粒子轨迹的喷管入口起始网格的面积；

步骤4042、粒子数密度计算：所述数据处理器根据公式计算得出当前所计算粒子轨迹在所途经制导信号影响网格内的粒子数密 度；式中，dt为当前所计算粒子轨迹在所途经制导信号影响网格内的驻留 时间，且dt为当前所计算粒子轨迹在所途经制导信号影响网格内的所有轨 迹点的时间步长之和；为步骤4041中计算得出的当前所计算粒子轨迹 的喷管入口起始网格内的粒子质量流率；

步骤4043、粒子平均直径计算：所述数据处理器根据公式 计算得出当前所计算粒子轨迹在所途经制导信号影响网 格内的粒子平均直径D_p；式中，k3为正整数，且k3=1,2，…，K，其中K 为当前所计算粒子轨迹在所途经制导信号影响网格内的轨迹点总数量，且 D_pk3为K个轨迹点中第k3个轨迹点处的粒子直径。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤4043中对当前所计算粒子轨迹在所途经制导信号影响网格内的 粒子平均直径D_p进行计算之前，先根据当前所计算粒子轨迹所途经制导信 号影响网格的上下边界径向坐标和左右边界轴向坐标，并结合步骤4041 中所找出当前所计算粒子轨迹的粒子轨迹数据中各轨迹点的轴向坐标和 径向坐标，对当前所计算粒子轨迹在所途经制导信号影响网格内的轨迹点 总数量K和K个轨迹点的轴向坐标与径向坐标分别进行确定；同时，结合 步骤4041中所找出当前所计算粒子轨迹的粒子轨迹数据中粒子直径随轨 迹变化的文件内所存储内容，对K个轨迹点中各轨迹点的粒子直径进行确 定。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤406中获得入射方向上可见光制导信号的衰减系数κ后，所述数 据处理器还需根据公式κ_t=c×κ对所获得的衰减系数κ进行修正，其中c 为修正系数且c=1.05～1.15。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤一中进行羽流场计算时，所述数据处理器先根据预先设定的喷管 几何参数与喷流计算域范围，建立对所述发动机喷管的内外羽流场进行数 值计算的二维轴对称模型；之后，所述数据处理器调用CFD前处理器，生 成所试验固体推进剂的羽流场计算域网格图，且所述CFD前处理器为 GAMBIT软件；然后，所述数据处理器调用FULENT软件，并结合预先设计 的喷管几何参数、喷流计算域范围和燃烧室工作参数，对所试验固体推进 剂进行羽流场计算，并输出羽流场计算结果；步骤一中对所试验固体推进 剂进行羽流场计算时，所采用凝聚相产物的控制方程为拉格朗日颗粒轨道 模型。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤一中所试验固体推进剂为含铝推进剂；所试验固体推进剂燃烧后 所产生的燃烧产物包括气相产物和凝聚相产物两种类型；步骤一中进行羽 流场计算时，需先输入所试验固体推进剂燃烧后的燃烧产物平衡组成；所 输入的燃烧产物平衡包括M个气相产物的平衡组分和Al₂O₃颗粒的平衡组 分；其中，M为所试验固体推进剂燃烧后所产生气相燃烧产物的数量；所 设定的Al₂O₃颗粒的平衡组分为所试验固体推进剂燃烧后所产生的所有凝 聚相产物的平衡组分之和；步骤404、步骤405和步骤406中所述的M× （N_Y入-1）个粒子轨迹均为Al₂O₃颗粒的轨迹；步骤405中对任一条粒子轨 迹的衰减系数进行换算之前，先通过所述参数输入单元输入Al₂O₃颗粒的 复折射率。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤301中M=8，且8个不同粒径粒子的粒径D_nr分别为D_n1、D_n2、D_n3、 D_n4、D_n5、D_n6、D_n7和D_n8，其中，D_n1＜D_n2＜D_n3＜D_n4＜D_n＜D_n5＜D_n6＜D_n7＜D_n8， 其中，D_n为通过所述参数输入单元预先设定的Al₂O₃颗粒的平均粒径。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤一中进行羽流场计算之前，还需进行能量特性参数计算，其计算 过程如下：

步骤Ⅰ、初始参数设定与存储：

首先，通过所述参数输入单元输入配制所试验固体推进剂所用的组分 数量H和各组分的组分信息，并将所输入的信息同步存储至与所述数据处 理器相接的数据存储单元内；其中，各组分的组分信息均包括化学式和质 量配比m_i，i为正整数，且i=1，2，…，H，H为配制所试验固体推进剂所 用组分的数量；其中，0＜m_i＜100，m₁+m₂+…+m_H=100，N≥2；

之后，通过所述参数输入单元自预先建立的燃烧产物数据库中，选出 所试验固体推进剂燃烧后所产生的所有燃烧产物；所述燃烧产物数据库中 存储有多种燃烧产物的属性信息；其中，各燃烧产物的属性信息均包括化 学式、相对分子质量和相态，其中相态为气相或凝聚相；同时，通过所述 参数输入单元对所试验固体推进剂燃烧后所产生燃烧产物的数量m和气 相产物的数量M进行设定；

步骤Ⅱ、建模：根据最小自由能原理，建立最小自由能数学模型和燃 烧室温度计算模型；

步骤Ⅲ、平衡组成计算：所述数据处理器调用步骤Ⅱ中所建立的最小 自由能数学模型，并结合步骤Ⅰ中所设定的初始参数，计算得出所试验固 体推进剂燃烧后的燃烧产物平衡组成，且计算得出的燃烧产物平衡组成包 括所试验固体推进剂燃烧后所产生的m个燃烧产物的平衡组分；

步骤Ⅳ、绝热燃烧温度计算：所述数据处理器调用步骤Ⅱ中所建立的 燃烧室温度计算模型，计算得出所试验固体推进剂燃烧后处于化学平衡状 态时的绝热燃烧温度。

上述一种固体推进剂烟雾影响可见光制导信号虚拟试验方法，其特征 是：步骤201中对所试验固体推进剂的喷流区域内所有非结构网格节点的 气相产物流场数据进行读取时，首先建立F个动态一维数组，其中 F=F1+F2，其中F1=7，F2为所试验固体推进剂燃烧后所产生气相产物的数 量；F个动态一维数组分别用于存储N个非结构网格点的F个流场变量数 据，并且F个流场变量数据分别为轴向坐标、径向坐标、温度、压强、F2 个气相产物的平衡组分、凝聚相浓度、气相轴向速度和气相密度；其中， N为喷流区域内所有非结构网格节点的数量；步骤201中所述的气相产物 羽流场计算结果为步骤一中羽流场计算完成后存储在*_tec.dat文件内的 计算结果；步骤302中粒子直径随轨迹变化的文件为FULENT软件计算完 成后自动保存的*_diam.fvp文件，粒子轨迹的时间步长文件为FULENT软 件计算完成后自动保存的*_time.fvp文件。

本发明与现有技术相比具有以下优点：

1、方法步骤简单、设计合理且实现方便。

2、投入成本低且使用操作简便，明显简化了固体推进剂的配方研制 过程，大大缩短了配方研制周期，大幅减小了配方研制成本。

3、固体推进剂配方设计过程中，只需输入配制固体推进剂所用各组 分的组分信息，之后采用数据处理器便能自动对固体推进剂进行羽流场模 拟计算与可见光影响计算，因而能简便、快速且准确地对所试验固体推进 剂烟雾对可见光制导性能的影响进行预估。之后，根据烟雾对可见光制导 性能的影响预估结果，对所试验固体推进剂的配方进行相应调整，从而为 固体推进剂的配方设计提供了极大的方便。

4、可见光制导信号影响虚拟测试简单、计算量小、实现方便且计算 结果准确，待完成羽流场计算后，先自动对羽流场计算结果进行结构网格 化处理，且结构网格化处理后，可利用结构网格化处理结果自动计算得出 所试验固体推进剂烟雾对可见光制导性能的影响结果。这样，通过对羽流 场计算结果进行结构网格化处理实现了离散化的衰减系数计算过程，减少 了数据拟合带来的计算误差，因而大幅度提高了可见光制导性能影响结果 的计算精度，并且计算过程简单、实现方便，能快速获得固体推进剂烟雾 对可见光制导信号的影响结果。

5、实用价值高，能对低特征信号固体推进剂的配方研制过程带来极 大便利，实际操作过程中，通过对所试验固体推进剂自动进行羽流场模拟 计算并相应对羽流场计算结果进行结构网格化处理后，能简便、快速且准 确获得固体推进剂烟雾对可见光制导信号的影响结果。实际操作过程中， 只需调整配制所试验固体推进剂所用各组分的重量配比，数据处理器便可 自动完成羽流场模拟计算、结构网格化处理与衰减系数计算过程，因此能 简便、直观且准确得出配制所试验固体推进剂所用各组分的重量配比对所 试验固体推进剂烟雾光学制导特性的影响，因而大大缩短了固体推进剂的 配方研制周期，大幅减小了配方研制成本，为固体推进剂的配方设计提供 了极大的方便。

6、使用效果好且推广应用前景广泛、适用面广，能简便投入推进剂 配方的优化设计过程中，可经济、高效地完成推进剂配方优化设计过程， 并且具有处理结果准确度高、处理速度快等诸多优点，只需几分钟便可完 成结构网格化处理过程，实用性非常强。同时，本发明的推广应用面广， 能有效推广适用至红外、紫外、激光等制导信号。

综上所述，本发明方法步骤简单、设计合理且实现方便、使用效果好， 能简便、快速且准确获得固体推进剂烟雾对可见光制导信号的影响结果。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明的方法流程框图。

图2为本发明羽流场计算时所采用发动机喷管与喷流计算域的结构示 意图。

具体实施方式

如图1所示的一种固体推进剂烟雾影响可见光制导信号虚拟试验方法， 包括以下步骤：

步骤一、羽流场计算及数据存储：采用数据处理器调用FULENT软件 对所试验固体推进剂进行羽流场计算后，将羽流场计算结果自动输出并存 储至与所述数据处理器相接的数据存储单元内；其中，所输出的羽流场计 算结果中包括气相产物羽流场计算结果和凝聚相产物羽流场计算结果。

本实施例中，步骤一中进行羽流场计算时，所述数据处理器先根据预 先设定的喷管几何参数与喷流计算域范围，建立对所述发动机喷管的内外 羽流场进行数值计算的二维轴对称模型；之后，所述数据处理器调用CFD 前处理器，生成所试验固体推进剂的羽流场计算域网格图，且所述CFD前 处理器为GAMBIT软件；然后，所述数据处理器调用FULENT软件，并结合 预先设计的喷管几何参数、喷流计算域范围和燃烧室工作参数，对所试验 固体推进剂进行羽流场计算，并输出羽流场计算结果。

所述发动机喷管的几何参数包括入口半径r₁、喉部半径r₂、出口半径r₃、 收敛段长度d1、喉部圆柱段长度d2、扩张角α和扩张段长度d3，其中入口 半径r₁为发动机喷管的入口半径，喉部半径r₂为发动机喷管的喉部半径， 出口半径r₃为发动机喷管的出口半径，收敛段长度d1为发动机喷管入口到 喷管喉部前端的长度，喉部圆柱段长度d2为发动机喷管的喉部长度，扩张 角α为发动机喷管扩张段的壁面与其轴线之间的夹角，且扩张段长度d3为 发动机喷管喉部末端到喷管出口的长度；喷流计算域范围包括喷流轴向长 度x_m和喷流径向长度y_m，其中喷流轴向长度x_m为发动机喷管出口到喷流下 游出口的长度d4，喷流径向长度y_m为发动机喷管轴线到喷流径向外边界的 长度，详见图2。

本实施例中，步骤一中所试验固体推进剂为含铝推进剂，所述含铝推 进剂为复合推进剂、双基推进剂或改性双基推进剂。

所试验固体推进剂燃烧后所产生的燃烧产物包括气相产物和凝聚相 产物两种类型，所输出的羽流场计算结果中包括气相产物羽流场计算结果 和凝聚相产物羽流场计算结果。

对燃烧室工作参数进行设定时，其设定过程如下：首先，通过所述参 数输入单元对发动机的燃烧室内压强P_c、燃烧室绝热温度、环境压力和环 境温度T_环分别进行设定；之后，再对燃烧产物平衡组成进行设定；然后， 通过所述参数输入单元对所试验固体推进剂燃烧后所产生Al₂O₃颗粒的平 均粒径D_n进行设定。

对燃烧产物平衡组成进行设定时，所设定的燃烧产物平衡组成中包括 所试验固体推进剂燃烧后所产生的M个气相产物的平衡组分和Al₂O₃颗粒 的平衡组分；其中，M为所试验固体推进剂燃烧后所产生气相燃烧产物的 数量；所设定的Al₂O₃颗粒的平衡组分为所试验固体推进剂燃烧后所产生 的所有凝聚相产物的平衡组分之和。

本实施例中，所采用的FLUENT软件为ANSYS FLUENT软件。

本实施例中，步骤一中进行发动机羽流场计算时，所采用气相产物的 的控制方程为湍流模型，且所述湍流模型采用两方程的修正的k-ε模型（即 Realizable k-εturbulent model）。GAMBIT软件生成的所设计固体推进 剂的羽流场计算域网格图中所采用的网格为矩形网格。

并且，步骤一中对所设计固体推进剂进行羽流场计算时，所采用凝聚 相产物的控制方程为拉格朗日颗粒轨道模型。

步骤二、气相产物羽流场数据结构网格化处理，其处理过程如下：

步骤201、喷流区域气相产物流场数据读取：采用所述数据处理器从 步骤一中所输出的气相产物羽流场计算结果中，读取喷流区域内所有非结 构网格节点的气相产物流场数据；所述喷流区域为步骤一中羽流场计算时 所用发动机喷管出口后方的矩形区域。

本实施例中，所设定的入射位置处于所述发动机喷管的中心轴线上。

本实施例中，步骤一中所述气相产物羽流场计算结果为存储在 *_tec.dat文件内的计算结果，且*_tec.dat文件为FULENT软件计算完成 后自动保存的Tecplot格式的文本文件。

例如，当所设计固体推进剂燃烧后所产生的燃烧产物中，包含O₂、H₂O、 CO、CO₂、H、H₂、O、OH、N₂等9个气相组分和凝相组分Al₂O₃时，*_tec.dat 文件的前3行内容如下所示：

TITLE="title"

VARIABLES=X，Y，"temperature"，"pressure"，"molef-OH"， "molef-H₂"，"molef-CO"，"molef-CO₂"，"molef-H₂O"，"molef-O₂"， "molef-O"，"molef-H"，"molef-N₂"，"dpm-concentration"， "axial-velocity"，"radial-velocity"，"mach-number"，"density"， "turb-kinetic-energy"，"turb-diss-rate"，"n2-src"，"oh-src"， "o-src"，"h2-src"，"h-src"，"co2-src"，"co-src"，"h2o-src"，"o2-src"

ZONE T="Rampant"，N=17841，E=17472，ET=QUADRILATERAL，F=FEBLOCK

第1行以TITLE开头，是固定的文件头。

第2行以VARIABLES开头，其后是文件中所保存的29个流场变量名 称，不同变量之间以逗号分隔。这29个流场变量依次为：轴向坐标，径 向坐标，温度，压强，OH组分摩尔分数，H₂组分摩尔分数，CO组分摩尔 分数，CO₂组分摩尔分数，H₂O组分摩尔分数，O₂组分摩尔分数，O组分摩 尔分数，H组分摩尔分数，N₂组分摩尔分数，离散相浓度，气相轴向速度， 气相径向速度，马赫数，气相密度，湍动能，湍流扩散速率，N₂组分源项， OH组分源项，O组分源项，H₂组分源项，H组分源项，CO₂组分源项，CO 组分源项，H₂O组分源项，O₂组分源项。

第3行以VARIABLES开头，其后的N=17841是步骤301中唯一关心的 信息，代表羽流场中共有17841个非结构网格点。

从第4行开始是连续的29个流场变量数据块，首先是N个非结构网 格点的轴向坐标数据（依次对应第1个到第N个非结构网格点），接下来 是N个非结构网格点的径向坐标数据，……，最后是N个O₂组分源项数据。 文本数据之间以单个空格分隔，且为具有小数点后5位有效数字的指数格 式。

本实施例中，步骤201中读取出的读取喷流区域内所有非结构网格节 点的气相产物流场数据，为*_tec.dat文件中所存储的数据。

并且，所读取的流场数据包括非结构网格点数量N以及分别与N个非 结构网格点相对应的多个流场变量数据块。本实施例中，所读取的流场数 据包括非结构网格点数量N=17841和分别与17841个非结构网格点相对应 的29个流场变量数据块。

本实施例中，步骤201中对所试验固体推进剂的喷流区域内所有非结 构网格节点的气相产物流场数据进行读取时，首先建立F个动态一维数组， 其中F=F1+F2，其中F1=7，F2为所试验固体推进剂燃烧后所产生气相产物 的数量；F个动态一维数组分别用于存储N个非结构网格点的F个流场变 量数据，并且F个流场变量数据分别为轴向坐标、径向坐标、温度、压强、 F2个气相产物的平衡组分、凝聚相浓度、气相轴向速度和气相密度；其中， N为喷流区域内所有非结构网格节点的数量。

本实施例中，F2=m-Q=9。且F=F1+F2=16。并且，16个动态一维数组 分别用于存储N个非结构网格点的轴向坐标，径向坐标，温度，压强，OH 组分摩尔分数，H₂组分摩尔分数，CO组分摩尔分数，CO₂组分摩尔分数， H₂O组分摩尔分数，O₂组分摩尔分数，O组分摩尔分数，H组分摩尔分数， N₂组分摩尔分数，凝聚相浓度，气相轴向速度和气相密度；16个动态一维 数组创建后，从*_tec.dat文件中读取对应数据进行存储。并且，16个动 态一维数组中所存储数据的数量均为N个。

步骤202、轴向坐标轴上非结构网格点提取：采用所述数据处理器从 步骤201中所述喷流区域内的所有非结构网格节点中，提取位于轴向坐标 轴上的所有非结构网格点，本步骤中所提取的位于轴向坐标轴上的非结构 网格点总数量为N_X；其中，轴向坐标轴为所述发动机喷管的中心轴线所 在的横坐标轴，位于轴向坐标轴上的非结构网格点的径向坐标y_h=0且其轴 向坐标x_h≥0，其中h为正整数，且h=1，2，…，N_X。

步骤203、径向坐标轴上非结构网格点提取：采用所述数据处理器从 步骤201中所述喷流区域内的所有非结构网格节点中，提取位于径向坐标 轴上的所有非结构网格点，且所提取的位于径向坐标轴上的非结构网格点 总数量为N_Y出；其中，径向坐标轴为所述发动机喷管出口处所在的纵坐标 轴且发动机喷管出口处的轴向坐标值为0，位于径向坐标轴上的非结构网 格点的轴向坐标x_k1=0且其径向坐标y_k1≥0，其中k1为正整数，且k1=1， 2，…，N_Y出。

其中，横坐标轴和纵坐标轴形成一个二维直角坐标系。

步骤204、构建结构网格图：将N_X条直线x=x_h和N_Y出条直线y=y_k1正 交后，构建出一个包含（N_X－1）×（N_Y出－1）个矩形网格的结构网格图。

步骤205、气相产物羽流场数据结构网格化处理：采用所述数据处理 器对步骤204中所构建结构网格图中各矩形网格的四个顶点的气相产物流 场数据分别进行重新赋值；所有矩形网格中各顶点的气相产物流场数据重 新赋值方法均相同，其中对于所构建结构网格图中任一个矩形网格的任一 顶点的气相产物流场数据进行重新赋值时，所述数据处理器先自步骤201 中所述喷流区域内的所有非结构网格节点中找出与当前被赋值顶点距离 最近的非结构网格节点，并将所找出的非结构网格节点的气相产物流场数 据赋给当前被赋值的顶点。

本实施例中，步骤202中提取出位于轴向坐标轴上的所有非结构网格 点后，所述数据处理器将所提取的位于轴向坐标轴上的所有非结构网格点 按照轴向坐标由小至大的顺序进行排列，并由左至右标注于所述轴向坐标 轴上；步骤203中提取出位于径向坐标轴上的所有非结构网格点后，所述 数据处理器将所提取的位于径向坐标轴上的所有非结构网格点按照径向 坐标由小至大的顺序进行排列，并由下至上标注于所述径向坐标轴上。

本实施例中，所述喷流区域为一个矩形区域，且其由GAMBIT软件生 成的羽流场计算域网格图中轴向坐标≥0的多个矩形网格组成。由于转换 到结构网格图后的喷流区域的外边界与GAMBIT软件生成的羽流场计算域 网格图中喷流区域的外边界相互重合。因此，首先确定横坐标轴上（径向 坐标=0，轴向坐标≥0）的所有非结构网格点，上述非结构网格点通过遍 历所述喷流区域中的所有非结构网格节点，并提取出满足径向坐标=0，轴 向坐标≥0的所有非结构网格节点并按照轴向坐标从小到大排序后得到； 之后，确定纵坐标轴（径向坐标≥0，轴向坐标=0）上的所有非结构网格 点，上述非结构网格点通过遍历所述喷流区域中的所有非结构网格节点， 并提取出满足径向坐标≥0，轴向坐标=0的所有非结构网格节点并按照径 向坐标从小到大排序以得到。

待横坐标轴上的所有非结构网格点和纵坐标轴上的所有非结构网格 点均确定后，按照直线正交形成矩形网格的形式便可构建出结构网格图， 所构建结构网格图中矩形网格的数量为（N_X－1）×（N_Y出－1）个，并且 每一个矩形网格的四个顶点均为结构网格点，且所构建结构网格图中结构 网格点的数量为N_X×N_Y出个。

本实施例中，步骤205中进行气相产物羽流场数据结构网格化处理时， 需对所构建结构网格图中N_X×N_Y出个结构网格点分别进行重新赋值。并 且，对各结构网格点进行重新赋值的原则是：从所述喷流区域中的所有非 结构网格节点中筛选出离该结构网格点距离最近的非结构网格节点，并将 所烧选出的非结构网格节点的流场数据赋给该结构网格点。

本实施例中，步骤205中完成气相产物羽流场数据结构网格化处理后， 还需将处理结果同步存储至与所述数据处理器相接的数据存储单元内。

步骤三、凝聚相产物羽流场数据结构网格化处理：所述数据处理器利 用步骤204中所构建的结构网格图对M个不同粒径粒子的粒子轨迹数据分 别进行结构网格化处理，过程如下：

步骤301、结构网格化处理初始参数设定：采用所述参数输入单元对 M的取值和M个不同粒径粒子的粒径D_nr分别进行设定；其中，r为正整数， 且r=1,2，…，M。

本实施例中，步骤301中M=8，且8个不同粒径粒子的粒径D_nr分别为 D_n1、D_n2、D_n3、D_n4、D_n5、D_n6、D_n7和D_n8，其中，D_n1＜D_n2＜D_n3＜D_n4＜D_n＜D_n5＜ D_n6＜D_n7＜D_n8，其中，D_n为通过所述参数输入单元预先设定的Al₂O₃颗粒的 平均粒径。

实际使用时，可以根据实际需要，对M的取值大小进行相应调整，并 且M的取值越大，粒子轨迹网格化数据处理结果越精确。

实际对M个不同粒径粒子的粒径D_nr进行设定时，参照步骤302中所 设定的所设计固体推进剂燃烧后所产生Al₂O₃颗粒的平均粒径D_n进行设 定，具体是在平均粒径D_n的左右两侧分别对称选取多个粒径。也就是说， 在采用FLUNT软件计算得出的粒径分布场中，在平均粒径D_n的左右两侧分 别对称选取多个粒径。

步骤302、粒子轨迹数据读取：采用所述数据处理器从步骤一中所输 出的凝聚相产物羽流场计算结果中读取所设计固体推进剂的所有粒子轨 迹数据；其中，所读取的凝聚相产物羽流场计算结果包括粒子直径随轨迹 变化的文件和粒子轨迹的时间步长文件。

本实施例中，步骤302中所述粒子直径随轨迹变化的文件和粒子轨迹 的时间步长文件，分别为FULENT软件计算完成后自动保存*_diam.fvp和 *_time.fvp文件。

同时，步骤302中所读取的凝聚相产物羽流场计算结果还包括粒子温 度随轨迹变化的文件，即*_temp.fvp，其中*_temp.fvp为文本格式的文件， 且*_diam.fvp、*_temp.fvp和*_time.fvp文件的内容存储格式均相同。 以时间步长文件*_time.fvp为例，其前7行为文件头，内容如下：

以时间步长文件*_time.fvp为例，其前7行为文件头，内容如下：

FVPARTICLES21

Tag Names

0

Variable Names

2

Particle Time Step

particle_id

除了第6行代表变量名称（Particle Time Step）而有所不同外，其 他3个*.fvp（即*_temp.fvp和*_diam.fvp）的文件头完全相同。从第8 行开始记录粒子轨迹数据，内容如下：

1141

-0.04967380.0248277000

-0.04958690.024827307.74289e-0060

-0.04915270.024825503.87144e-0050

-0.04912650.024825402.3338e-0060

-0.04869230.024823503.87165e-0050

第8行是1个整数，代表该条轨迹上的总轨迹点数NL(即1141)。那 么从第9行开始的1141行记录了1141个轨迹点的数据，每行由5个字段 构成。前3个字段分别代表该点的x，y，z坐标，对于二维情况z=0，第 4个字段是该轨迹点处的时间步长，第5个字段为该条轨迹的序号(从0 开始)。第1条轨迹结束后接着是第2条轨迹的数据段，且每条轨迹的数 据段均以只包含一个整数(NL)的行开始，代表从下1行开始的NL行为同 一条轨迹的数据。

步骤303、发动机喷管入口矩形网格数量获取及各喷管入口矩形网格 的上下边界确定：首先，采用所述数据处理器从步骤一中输出的气相产物 羽流场计算结果中读取发动机喷管区域内所有非结构网格节点的流场数 据；之后，采用所述数据处理器从所读取的发动机喷管区域内的所有非结 构网格节点中提取位于直线x=-Δd上的所有非结构网格节点，所提取的位 于直线x=-Δd上的非结构网格点总数量为N_Y入；其中，位于直线x=-Δd上 的非结构网格节点的轴向坐标x_k2=-Δd且其径向坐标y_k2≥0，其中k2为正 整数，且k2=1，2，…，N_Y入；Δd为所述发动机喷管入口至喷管出口之间 的距离；步骤204中所构建结构网格图中位于所述发动机喷管入口处的矩 形网格为喷管入口矩形网格，所述喷管入口矩形网格数量为（N_Y入-1）个， 各喷管入口矩形网格的上下边界分别为上下相邻的两条直线y=y_k2。

步骤四、可见光制导信号影响虚拟测试，其测试过程如下：

步骤401、初始参数设定：采用所述参数输入单元输入试验用可见光 的波长λ。

本实施例中，λ=0.7μm～0.9μm。

步骤402、制导信号入射方向确定：先采用所述参数输入单元对可见 光的入射位置进行设定，再根据所设定的入射位置对可见光的入射方向 x=x_入进行确定，并采用所述数据处理器找出位于直线x=x_入上的所有矩形 网格；其中，x_入为所设定入射位置与所述发动机喷管出口之间的间距；位 于直线x=x_入上的所有矩形网格均为制导信号影响网格。

当x=x_入与相邻两列矩形网格之间的交界线重合时，位于直线x=x_入左 侧的一列矩形网格为位于直线x=x_入上的所有矩形网格。

步骤403、粒子轨迹途经制导信号影响网格确定：采用所述数据处理 器对以步骤303中所述的（N_Y入-1）个喷管入口矩形网格作为起始网格的 所有粒子轨迹途经的制导信号影响网格进行确定；其中，以（N_Y入-1）个 喷管入口矩形网格中任一个喷管入口矩形网格作为起始网格的粒子轨迹 轨迹均包括M个不同粒径粒子轨迹，以（N_Y入-1）个喷管入口矩形网格作 为起始网格的所有粒子轨迹数量为M×（N_Y入-1）个。

由于发动机喷管羽流场中所有的粒子轨迹均起源于发动机喷管入口， 并且各喷管入口矩形网格中的粒子轨迹起始点均为该喷管入口矩形网格 左边界线的中点，其中各喷管入口矩形网格的左边界线均位于纵坐标轴 上。本实施例中，以各喷管入口矩形网格作为起始网格的粒子轨迹的数量 均为M个。

步骤403中M个不同粒径粒子轨迹分别为粒径D_n1的粒子轨迹、粒径 为D_n2的粒子轨迹、…、和粒径为D_nM的粒子轨迹。

本实施例中，步骤403中对以（N_Y入-1）个喷管入口矩形网格作为起 始网格的所有粒子轨迹途经的制导信号影响网格进行确定时，所有粒子轨 迹途经的制导信号影响网格确定方法均相同；对于任一个粒子轨迹途经的 制导信号影响网格进行确定时，根据步骤302中所读取的当前所处理粒子 轨迹的粒子轨迹数据中各轨迹点的轴向坐标和径向坐标，并结合步骤204 中所构建的结构网格图，找出当前所处理粒子轨迹的途经的制导信号影响 网格。

本实施例中，步骤403中进行初始参数设定时，还需采用所述参数输 入单元输入试验用可见光的入射光程长L，其中L=2×y_m，y_m为喷流径向 长度且其为发动机喷管轴线到喷流径向外边界的长度；步骤406中获得入 射方向上可见光制导信号的衰减系数κ后，所述数据处理器还需根据公式 T=e^-κ×L计算得出试验用可见光的光谱透射比T。

步骤404、粒子轨迹网格化数据计算：采用所述数据处理器对M× （N_Y入-1）个粒子轨迹的网格化数据分别进行计算，所有粒子轨迹的网格 化数据计算方法均相同；其中，对M×（N_Y入-1）个粒子轨迹中任一条粒 子轨迹的网格化数据进行计算时，采用所述数据处理器计算得出当前所计 算粒子轨迹在所途经制导信号影响网格内的粒子数密度和粒子平均直径。

本实施例中，步骤404中对M×（N_Y入-1）个粒子轨迹中任一条粒子 轨迹的网格化数据进行计算时，过程如下：

步骤4041、发动机喷管入口起始网格确定及粒子质量流率计算： 采用所述数据处理器从步骤302中所读取的所有粒子轨迹数据中找出当前 所计算粒子轨迹的粒子轨迹数据，并根据所找出的粒子轨迹数据中的发动 机喷管入口处轨迹点的径向坐标和步骤303中所确定的各喷管入口矩形网 格的上下边界，对当前所计算粒子轨迹的喷管入口起始网格进行确定；

待当前所计算粒子轨迹的喷管入口起始网格确定后，所述数据处理器 找出所确定喷管入口起始网格的顶点一和顶点二的气相产物流场数据，并 从顶点一的气相产物流场数据中找出顶点一处的气相密度ρ_g1和气相轴向 速度u_g1，且从顶点二的气相产物流场数据中找出顶点二处的气相密度ρ_g2和气相轴向速度u_g2，其中顶点一为当前所计算粒子轨迹的喷管入口起始网 格的左上方顶点，且顶点二为当前所计算粒子轨迹的喷管入口起始网格的 左下方顶点；之后，所述数据处理器根据公式计算得出当前 所计算粒子轨迹的喷管入口起始网格内的粒子质量流率式中，M为步 骤301中所设定的需进行处理的不同粒径粒子数量；f_ptog为凝聚相与气相 的流率比且式中n_x为所试验固体推进剂燃烧后所产生的所有 凝聚相产物的平衡组分之和；为当前所处理粒子轨迹的喷管入口起始网 格内的气相流率，且${\stackrel{·}{m}}_g=\frac{{\stackrel{·}{m}}_{g1}+{\stackrel{·}{m}}_{g2}}{2},$式中其中S_网格为当前所处理粒子轨迹的喷管入口起始网格的面积；

步骤4042、粒子数密度计算：所述数据处理器根据公式计算得出当前所计算粒子轨迹在所途经制导信号影响网格内的粒子数密 度；式中，dt为当前所计算粒子轨迹在所途经制导信号影响网格内的驻留 时间，且dt为当前所计算粒子轨迹在所途经制导信号影响网格内的所有轨 迹点的时间步长之和；为步骤4041中计算得出的当前所计算粒子轨迹 的喷管入口起始网格内的粒子质量流率；

步骤4043、粒子平均直径计算：所述数据处理器根据公式 计算得出当前所计算粒子轨迹在所途经制导信号影响网 格内的粒子平均直径D_p；式中，k3为正整数，且k3=1,2，…，K，其中K 为当前所计算粒子轨迹在所途经制导信号影响网格内的轨迹点总数量，且 D_pk3为K个轨迹点中第k3个轨迹点处的粒子直径。

本实施例中，步骤4043中对当前所计算粒子轨迹在所途经制导信号 影响网格内的粒子平均直径D_p进行计算之前，先根据当前所计算粒子轨迹 所途经制导信号影响网格的上下边界径向坐标和左右边界轴向坐标，并结 合步骤4041中所找出当前所计算粒子轨迹的粒子轨迹数据中各轨迹点的 轴向坐标和径向坐标，对当前所计算粒子轨迹在所途经制导信号影响网格 内的轨迹点总数量K和K个轨迹点的轴向坐标与径向坐标分别进行确定； 同时，结合步骤4041中所找出当前所计算粒子轨迹的粒子轨迹数据中粒 子直径随轨迹变化的文件内所存储内容，对K个轨迹点中各轨迹点的粒子 直径进行确定。

实际对当前所计算粒子轨迹在所途经制导信号影响网格内的轨迹点 总数量K进行确定时，遍历当前所计算粒子轨迹上的轨迹点，并找出轴向 坐标在所途经制导信号影响网格的左右边界轴向坐标之间且径向坐标在 当前所途经制导信号影响网格的上下边界径向坐标之间的轨迹点。K个轨 迹点将所途经制导信号影响网格内的当前所计算粒子轨迹分成多个轨迹 段，且多个所述轨迹段的时间步长之和（即K个轨迹点的时间步长之和） 便为dt。

步骤405、各粒子轨迹衰减系数换算：采用所述数据处理器对M× （N_Y入-1）个粒子轨迹的衰减系数分别进行换算，所有粒子轨迹的衰减系 数换算方法均相同；其中，对M×（N_Y入-1）个粒子轨迹中任一条粒子轨 迹的衰减系数进行换算时，所述数据处理器根据公式κ₀=πr²×Q×n_p进行换 算，式中κ₀为当前所换算粒子轨迹的衰减系数，d为步骤404中计 算得出的当前所换算粒子轨迹在所途经制导信号影响网格内的粒子平均 直径，n_p为步骤404中计算得出的当前所换算粒子轨迹在所途经制导信号 影响网格内的粒子数密度；Q为衰减效率且式中 a_n和b_n为入射波长为λ时半径为r的球形粒子的Mie散射系数， Re(a_n+b_n)表示取a_n+b_n的实部，n为正整数。

本实施例中，对入射波长为λ时半径为r的球形粒子的Mie散射系数a_n和b_n进行确定时，参照以下四个文献进行确定。其中，文献一为 W.J．Lentz.Generating Bessel functions in mie scattering  calculations using continued fractions[J].Applied Optics，1976， 4（3）。文献二为W.J．Wiscombe.Improved mie scattering  algorithms[J].Applied Optics，1980，10（6）。文献三为 Y.Xu.Electromagnetic scattering by and aggregate of  spheres[J].Applied Optics，1985，8（3）。文献四为W.Yang．Improved  recursive algorithm for light scattering by a multilayered  sphere[J].Applied Optics，2003，11（5）。

同时，也可参照2007年《固体火箭技术》第30期第5卷出版的《固 体火箭发动机排气羽流对红外制导信号衰减的计算方法》一文，作者：张 硕、王宁飞、张平。

步骤406、入射方向上可见光制导信号衰减系数计算：采用所述数据 处理器将步骤405中计算得出的M×（N_Y入-1）个粒子轨迹的衰减系数叠 加后，获得入射方向上可见光制导信号的衰减系数κ。

本实施例中，步骤406中获得入射方向上可见光制导信号的衰减系数 κ后，所述数据处理器还需根据公式κ_t=c×κ对所获得的衰减系数κ进行修 正，其中c为修正系数且c=1.05～1.15。

实际使用时，可以根据具体需要，对c的取值大小进行相应调整。本 实施例中，c=1.1。

本实施例中，步骤一中所试验固体推进剂为含铝推进剂；所试验固体 推进剂燃烧后所产生的燃烧产物包括气相产物和凝聚相产物两种类型；步 骤一中进行羽流场计算时，需先输入所试验固体推进剂燃烧后的燃烧产物 平衡组成；所输入的燃烧产物平衡包括M个气相产物的平衡组分和Al₂O₃颗粒的平衡组分；其中，M为所试验固体推进剂燃烧后所产生气相燃烧产 物的数量；所设定的Al₂O₃颗粒的平衡组分为所试验固体推进剂燃烧后所 产生的所有凝聚相产物的平衡组分之和；步骤404、步骤405和步骤406 中所述的M×（N_Y入-1）个粒子轨迹均为Al₂O₃颗粒的轨迹；步骤405中对 任一条粒子轨迹的衰减系数进行换算之前，先通过所述参数输入单元输入 Al₂O₃颗粒的复折射率。

实际使用过程中，对Al₂O₃颗粒的复折射率进行确定时，参照本领域 技术人员能获得的有关Al₂O₃颗粒复折射率的参考文献进行设定。

本实施例中，参照2008年哈尔滨工业大学博士论文《三氧化二铝粒子 复折射率的实验研究》，作者：齐宏，对Al₂O₃颗粒的复折射率进行设定。

本实施例中，步骤405中对任一条粒子轨迹的衰减系数进行换算之前， 所述数据处理器还需根据公式计算得出当前所计算粒子轨 迹在所途经制导信号影响网格内的平均温度T_p，所输入Al₂O₃颗粒的复折 射率为Al₂O₃颗粒在温度T_p时的复折射率，T_pk3为步骤4043中K个轨迹点中 第k3个轨迹点处的粒子温度（从*_temp.fvp文件中读取）。

本实施例中，步骤一中进行羽流场计算之前，还需进行能量特性参数 计算，其计算过程如下：

步骤Ⅰ、初始参数设定与存储：

首先，通过所述参数输入单元输入配制所试验固体推进剂所用的组分 数量H和各组分的组分信息，并将所输入的信息同步存储至与所述数据处 理器相接的数据存储单元内；其中，各组分的组分信息均包括化学式和质 量配比m_i，i为正整数，且i=1，2，…，H，H为配制所试验固体推进剂所 用组分的数量；其中，0＜m_i＜100，m₁+m₂+…+m_H=100，N≥2；

之后，通过所述参数输入单元自预先建立的燃烧产物数据库中，选出 所试验固体推进剂燃烧后所产生的所有燃烧产物；所述燃烧产物数据库中 存储有多种燃烧产物的属性信息；其中，各燃烧产物的属性信息均包括化 学式、相对分子质量和相态，其中相态为气相或凝聚相；同时，通过所述 参数输入单元对所试验固体推进剂燃烧后所产生燃烧产物的数量m和气 相产物的数量M进行设定；

步骤Ⅱ、建模：根据最小自由能原理，建立最小自由能数学模型和燃 烧室温度计算模型；

步骤Ⅲ、平衡组成计算：所述数据处理器调用步骤Ⅱ中所建立的最小 自由能数学模型，并结合步骤Ⅰ中所设定的初始参数，计算得出所试验固 体推进剂燃烧后的燃烧产物平衡组成，且计算得出的燃烧产物平衡组成包 括所试验固体推进剂燃烧后所产生的m个燃烧产物的平衡组分；

步骤Ⅳ、绝热燃烧温度计算：所述数据处理器调用步骤Ⅱ中所建立的 燃烧室温度计算模型，计算得出所试验固体推进剂燃烧后处于化学平衡状 态时的绝热燃烧温度。

步骤一中进行羽流场计算时，结合能量特性参数计算结果，对燃烧室 工作参数进行设定。

本实施例中，对燃烧产物平衡组成进行设定时，所设定的燃烧产物平 衡组成为根据吉布斯最小自由能原理，计算得出燃烧产物的化学平衡组 成。

根据热力学原理，在高温条件下固体推进剂的燃烧产物可视为理想气 体，则整个系统的自由能就等于该系统各组分自由能的总和，己知物质的 自由能是压力、温度和浓度的函数，当该体系达到化学平衡时，体系的自 由能最小。因此，在一定的压力和温度条件下，求出既能使物系自由能最 小又符合质量守恒定律的一组分值，则这组组分值即为该条件下物系的燃 烧产物平衡组成。根据系统达到化学平衡时其自由能函数之总和为最小的 原理，采用迅速收敛的数学方法，通过重复迭代可解任何复杂系统的化学 平衡组成。

本实施例中，步骤Ⅱ中所建立的最小自由能数学模型为 (1)，式(1)中：j为正整数，且j=1，2，…， B，B为固体推进剂所含化学元素的种类；s为正整数，且s=1，2，…，m， m为固体推进剂燃烧后处于化学平衡状态时所含燃烧产物的种数； 其中μ_s为预先通过所述参数输入单元输入的第s种燃烧产物 的化学势（KJ/mol），n_s为1000g固体推进剂燃烧后处于化学平衡状态时 所含第s种燃烧产物的摩尔数（mol/Kg）且n_s≥0，a_sj为1mol第s种燃烧 产物中所含第j种化学元素的原子数；b_j为1000g固体推进剂中所含第j种 化学元素的原子数，π_j为拉格朗日乘子。

实际使用时，最小自由能法就是求出在满足式（1）条件下的一组ns值 使得体系自由能最小，其中s=1、2…m，且n_s≥0，这是一个多元函数的 条件极值问题，用拉格朗日乘数法则可求解。

实际对化学平衡组成进行求解时，也可以采用河南科学技术出版社于 1999年出版的由田德余、刘剑洪著的《化学推进剂计算能量学》一书中， 第六章“能量特性计算的基本方程”中第6.3节“含凝聚相燃烧产物的 化学平衡组成”中所记载的化学平衡组成计算方法进行计算。

所建立的燃烧室温度计算模型为绝热温度计算模型，且所建立的绝热 温度计算模型为(2)，式(2)中H_c1为绝热温度 T=T₁时的1000g燃烧产物的焓，H_c2为绝热温度T=T₂时的1000g燃烧产物 的焓，H_c1＜H_p＜H_c2，并且T_c1和T_c2均根据公式其中n_s为 1000g固体推进剂燃烧后处于化学平衡状态时所含第s种燃烧产物的摩尔 数（mol/Kg），H_cs为1mol第s种燃烧产物在绝热温度为T时的焓， H_cs=RT(α_s1+α_s2T2+α_s3T²3)+α_s4T³4+α_s5T⁴5+α_s6T⁵6），其中，R为通用气 体常数（Kg·m/mol·K），T为绝热温度，α_s1、α_s2、α_s3、α_s4、α_s5和α_s6为 预先通过所述参数输入单元输入的第s种燃烧产物的热力学函数温度系 数；其中M_i为配制固体推进剂所用第i种组分的相对分 子质量，H_i为预先通过所述参数输入单元输入的1mol第i种组分的焓，W_i为配制固体推进剂所用第i种组分的质量百分数。

本实施例中，步骤Ⅲ中进行平衡组成计算时，所述数据处理器调用参 数计算模块，根据步骤Ⅰ中所输入的配制所设计固体推进剂所用各组分的 化学式和质量配比m_i，并结合各组分的相对分子质量，对a_sj和b_j进行计算； 之后，所述数据处理器结合预先输入的μ_s和π_j，且根据公式 (1)计算得出n_s，此时便获得所设计固体推进 剂燃烧后的燃烧产物平衡组成；

步骤Ⅳ中进行绝热燃烧温度计算时，所述数据处理器结合步骤Ⅲ中计 算得出的所设计固体推进剂燃烧后的燃烧产物平衡组成，并根据公式 (2)，计算得出燃烧室温度T_c；其中， 其中M_i为配制所设计固体推进剂所用第i种组分的相对 分子质量，H_i为1mol第i种组分的焓，W_i为配制所设计固体推进剂所用第 i种组分的质量百分数。其中，W_i与步骤Ⅰ中输入配制所设计固体推进剂 所用各组分的质量配比m_i一致。

另外，步骤一中进行羽流场计算之前，还需先通过所述参数输入单元 输入燃烧产物的有限速率化学反应模型。所输入的有限速率化学反应模 型，与步骤一中通过所述参数输入单元自预先建立的燃烧产物数据库中， 选出的所评价固体推进剂燃烧后所产生的燃烧产物相对应。

实际进行羽流场计算时，输入燃烧产物的有限速率化学反应模型时， 参照1991年出版的王荪源主编的《火箭装药热力性能手册》、2002年科 学出版社出版的赵坚行主编的《燃烧的数值模拟》或2008年09月国防科 技大学出版社出版的杨涛、方丁西、唐乾刚主编的《火箭发动机燃烧原理》 等文献对所输入的有限速率化学反应模型进行确定。

例如，当固体推进剂燃烧后所产生的燃烧产物中，包含O₂、H₂O、CO、 CO₂、H、H₂、O、OH、N₂等9个气相组分和凝相组分Al₂O₃时，化学反应模 型为9组分10反应的化学动力学模型，且所采用的9组分10反应的化学 动力学模型详见表2：

表29组分10反应的化学动力学模型反应机理数据表

化学反应式 A(cm3/mol·s) n E(j) CO+O+M=CO₂+M 8.310E-12 0 -9.70E+03 CO+OH=CO₂+H 6.323E+06 -1.5 -2.08E+03 H₂+OH=H₂O+H 1.024E+08 1.6 1.38E+04 H₂+O=OH+H 5.119E+04 2.67 2.63E+04 H+O₂=OH+O 1.987E+14 0 7.03E+04 OH+OH=H₂O+O 1.506E+09 1.14 4.14E+02 H+H+M=H₂+M 1.493E-06 -1 0 O+O+M=O₂+M 2.409E-07 -1 0 O+H+M=OH+M 7.829E-06 -1 0 H+OH+M=H₂O+M 3.673E-02 -2 0

当固体推进剂燃烧后所产生的燃烧产物中，包含O₂、H₂O、CO、CO₂、H、 H₂、O、OH、N₂、HCL、CL、CL₂等12个气相组分和凝相组分Al₂O₃时，化学 反应模型为12组分17的化学动力学模型，且所采用的12组分17反应的 化学动力学模型详见表3：

表312组分17反应的化学动力学模型反应机理数据表

以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是 根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构 变化，均仍属于本发明技术方案的保护范围内。