1、78|电子制作 2023 年 3 月信息工程0 引言 我国探月工程二期月球着陆探测器采用具有大范围变推力工作能力的液体变轨发动机作为其主发动机，用于完成月球着陆探测器中途修正、近月制动、动力下降等软着陆任务。发动机的推力由燃料燃烧产生的喷流与喷管内型面作用而产生，推力线理论上为喷管内型面的几何中心线即为喉部横截面中心到喷管出口中心的连线。该型发动机在热试车前及交付总体前要进行推力线参数测量，为总体装配确定发动机推力线在舱体坐标系下的位置提供坐标数据。一般来说，火箭发动机推力线参数要求推力线与对接坐标的偏移量、偏斜量，但是探月系列发动机推力线参数要求较为特殊，要求参数较多，有些参数商业软件无法直
2、接给出结果，需要检测人员自行根据实际的推力方向进行换算，费时费力，存在质量隐患，本文通过转换思路，从基础入手，建立推力线，坐标系的数学模型，运用高等数学，计算机图形学等方法开发一套计算软件进行自动化计算推力线参数，该方法避免人为出现的计算错误，并极大提高了计算效率。1 推力线参数及计算现状推力线参数计算前需要确定推力线，推力线一般由热标试车(实际推力线)或数据者发动机测量(理论推力线)数据确定。1.1 推力线参数发动机不同状态的推力线测量，均要提供以下参数，如图推力线参数如图 1 所示。其中坐标系()OXYZ O X Y Z为发动机对接面（精测镜）坐标系，坐标系11 11O X YZ为坐标系(
3、)OXYZ O X Y Z的平移坐标系。对于对接面坐标系，1O点为推力向量与 YOZ 平面的交点，为推力向量与 X1 轴夹角，为推力向量 Y1 轴的夹角，为推力向量与 Z1 轴的夹角。角度为交点1O在平面 YOZ 中的方位角。1O B为推力线在平面 YOZ 的投影，角度为1O B在平面 YOZ 中的方位角。其中角度、较为特殊，其数值分别为 Y 轴逆时针转动到1OO及1O B的角度。同样对于精测镜坐标系，2O点，22222(,)，2O O的意义相同，不再赘述。(1)推力线平移量：发动机喷管喉部中心点在对接面的垂点与对接面理论中心距离；(2)推力线偏斜量：发动机推力线与对接面垂线的夹角；(3)对接
4、面坐标系下推力线与对接面交点空间坐标值1O点坐标；(4)发 动 机 实 际 推 力 线 与 对 接 面 坐 标 轴 的 夹 角11111(,)及1OO距离；(5)精测镜坐标系下推力线与对接面交点空间坐标值(,)oooxyz；(6)精测镜坐标系下，推力线与精测镜坐标系Y O Z平面交点2O点空间坐标值；(7)发 动 机 实 际 推 力 线 与 精 测 镜 坐 标 轴 的 夹 角22222(,)及2O O距离。XYZX1Y1Z1OO1推力线AB(X)(Y)(Z)(O)(O2)图 1 坐标系下推力线参数示意图 1.2 推力线两点坐标获取热标前的发动机测量：通过测量发动机实体几何特征，获取发动机喷口内
5、圆中心及喉部外部中心在对接面坐标系、精测坐标系下的坐标值。喷口中心和喉部中心为推力线上的两个点。热标试车后的发动机测量：通过发动机热标试车，获取发动机推力线部分参数，热标数据为1OO距离及111,角度值。通过空间几何关系：111111222111cossinsincossincoscoscoscos1=+=某型探月发动机推力线参数自动化计算方法研究 陈晓航1，马娟1，郭华2，李雄飞3（1.西安航天发动机有限公司，陕西西安，710100；2.西安理工大学艺术与设计学院，陕西西安，710048；3.陕西师范大学计算机科学学院，陕西西安，710062）摘要：探月系列发动机推力线参数要求较为特殊，有些
6、参数商业软件无法直接给出结果，需要检测人员自行根据实际的推力方向进行换算，费时费力，存在质量隐患，本文通过建立推力线，坐标系的数学模型，基于C#语言编写推力线参数计算需要的空间几何算法，实现了推力线参数自动化计算。关键词：发动机；推力线；自动化计算DOI:10.16589/11-3571/|79信息工程可以得到：11,值。根据空间推力线参数方程：010101coscoscosxxAyyAzzA=+=+=+其中：(,)oooxy z为1O点坐标。00110110cossinxyOOzOO=令 A 为常量定为 100。可得到 x，y，z 值，从而得到推力线上1O,A两点坐标值。1.3 推力线参数计
7、算目前推力线参数计算流程如图 2 所示。根据以上流程就可以计算出推力线参数，但是在计算，角时均需要换算，具体换算方法为：B 点为推力线上A 点在 YOZ 平面的投影点。O1点在坐标系第一象限和第二象限时，=,O1点在坐标系第三象限和第四象限时，360=，B 点在平移坐标系11 11O X YZ第一象限和第二象限时，=，B 点在平移坐标系第三象限和第四象限时，360=，如图 3 所示。1.4 存在的问题现有推力线参数计算方法过于复杂，在计算，角时，由于商业软件无法直接给出所需的值，必须经过换算。尤其在换算 角时，必须在平移坐标系下进行，极易出现换算错误。为了尽量避免出现计算错误，目前在进行推力线
8、参数计算时，均是由 2 组测量人员共 4 人分别计算，效率低下，在这种情况下依然出现多次计算错误，有很大的安全隐患。YZOO1BY1Z1YZOO1BY1Z1实际需要的,角商业软件算出的的,角图 3 O1点在第三象限且 B 点在平移坐标系第四象限时，换算示意图2 推力线参数自动化计算实现推力线参数计算具有输入量少，运算复杂，易出错，效率低等特点，非常适合设计一款专用计算软件进行计算。专用软件仅需要输入推力线两点坐标值，对接面坐标系下精测镜坐标系的轴向和原点数据，点击计算按钮直接计算出结果，流程如图 4 所示，软件界面如图 5 所示。热标试车数据实体推力测量点击计算按钮专业软件内输入对接面坐标系下
9、精测镜坐标系轴向及原点数据计算推力线两点坐标专用软件内输入推力线两点坐标结束 图 4 推力线参数自动化计算流程图测量对接面，建立对接面坐标系，抽取坐标轴将推力矢量两点坐标构造推力线构造OO1直线计算并换算1测量精测镜，建立精测镜坐标系，并抽取坐标轴计算O1点在精测镜坐标下坐标值计算1，1，1将推力线投影至YOZ平面，计算并换算1计算2，2，2构造交点O1，并计算OO1构造交点O2，并计算OO2构造OO2直线计算并换算2将推力线投影至YOZ平面，计算并换算2图 2 推力线参数计算流程示意图80|电子制作 2023 年 3 月信息工程推力线参数计算过程中用到了点、直线、平面构造，坐标系建立，坐标系
10、转换，直线和面求交点，直线投影至平面，直线夹角计算，角换算等空间几何计算。以下详述几种重要的算法实现。图 5 推力线参数自动化计算软件界面 2.1 坐标系建立及坐标系转换的实现笛卡尔坐标系是由三个相互垂直的向量及原点构成，但是实际测量中没有完全垂直的两个特征，总会有误差，所以实际测量中的坐标系建立和理论几何中的坐标系建立有所不同。理论上两个相互垂直的向量和一个原点才可建立坐标系，实际测量中的特征向量不可能完全垂直，这就造成了两个不垂直的向量无法建立坐标系，实际工程测量中均是通过算法调整一个向量使其与另一个向量垂直后再建立坐标系。算法如下。假定由于用于建立坐标系的两个向量为i，j。i与j相互垂直
11、时，直接通过ijk=即可得到第三个向量。i与j不垂直时：ijk=向量k?同时垂直于i、jikj=i、j、k三个向量相互垂直坐标系所需的三个向量就变成了i、j、k，可以看到我们优先保证了i?不变，调整j?至 j?，所以实际测量i、j的选择顺序不同，建立坐标系也不同。推力线两点坐标是在对接面坐标系下获得的，所以在计算精测镜坐标下的推力线参数时，需要建立精测精坐标系，之后再将推力线两点坐标转换到精测精坐标系下，然后再进行计算。按照上述方法建立精测镜坐标系，得到对接面坐标系下精测镜坐标轴向量和原点的值分别为,i j k和(,)oooxy z，设某点 P 在对接面坐标下坐标值为111(,)x y z，该
12、点在精测精坐标下坐标为222(,)xy z,其中,i j k均为 31矩阵。根据公式(1)即可得出 P 点在精测镜坐标系下坐标值222(,)xy z。222101010(,)(,),xy zxxyy zzi j k=式(1)2.2 直线投影至平面算法实现由于直线是由一个向量和一个点确定的，向量确定直线方向，点确定直线位置。将直线投影至平面的算法分为两部实现，将直线向量投影至平面和将直线上的一点投影至平面。已知直线l的向量为m(,)lmni jk，平面p法线向量为n(,)pppijk，直线上一点为lP(,)lllx y z，平面上一点为pP(,)pppxyz。mns=向量s同时垂直于m、n ns
13、m=m为直线投影后的向量设点lP(,)lllx y z投影至平面上的点为lP(,)lllxyz，直线llPP垂直于平面p,则llPP向量(,)llllllxxyyzz平行于平面p法线向量为n，将lP(,)lllxyz在平面p上，将其坐标带入平面一般方程可得方程组：()()()0llllllllllpplpplppxxyyzzijkxxiyyjzzk=|+=式(2)解该方程组即可得到lP(,)lllxyz。投影后的点和投影后的向量确定了投影后的直线。2.3 ，角计算及换算直线的夹角计算使用公式(3)计算的。其中m、n为两个向量，为两向量夹角。cosarccos()m nm n=式(3)由于对接面
14、坐标系及精测镜坐标系下，角度计算方法相同，介绍计算方法时不再区分。点1O坐标可以通过直线与平面求交点计算出，1O B的向量可以通过推力线投影计算得出。如图 6 所示，、角为通过向量直接计算的向量夹角，为我们所求角度。Y 轴向量为(0,1,0)，设1O B向量为(0,j,k)，与的换算关系如式（4）。0,0,360kk=时时 式(4)由该式在设计计算软件时，仅需要加入一个判断0k 时语句就可以换算出，。YZBYZOBOO1向量指向第三象限O1B向量指向第四象限O1O1Y1Y1O()OO1向量指向第三象限O1B向量指向第二象限图 6 向量夹角与，换算关系示意图|81信息工程3 软件验证及计算偏差控
15、制软件编制完成后，进行了验证，选取了 20 组自定数据分别用本软件和商业软件计算，计算结果对比显示，软件输入数据与商业软件设置完全一致时，得出的结果也完全一致。但是，推力线两点数据以及精测镜坐标轴数据一般都是无理数，输入数据时，必须进行数据处理，即对某位小数进行四舍五入。在保留 7 位小数时，1O点坐标、11111、11111、11111、11111、11111、1OO、精测镜坐标系下2O点坐标、2OO等值差异一般在第五位小数，可以忽略不计。当2角比较小时，精测镜轴向向量对角度影响比较大，在保留7位小数时，当2角在01之内时，22222(,)的差异一般在第一或者第二小数，在保留 10 位小数时
16、当2角在 0 1，22222(,)，的差异一般在第 3 到第 6 小数。虽然对于坐标值及尺寸，保留 7 位小数就足够精确了，但是对于角度保留小数位数越多就越好。在输入一致情况下，本软件的与商业软件的计算结果完全一致。输入数据经过处理，就导致本软件输入数据与商业软件数据不一致，导致计算结果有偏差，为了尽量减小这种偏差，应将输入数据保留尽可能多的小数位。若能满足使用者或者设计要求精度下，也可以适当减少保留位数。4 软件实际应用该软件于 2018 年 6 月正式编制完成，为了进一步确认该软件计算安全可靠,自编制完成之日起，凡是探月着陆器发动机的推力线计算均由人工利用商业软件(简称人工计算)和用该软件分别计算，包括了推力室及发动机状态下推力线参数计算，共计 20 余次。其中有两次人工计算结果与该软件计算结果不一致，经过分析发现，一次是人工计算在角度换算时出错，一次是人工计算在构造推力线时出错，该软件并未出现过计算错误。从实际应用情况来看，该软件确实比人工计算更安全更可靠。5 结论由于该系列发动机推力线参数计算输入参数少，流程复杂，计算困难易出错，目前方法计算效率低且很难保证准确性，所以基于 C