某型探月發(fā)動(dòng)機(jī)推力線參數(shù)自動(dòng)化計(jì)算方法研究

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

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

0 引言

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

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

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

1.1 推力線參數(shù)

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

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

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

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

(3)對接面坐標(biāo)系下推力線與對接面交點(diǎn)空間坐標(biāo)值O1點(diǎn)坐標(biāo)；

(4) 發(fā)動(dòng)機(jī)實(shí)際推力線與對接面坐標(biāo)軸的夾角(α1,β1,γ1,ε1,θ1)及OO1距離；

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

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

1.2 推力線兩點(diǎn)坐標(biāo)獲取

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

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

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

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

令A(yù) 為常量定為100?？傻玫絰，y，z 值，從而得到推力線上O1,A兩點(diǎn)坐標(biāo)值。

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

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

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

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

圖3 O1 點(diǎn)在第三象限且B 點(diǎn)在平移坐標(biāo)系第四象限時(shí)ε，θ 換算示意圖

1.4 存在的問題

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

2 推力線參數(shù)自動(dòng)化計(jì)算實(shí)現(xiàn)

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

圖4 推力線參數(shù)自動(dòng)化計(jì)算流程圖

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

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

2.1 坐標(biāo)系建立及坐標(biāo)系轉(zhuǎn)換的實(shí)現(xiàn)

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

假定由于用于建立坐標(biāo)系的兩個(gè)向量為i，j。

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

i與j不垂直時(shí)：

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

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

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

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

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

2.3 ε，θ 角計(jì)算及換算

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

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

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

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

3 軟件驗(yàn)證及計(jì)算偏差控制

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

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

4 軟件實(shí)際應(yīng)用

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

5 結(jié)論

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