陳曉航，馬娟，郭華，李雄飛

（1.西安航天發(fā)動機有限公司，陜西西安，710100；2.西安理工大學藝術(shù)與設(shè)計學院，陜西西安，710048；3.陜西師范大學計算機科學學院，陜西西安，710062）

0 引言

我國探月工程二期月球著陸探測器采用具有大范圍變推力工作能力的液體變軌發(fā)動機作為其主發(fā)動機，用于完成月球著陸探測器中途修正、近月制動、動力下降等軟著陸任務。發(fā)動機的推力由燃料燃燒產(chǎn)生的噴流與噴管內(nèi)型面作用而產(chǎn)生，推力線理論上為噴管內(nèi)型面的幾何中心線即為喉部橫截面中心到噴管出口中心的連線。該型發(fā)動機在熱試車前及交付總體前要進行推力線參數(shù)測量，為總體裝配確定發(fā)動機推力線在艙體坐標系下的位置提供坐標數(shù)據(jù)。一般來說，火箭發(fā)動機推力線參數(shù)要求推力線與對接坐標的偏移量、偏斜量，但是探月系列發(fā)動機推力線參數(shù)要求較為特殊，要求參數(shù)較多，有些參數(shù)商業(yè)軟件無法直接給出結(jié)果，需要檢測人員自行根據(jù)實際的推力方向進行換算，費時費力，存在質(zhì)量隱患，本文通過轉(zhuǎn)換思路，從基礎(chǔ)入手，建立推力線，坐標系的數(shù)學模型，運用高等數(shù)學，計算機圖形學等方法開發(fā)一套計算軟件進行自動化計算推力線參數(shù)，該方法避免人為出現(xiàn)的計算錯誤，并極大提高了計算效率。

1 推力線參數(shù)及計算現(xiàn)狀

推力線參數(shù)計算前需要確定推力線，推力線一般由熱標試車(實際推力線)或數(shù)據(jù)者發(fā)動機測量(理論推力線)數(shù)據(jù)確定。

1.1 推力線參數(shù)

發(fā)動機不同狀態(tài)的推力線測量，均要提供以下參數(shù)，如圖推力線參數(shù)如圖1 所示。其中坐標系OXYZ(O′X′Y′Z′)為發(fā)動機對接面（精測鏡）坐標系，坐標系O1X1Y1Z1為坐標系OXYZ(O′X′Y′Z′)的平移坐標系。對于對接面坐標系，O1點為推力向量與YOZ 平面的交點，α為推力向量與X1 軸夾角，β為推力向量Y1 軸的夾角，γ為推力向量與Z1 軸的夾角。角度ε為交點O1在平面YOZ 中的方位角。O1B為推力線在平面YOZ 的投影，角度θ為O1B在平面YOZ 中的方位角。其中角度ε、θ較為特殊，其數(shù)值分別為Y 軸逆時針轉(zhuǎn)動到OO1及O1B的角度。同樣對于精測鏡坐標系，O2點，(α2,β2,λ2,ε2,θ2)，O′O2的意義相同，不再贅述。

圖1 坐標系下推力線參數(shù)示意圖

(1)推力線平移量：發(fā)動機噴管喉部中心點在對接面的垂點與對接面理論中心距離；

(2)推力線偏斜量：發(fā)動機推力線與對接面垂線的夾角；

(3)對接面坐標系下推力線與對接面交點空間坐標值O1點坐標；

(4) 發(fā)動機實際推力線與對接面坐標軸的夾角(α1,β1,γ1,ε1,θ1)及OO1距離；

(6)精測鏡坐標系下，推力線與精測鏡坐標系Y′O′Z′平面交點O2點空間坐標值；

(7) 發(fā)動機實際推力線與精測鏡坐標軸的夾角(α2,β2,γ2,ε2,θ2)及O′O2距離。

1.2 推力線兩點坐標獲取

熱標前的發(fā)動機測量：通過測量發(fā)動機實體幾何特征，獲取發(fā)動機噴口內(nèi)圓中心及喉部外部中心在對接面坐標系、精測坐標系下的坐標值。噴口中心和喉部中心為推力線上的兩個點。

熱標試車后的發(fā)動機測量：通過發(fā)動機熱標試車，獲取發(fā)動機推力線部分參數(shù)，熱標數(shù)據(jù)為OO1距離及α1,ε1,θ1角度值。通過空間幾何關(guān)系：

