翼柱形药柱燃面退移过程的变量化设计方法①

蔡 强，鲍福廷，丁 林，刘 旸

(西北工业大学 航天学院，西安 710072)

0 引言

固体火箭发动机翼柱形药柱具有三维复杂构型，其燃面面积很难用解析公式表达。实体造型方法［1］是目前采用的主要方法，通常对商用CAD软件进行二次开发，避免了复杂的数理方程推导，却容易因特征定义不严格，而出现特征消失或自相交，导致几何模型拓扑畸变而造型失败，且商用CAD软件计算大型药柱的时间较长，不利于方案的验证与药柱的优化设计。文献［2］中的方法仅能处理具有1组翼的药柱，其采用基于造型历史追溯的燃面算法，容易发生统计疏漏和重复。

为了解决上述问题，本文结合翼柱形药柱构型的特点，引入变量化设计方法，在几何造型内核上自主开发。通过定义统一的翼形体特征函数，给出翼形体的关键设计变量，严格控制模型参数的变化，避免了造型失败;调节三维图形变换矩阵中的比例变换因子，缩小几何模型的尺寸，显著缩短了大型翼柱形药柱的仿真时间;本文还研究新的燃面面积算法，避免统计燃面出错，并针对后翼、前后翼和大小翼等各种类型的翼柱形药柱，总结了统一的设计流程。

1 翼柱形药柱变量化仿真建模

1.1 假设条件

(1)整个翼柱形药柱表面光滑均一，没有裂纹、尖点等缺陷。

(2)翼柱形药柱的燃面同时被点燃，燃面各点的条件相同。

(3)燃面退移满足平行层燃烧定律，即燃面各点都以相同燃速向药柱内推进。

(4)翼柱形药柱的各组翼形体周向均匀分布，满足严格的循环对称条件。

1.2 建模仿真过程

本文将翼柱形药柱分为外轮廓、回转体和翼形体3种特征形体，其变量化设计过程可分为变量化约束草图绘制、特征形体定义与生成和燃面退移仿真3步，如图1所示。首先，交互式绘制外轮廓、回转体和翼形体的二维截面草图，并对草图施加约束，定义设计变量与变量方程。然后，根据药柱相应特征的二维草图，采用几何造型算法生成外轮廓、回转体和翼形体的三维几何实体，外轮廓依次“减去”回转体和翼形体，就得到药柱的初始形状，并调用统计函数获得初始燃面面积及质量特性信息。最后，仿真燃面退移时，按照平行层燃烧定律假设，每次退移一定的肉厚，通过改变回转体和翼形体的特征变量，使得变化后的各形体表面与初始表面等距，外轮廓“减去”变化后的回转体和翼形体，得到一系列退移后的药柱形状，计算燃烧面积和质量特性，得到燃面面积、质心、惯量、形状的变化规律。

图1 翼柱形药柱变量化设计过程Fig.1 Variational design principle for finocyl grain

2 翼柱形药柱变量化设计关键技术

在翼柱形药柱的实际应用过程中，采用在药柱前段增加1组轴向前翼，或在后段周向增加另1组小翼的方法来增大初始燃面，从而为助推段提供更大推力，前者称为前后翼，后者称为大小翼。从药柱二维截面生成三维特征形体以及燃面退移时，为了简便，通常只处理循环对称药柱的1个对称单元，然而前后两组翼通常不在同一对称平面内，大小2组翼是周向阵列而成，因此需采取措施，确保所有翼形体都处在一份对称单元中，从而得到合理的布尔运算结果。另外，还要解决翼形体的通用描述方法、稳健准确的燃面面积算法和缩短仿真时间等关键技术。

2.1 图形比例变换

本文采用的是支持边界表达方法的精确几何造型内核［3］，药柱的尺寸规模越大，耗费机时越长。因此，一个直观的想法就是在保证造型精度的条件下，对外轮廓、回转体和翼形体的相应二维截面进行图形比例变换，以缩小尺寸［4］。这样生成的三维特征形体尺寸较小，能显著缩短后续的布尔运算和等距偏移运算所耗费的时间。

式(1)为三维图形比例变换矩阵:

式中 (x，y，z，1)为原始截面各点的齐次坐标;(x*，y*，z*，1)为缩放后的特征截面相应点的齐次坐标;(xf，yf，zf)为比例变换的参考点坐标;Sx、Sy和 Sz分别为x、y和z坐标分量的比例变换因子，比例变换因子的取值以不损失造型精度的前提下越小越好。

