某型探月發動機推力線參數自動化計算方法研究

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

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

0 引言

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

1 推力線參數及計算現狀

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

1.1 推力線參數

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

圖1 坐標系下推力線參數示意圖

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

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

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

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

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

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

1.2 推力線兩點坐標獲取

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

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

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

根據空間推力線參數方程：

令A 為常量定為100。可得到x，y，z 值，從而得到推力線上O1,A兩點坐標值。

1.3 推力線參數計算

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

圖2 推力線參數計算流程示意圖

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

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

1.4 存在的問題

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

2 推力線參數自動化計算實現

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

圖4 推力線參數自動化計算流程圖

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

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

2.1 坐標系建立及坐標系轉換的實現

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

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

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

i與j不垂直時：

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

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

2.2 直線投影至平面算法實現

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

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

2.3 ε，θ 角計算及換算

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

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

圖6 向量夾角與ε，θ 換算關系示意圖

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

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

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

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

4 軟件實際應用

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

5 結論

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