可以得到：β1,γ1值。

根據(jù)空間推力線參數(shù)方程：

令A 為常量定為100?？傻玫絰，y，z 值，從而得到推力線上O1,A兩點坐標值。

1.3 推力線參數(shù)計算

目前推力線參數(shù)計算流程如圖2 所示。

圖2 推力線參數(shù)計算流程示意圖

根據(jù)以上流程就可以計算出推力線參數(shù)，但是在計算ε，θ角時均需要換算，具體換算方法為：B 點為推力線上A 點在YOZ 平面的投影點。O1點在坐標系第一象限和第二象限時，ε=ε′,O1點在坐標系第三象限和第四象限時，ε=360° -ε′，B 點在平移坐標系O1X1Y1Z1第一象限和第二象限時，θ=θ′，B 點在平移坐標系第三象限和第四象限時，θ=360° -θ′，如圖3 所示。

圖3 O1 點在第三象限且B 點在平移坐標系第四象限時ε，θ 換算示意圖

1.4 存在的問題

現(xiàn)有推力線參數(shù)計算方法過于復雜，在計算ε，θ角時，由于商業(yè)軟件無法直接給出所需的值，必須經(jīng)過換算。尤其在換算θ角時，必須在平移坐標系下進行，極易出現(xiàn)換算錯誤。為了盡量避免出現(xiàn)計算錯誤，目前在進行推力線參數(shù)計算時，均是由2 組測量人員共4 人分別計算，效率低下，在這種情況下依然出現(xiàn)多次計算錯誤，有很大的安全隱患。

2 推力線參數(shù)自動化計算實現(xiàn)

推力線參數(shù)計算具有輸入量少，運算復雜，易出錯，效率低等特點，非常適合設(shè)計一款專用計算軟件進行計算。專用軟件僅需要輸入推力線兩點坐標值，對接面坐標系下精測鏡坐標系的軸向和原點數(shù)據(jù)，點擊計算按鈕直接計算出結(jié)果，流程如圖4 所示，軟件界面如圖5 所示。

圖4 推力線參數(shù)自動化計算流程圖

圖5 推力線參數(shù)自動化計算軟件界面

推力線參數(shù)計算過程中用到了點、直線、平面構(gòu)造，坐標系建立，坐標系轉(zhuǎn)換，直線和面求交點，直線投影至平面，直線夾角計算，ε，θ角換算等空間幾何計算。以下詳述幾種重要的算法實現(xiàn)。

2.1 坐標系建立及坐標系轉(zhuǎn)換的實現(xiàn)

笛卡爾坐標系是由三個相互垂直的向量及原點構(gòu)成，但是實際測量中沒有完全垂直的兩個特征，總會有誤差，所以實際測量中的坐標系建立和理論幾何中的坐標系建立有所不同。理論上兩個相互垂直的向量和一個原點才可建立坐標系，實際測量中的特征向量不可能完全垂直，這就造成了兩個不垂直的向量無法建立坐標系，實際工程測量中均是通過算法調(diào)整一個向量使其與另一個向量垂直后再建立坐標系。算法如下。

假定由于用于建立坐標系的兩個向量為i，j。

i與j相互垂直時，直接通過i×j=k即可得到第三個向量。

i與j不垂直時：

坐標系所需的三個向量就變成了i、、k，可以看到我們優(yōu)先保證了?不變，調(diào)整?至?，所以實際測量i、j的選擇順序不同，建立坐標系也不同。

推力線兩點坐標是在對接面坐標系下獲得的，所以在計算精測鏡坐標下的推力線參數(shù)時，需要建立精測精坐標系，之后再將推力線兩點坐標轉(zhuǎn)換到精測精坐標系下，然后再進行計算。按照上述方法建立精測鏡坐標系，得到對接面坐標系下精測鏡坐標軸向量和原點的值分別為[i,j,k]和(xo,yo,zo)，設(shè)某點P 在對接面坐標下坐標值為 (x1,y1,z1)，該點在精測精坐標下坐標為 (x2,y2,z2),其中i,j,k均為3×1矩陣。根據(jù)公式(1)即可得出P 點在精測鏡坐標系下坐標值(x2,y2,z2)。

2.2 直線投影至平面算法實現(xiàn)