翼柱形药柱特征形体原始截面的各点依次乘以比例变换矩阵，就得到特征形体缩放后截面的各点坐标。

2.2 翼形体特征函数定义

考虑到通用性，本文定义翼形体特征函数如下:

式中 Section为翼形体二维特征截面;w为翼形体厚度;r为翼形体顶端倒角半径(见图2);angle为该组翼形体起始角度;symN为药柱的循环对称数目;npc为每个外轮廓对称单元中该组翼的个数。

图2 翼形体示意图Fig.2 Parameters of fin entity

药柱的循环对称数目symN是前后2组翼个数(nf，na)或大小2组翼个数(nb，ns)的最大公约数，则

翼在外轮廓对称单元中的周向排列如图3所示。翼形体起始角度angle表示第一个翼距外轮廓对称单元中心线的角度，其余翼的角度为

图3 外轮廓对称单元中多个翼的角度分布Fig.3 Angles of multiple sets of fins defined in one cyclic symmetry cell

每组翼的第1个翼作为几何造型的样本，其余通过将样本的副本绕X轴进行三维图形旋转变换得到［4］，即

式中 (x，y，z，1)为样本点的齐次坐标;(x*，y*，z*，1)为副本点的齐次坐标。

经过如上定义以后，即可确保所有翼形体均处在外轮廓实体对称单元中，从而得到严格旋转对称的1/symN份药柱实体，然后周向阵列即可得到完整的药柱信息。

2.3 颜色标识的燃面面积计算方法

本文采用边界表达的几何模型，其数据结构完整记录了模型的拓扑、几何和属性信息［5］，颜色作为基本属性，附加在药柱外轮廓实体和内腔特征实体的几何模型上，并随着实体模型之间的布尔运算及图形变换，保持其合理性，具体体现在:

式中 药柱实体Sgrain的表面来源于2部分:由外轮廓实体Sout的几何模型表面分裂而成，记为SFa;由内腔特征实体Score的几何表面分裂而成，记为SFb。

本文用不同的颜色标记燃面和非燃面:

(1)外轮廓实体所有表面标识为蓝色，内腔特征实体表面标识为红色。因此，所有的SFa为蓝色，所有的SFb为红色。

(2)按照式(6)进行布尔差运算，得到的药柱实体的表面中，源自外轮廓实体的为蓝色，源自内腔特征实体的为红色。显然，红色的面即为翼柱形药柱的燃面，统计所有红色面的面积，即为药柱的初始燃面面积。

(3)内腔特征实体扩张以后，仍保持所有表面为红色。同理，与不发生变化的外轮廓实体进行布尔差运算之后，统计得到的扩张后药柱实体的红色表面面积，即为退移后的燃面面积。

这样，只需在初始时进行1次颜色标识，就可保证燃面统计不发生疏漏，得到精确的计算结果。

3 翼柱形药柱的变量化设计流程

在解决以上关键技术后，整个翼柱形药柱的变量化设计流程总结见图4。

图4 翼柱形药柱的变量化设计流程Fig.4 Variational design flow chart of finocyl grain

(1)在变量化设计环境中，交互式绘制外轮廓、回转形内孔、前后翼或大小翼的约束草图，得到各个特征实体的原始二维截面。

(2)按照2.1节的方法对特征实体的截面执行比例变换，缩小图形的尺寸规模，以节省求解时间。

(3)由特征截面分别生成外轮廓的对称单元、回转形内孔的对称单元、前后翼或大小翼的三维实体，并按照2.2节中的处理方法，确保所有翼形体在外轮廓对称单元中。

(4)按照2.3节中介绍的颜色标识的燃面面积计算方法，计算初始燃面面积(肉厚Web=0)和退移一定距离(肉厚Web＞0)以后的燃面面积，每退移一步依次计算燃面面积及质量特性。

(5)判断药柱是否完全燃烧，即判断当前药柱质量是否为0，如果完全燃烧，就输出燃面面积及质量特性随烧去肉厚的变化关系;否则，继续执行退移仿真。

4 实例验证

本文以LEDAS约束求解器［6］，构建变量化设计环境，实现了上述针对翼柱形药柱的处理方法。

4.1 可行性实例验证及结果分析

为了验证该方法的可行性，以1个有3组翼形体的药柱作为示例，截面尺寸如图5所示。3组翼分为前翼、后大翼、后小翼，其参数见表1。

