具体实施方式

为了使本领域技术人员更好的理解本发明的技术方案并能予以实施，下面结合附图和具体实施例对本发明进行详细说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

实施例1

本发明提供了一种固体火箭发动机燃面退移计算方法，包括以下步骤：

步骤1、将固体火箭发动机推进剂装药几何模型离散为三角面片模型，将所有三角面片分类为两个三角面片集合。具体的，本实施例中，三角面片分类方法具体为：点火时刻即开始燃烧的那些三角面片的集合记为B，B中所有三角面片组成一个薄片型的曲面。模型中所有的三角面片的集合记为G(这里的所有三角面片包括所有三角面片，B是G的子集)，G中所有三角面片组成一个封闭的几何体。如图1所示，粗线与虚线共同组成固体火箭发动机推进剂装药几何模型(集合G)，但“初始时燃烧的表面”组仅包含虚线(集合B)。

步骤2、如图2所示，对输入的固体火箭发动机所在的三维包络盒划分笛卡尔网格，该网格的每一个节点的节点坐标记为x。

步骤3、对于三角面片集合B，在笛卡尔网格的每一个节点上，使用AABB(Axis-aligned bounding boxes)快速距离场算法，计算从节点坐标x出发到三角面片集B中任一三角面片的最短距离，形成定义在笛卡尔网格所有节点上的最短距离场函数

步骤4、基于固体推进剂燃速模型，计算每一个节点坐标x处的推进剂燃速r，无需考虑节点属于燃烧表面还是非燃烧表面，形成笛卡尔网格上的燃速分布场r(x)。

其中，r(x)＝ap^n，a和n分别为节点坐标x处的推进剂的燃速系数和压强指数，由推进剂物性参数得到。p为当前时刻固体火箭发动机燃烧室内的压强。

步骤5、在每一个节点上，执行计算修改场函数在每个节点坐标x处的数值。

步骤6、使用Marching cubes等值面抽取算法，基于新的场函数生成用三角面片表达的等值面，将该等值面中的所有三角面片的集合记为I，I由大量三角面片构成，如图3所示。

步骤7、进行几何布尔运算，将等值面I处于封闭几何体G之外的部分剪切掉，更新I集合使其仅包含剪切后的等值面中的三角面片。如图4为已经剪切去除了推进剂药柱之外部分的等值面，同时也是当前时刻的燃烧的表面，此时，剪切过的I的面积即为当前时刻的推进剂燃烧表面面积，具体包括以下步骤：

遍历I中的每一个三角形t，判定t是否处在G内部。

如果t处于G内部则不做任何操作。

如果t处于G所表示的封闭几何体外部，则从I中删除t。

如果t与G中的若干三角形相交(相交判定过程中可以使用AABB(Axis-alignedbounding boxes)快速相交测试算法进行加速)，则将t沿交线分割为“处在G内部的部分”和“处在G外部的部分”，将“处在G内部的部分”分割为若干三角形(取决于具体形状)，并把这些三角形加入I，从I中删除t。

步骤8、组成I的每一个三角面片t的面积乘以该处的燃烧速度，再乘以推进剂密度，即为当前时刻t所在位置的燃气生成率，将所有三角面片对应的燃气生成率数据加和，即可得到该固体火箭发动机在当前时刻的总燃气生成率，记录燃烧表面面积和所有燃气生成率数据。

步骤9、根据燃烧表面面积和燃气生成率数据，基于内弹道曲线数值模型(例如准稳态或瞬态零维内弹道模型)计算固体火箭发动机燃烧室内当前时刻的压强p。

步骤10、将剪切过的I视作初始时燃烧的表面，用I集合取代B集合，在笛卡尔网格的每一个节点上，使用AABB快速距离场算法，计算节点到三角面片集B的最短距离，形成定义在笛卡尔网格所有节点上的一个新的最短距离场函数遍历所有网格节点，检查对应的是否小于零。如果是，则将的数值修改为最后，用修改后的取代原有的

步骤11、基于步骤10中得到的新的和步骤9中得到的新的p，循环执行上述步骤4-10。直至I中不包含任何三角面片。此时，推进剂药柱完全燃尽。

燃烧表面面积变化曲线、燃气生成率曲线，以及压强p随时间变化的曲线，即为燃面退移计算结果。

该方法整个计算流程不涉及数值差分运算，节省了算力，避免了数值差分运算在在外角平分线区域带来的燃烧表面形状畸变和燃速不连续点附近的数值耗散现象和抹平效应。

以上所述实施例仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换，均属于本发明的保护范围。