由于直線是由一個向量和一個點確定的，向量確定直線方向，點確定直線位置。將直線投影至平面的算法分為兩部實現(xiàn)，將直線向量投影至平面和將直線上的一點投影至平面。

已知直線l的向量為m(il,jm,kn)，平面p法線向量為n(ip,jp,kp)，直線上一點為Pl(xl,yl,zl)，平面上一點為Pp(xp,yp,zp)。

2.3 ε，θ 角計算及換算

直線的夾角計算使用公式(3)計算的。其中m、n為兩個向量，α為兩向量夾角。

由于對接面坐標系及精測鏡坐標系下ε，θ角度計算方法相同，介紹計算方法時不再區(qū)分。點O1坐標可以通過直線與平面求交點計算出，O1B的向量可以通過推力線投影計算得出。如圖6 所示，ε′、θ′角為通過向量直接計算的向量夾角，ε，θ為我們所求角度。Y 軸向量為(0,1,0)，設(shè)O1B向量為(0,j,k)，θ′與θ的換算關(guān)系如式（4）。

圖6 向量夾角與ε，θ 換算關(guān)系示意圖

由該式在設(shè)計計算軟件時，僅需要加入一個判斷k≥ 0時 語句就可以換算出ε，θ。

3 軟件驗證及計算偏差控制

軟件編制完成后，進行了驗證，選取了20 組自定數(shù)據(jù)分別用本軟件和商業(yè)軟件計算，計算結(jié)果對比顯示，軟件輸入數(shù)據(jù)與商業(yè)軟件設(shè)置完全一致時，得出的結(jié)果也完全一致。但是，推力線兩點數(shù)據(jù)以及精測鏡坐標軸數(shù)據(jù)一般都是無理數(shù)，輸入數(shù)據(jù)時，必須進行數(shù)據(jù)處理，即對某位小數(shù)進行四舍五入。在保留7位小數(shù)時，O1點坐標、α1、、β1、γ1、ε1、θ1、OO1、精測鏡坐標系下O2點坐標、OO2等值差異一般在第五位小數(shù)，可以忽略不計。當α2角比較小時，精測鏡軸向向量對角度影響比較大，在保留7位小數(shù)時，當α2角在0°～1°之內(nèi)時，(α2,β2,γ2,ε2,θ2)的差異一般在第一或者第二小數(shù)，在保留10位小數(shù)時當α2角在0°～1°，(α2,β2,γ2,ε2,θ2)，的差異一般在第3 到第6 小數(shù)。雖然對于坐標值及尺寸，保留7 位小數(shù)就足夠精確了，但是對于角度保留小數(shù)位數(shù)越多就越好。

在輸入一致情況下，本軟件的與商業(yè)軟件的計算結(jié)果完全一致。輸入數(shù)據(jù)經(jīng)過處理，就導致本軟件輸入數(shù)據(jù)與商業(yè)軟件數(shù)據(jù)不一致，導致計算結(jié)果有偏差，為了盡量減小這種偏差，應將輸入數(shù)據(jù)保留盡可能多的小數(shù)位。若能滿足使用者或者設(shè)計要求精度下，也可以適當減少保留位數(shù)。

4 軟件實際應用

該軟件于2018 年6 月正式編制完成，為了進一步確認該軟件計算安全可靠,自編制完成之日起，凡是探月著陸器發(fā)動機的推力線計算均由人工利用商業(yè)軟件(簡稱人工計算)和用該軟件分別計算，包括了推力室及發(fā)動機狀態(tài)下推力線參數(shù)計算，共計20 余次。其中有兩次人工計算結(jié)果與該軟件計算結(jié)果不一致，經(jīng)過分析發(fā)現(xiàn)，一次是人工計算在角度換算時出錯，一次是人工計算在構(gòu)造推力線時出錯，該軟件并未出現(xiàn)過計算錯誤。從實際應用情況來看，該軟件確實比人工計算更安全更可靠。

5 結(jié)論

由于該系列發(fā)動機推力線參數(shù)計算輸入?yún)?shù)少，流程復雜，計算困難易出錯，目前方法計算效率低且很難保證準確性，所以基于C#開發(fā)自動化計算軟件是解決此類問題的很好途徑，不僅避免人為出現(xiàn)的計算錯誤而且并極大提高了計算效率。