图5 大小翼和前后翼药柱二维截面草图Fig.5 2-D geometric constraint sketch of rear big fin and rear small fin

表1 3组翼的参数Table 1 Parameters of three sets of fins

药柱燃面退移过程及燃面面积-肉厚关系曲线如图6和图7所示。

图6 大小翼和前后翼药柱燃面退移过程仿真及燃面面积-肉厚曲线Fig.6 Grain burnback simulation and profile of burning surface area vs web

这个燃面退移过程分为3段，初始阶段燃面不断增大，直至肉厚Web=30 mm时，后大翼开始烧穿;Web=30～60 mm为减面燃烧，前翼和后小翼逐渐扩张，后大翼逐渐消失;Web=60 mm时，后大翼完全消失，前翼和后小翼开始烧穿，直至Web=85 mm药柱完全消失，这一阶段也为减面燃烧。仿真时，比例变换因子Sx=Sy=Sz=1，肉厚推进步长为1 mm，85步计算共耗时221.63 s，如果令 Sx=Sy=Sz=0.02，即缩小 50 倍，计算时间缩短至151.16 s(CPU主频2.8 GHz)。

图7 参考文献实例的燃面退移仿真过程及燃面面积-肉厚曲线Fig.7 Grain burnback simulation and profile of burning surface area vs web

4.2 计算精度实例验证及结果分析

为了验证本文计算方法的精度，对文献［7］中17#药柱进行燃面退移仿真，并结合推进剂参数和喷管参数，进行了内弹道计算。整个药柱肉厚为21.35 mm，燃烧初始阶段(Web=0～4 mm)内孔和翼形体同时扩张，为增面燃烧;在Web=4 mm时，后端翼形体开始烧穿，直至Web=12 mm翼形体完全烧完，这一阶段为减面燃烧;燃烧最后阶段(Web=12～21.35 mm)为近似等面燃烧过程。该例本文设置的比例变换因子为Sx=Sy=Sz=2/3，仿真肉厚步长与文献［7］相同，都为1 mm，22步仿真时间仅为7.344 s，而文献［7］中仿真耗费了将近30 min。

本文内弹道计算结果与文献［7］中压力-时间进行对比，如图8所示。由图8可看出，由于本文内弹道计算未考虑点火段和拖尾段的非定常效应，点火段和拖尾段与试验曲线差距较大，其余吻合相当好(最大误差仅3.1%)。

图8 压力-时间曲线计算结果与试验曲线对比Fig.8 Results of pressure vs time compared with experiment data

5 结论

(1)引入翼的起始角度，合理布置大小翼、前后翼的位置，使其在1个外轮廓实体对称单元中。调节三维比例变换矩阵的比例变换因子，以缩小药柱的尺寸规模，显著缩短求解时间。

(2)在基于造型历史追溯的燃面面积统计方法的基础上，提出改进的颜色标识燃面面积计算方法，该方法具有执行效率高、操作简单、不会统计疏漏等优点。

(3)归纳总结了翼柱形药柱的变量化设计流程，包括特征形体变量化约束草图的定义、比例变换、燃面退移过程仿真、统计燃面面积等步骤。并通过2个实例，对其有效性及计算精度进行了验证。

(4)文中研究工作为翼柱形药柱的设计、参数优化和燃面退移仿真提供一套实用的解决方法。同时，本文的研究思路也适用于其他复杂的任意三维药型。

［1］董新刚，陈林泉，侯晓.基于Pro/E平台下的固发装药CAD软件［C］//昆明:2002年中国宇航学会固体推进专业委员会年会论文集(上)，2002:109-114.

［2］刘旸，鲍福廷，蔡强，等.变量化约束草图驱动的药柱三维特征设计研究［J］.固体火箭技术，2010，33(1).

［3］蔡强，鲍福廷.基于ACIS几何造型平台的固体火箭发动机装药设计［J］.固体火箭技术，2008，31(3).

［4］孙家广，等.计算机图形学(第三版)［M］.北京:清华大学出版社，1998:358-373.

［5］詹海生，李广鑫，马志欣.基于ACIS的几何造型技术与系统开发［M］.北京:清华大学出版社，2002:30-49.

［6］LEDAS Ltd.LEDAS geometric solver 2D overview［R］.http://ledas.com/products/lgs2d.

［7］Püskülcü G，Ulas A.3-D grain burnback analysis of solid propellant rocket motors:Part 2—modeling and simulations［J］.Aerospace Science and Technology，2008(12):585-591.