張 驍,劉丙杰,王瑞臣
(海軍潛艇學院,山東 青島 266199)
在火箭或彈道導彈的發(fā)射過程中,其各級火箭發(fā)動機在燃料耗盡或者到達預定位置后會與本體分離,在自身分離速度、地球引力和空氣阻力等多種復雜因素的共同作用下,經(jīng)過一定時間后以高速沖擊地面。隨著航天和導彈技術的發(fā)展應用,產(chǎn)生的發(fā)動機殘骸越來越多,對地面人員和設施帶來的威脅也越來越大。當前各火箭發(fā)射基地在選址過程中均對發(fā)動機殘骸的落點位置進行了考慮,一般能夠避開人員密集區(qū)域,但仍然需要找到可靠方法計算臨時選址或者機動式發(fā)射前提下的火箭發(fā)動機殘骸落點范圍[1]。
在以往的發(fā)射任務中,指揮人員可以使用簡易模型與經(jīng)驗公式對殘骸落點進行預測。但隨著實際任務數(shù)量的增大、火箭射程的增加、目標點位的變化、發(fā)射機動范圍的擴大、發(fā)射條件的放寬,殘骸可能的落點范圍越發(fā)難以預測,給任務負責人員帶來了較大的決策壓力。因此為進一步減小對地面人員和設施的威脅,避免發(fā)生嚴事故,減小后續(xù)殘骸回收的難度,需探索建立一種更加精準、可靠、方便的殘骸落點計算方法[2]。
目前已有相關文獻闡述殘骸落點計算的方法,但算法側重點和面向的任務特點各不相同。文獻[3]提出了基于實時定位的火箭殘骸落點計算模型,精度較高,但需要殘骸具有自主定位和信號發(fā)射的能力,無法應用于當前大多數(shù)火箭發(fā)射任務;文獻[4]提出了非定常 CFD 數(shù)值方法,主要針對具有復雜形態(tài)的大型殘骸,特點是能夠?qū)崿F(xiàn)了對目標氣動發(fā)熱和分離解體過程的模擬,但在落點計算的結果上未能明顯降低誤差;文獻[5]提出了容積卡爾曼濾波(CKF)和時變自回歸(TVAR)模型融合的助推器落點預測方法,在誤差上比單一時間序列的預測方法要小,但精度仍不如數(shù)值積分方法。
本研究基于發(fā)射任務單位的實際信息獲取能力和精度要求,以能夠獲得殘骸與火箭分離點在發(fā)射坐標系中的具體坐標、各速度分量、發(fā)射點的經(jīng)、緯度、高程、射向、殘骸質(zhì)量(假定不變)和有效參考面積為前提,構建一種新的殘骸落點計算方法,用以解決任務單位的實際問題。
火箭發(fā)動機殘骸在與火箭主體分離后不再受主動力的影響,將完全處于被動段飛行,除地球引力、柯氏慣性力和牽引慣性力等必須考慮的因素外,空氣阻力在其被動段飛行的后半部分具有相當大的影響,是殘骸落地速度和位置的重要因素。在實際過程中,各力對殘骸的具體作用與殘骸的初始狀態(tài)都是十分復雜的,為兼顧模型精確度和簡潔實用性,首先對殘骸受力過程進行適當?shù)暮喕痆6-7]。
1.1.1 升降兩端等高度殘骸速度模型
如圖1所示,假設殘骸由點q運動至相同高度的對應點q′的情況下,根據(jù)理論力學相關知識可列公式:
(1)
式中,F為從q點到q′的所有外力,即引力mg和空氣阻力X的和,通過對作用距離s進行積分得到外力F對殘骸所做總功,其值應當?shù)扔趶膓到q′的動能變化量。由于q和q′位于等高度,即兩點地心矢徑的數(shù)值相同,因此引力勢能相同,即在此過程中引力做功總量為零。而空氣阻力X的方向一定與速度方向相反,對殘骸做功只能是負值。因此引力mg和空氣阻力X在這一過程中對殘骸的總做功為負值,即公式(1)右邊的積分一定為負值,由此可推出Vg′ 圖1 殘骸受力分析圖 1.1.2 殘骸速度極值點位置分析 在不考慮空氣阻力的情況下,只有引力對殘骸做功,而彈道頂點為引力勢能最大的點,因此動能一定最小,即速度最小。但考慮空氣阻力后,殘骸速度變化就復雜很多,不僅最小值點將很難確定,而且會出現(xiàn)包括速度極大值點在內(nèi)的多個極值點,但無論如何,極值點的切向加速度等于零為必要條件,即dv/dt=0。 殘骸在其運動的切線方向只受兩個力的影響,空氣阻力X和引力mg在切線方向的分量mgsinθ(θ為速度與水平方向的夾角)。切向加速度為零則意味著X與mgsinθ大小相等方向相反。根據(jù)圖1可知,在點q時的X和mgsinθ的方向均與速度方向相反,殘骸一直處于減速狀態(tài),只有在到達點q′時X和mgsinθ才可能出現(xiàn)方向相反的情況,因此速度極值點必然僅存在于下降段[8]。 一般情況下存在兩個殘骸速度極值點:點一是速度極小值點,位于殘骸剛經(jīng)過彈道頂點的時候,此時由于空氣極為稀薄,空氣阻力X的值與引力mg相比很小,因此只需要很小的sinθ值即可達到切向加速度為零的條件,即剛過彈道頂點時速度向下且與水平方向成某一小角度時。點二是速度極大值點,位于殘骸以極高速度進入稠密大氣層的時候,此時的空氣阻力X受殘骸速度和大氣密度的雙重影響而極速增大,直至在某一位置與mgsinθ達到平衡,即此時的殘骸雖然處于彈道的下降段,但停止加速,轉(zhuǎn)為減速。 上述兩個速度極值點一定存在,但在實際中很可能不僅僅有兩個極值點。殘骸在進入稠密大氣層后的減速階段后期,sinθ趨近于-1,即彈道接近垂直,在殘骸速度的減小的同時,大氣密度卻不再有明顯變化,空氣阻力X逐漸減小后與引力mgsinθ緩慢趨向于平衡,加速度趨近于0,可視為第三個速度極值點。如果考慮實際更復雜的情況,火箭發(fā)動機殘骸一般具有較為復雜且變化的氣動外形,其等效參考面積較大,更容易受到空氣阻力的影響,甚至在下降段可能存在空氣抬升力和下降力的復雜影響,因此有更多極值點出現(xiàn)的可能性。 1.1.3 殘骸射程影響因素分析 將式(1)中做功為零的mg去掉并適當變形可得: (2) 殘骸的墜落過程中,其有效面積S和質(zhì)量m并不是一成不變的,在超高速飛行時,空氣對殘骸產(chǎn)生的熱效應和剝離效應將直接影響S與m的值。對這一過程進行量化的最佳方法是建立S、m與殘骸速度V、殘骸所處大氣狀態(tài)C的函數(shù)關系,且必須考慮殘骸材料和形狀的影響。具體是否復雜化,需綜合考量模型使用單位對精度和效率的要求。 了解上述殘骸彈道特性后有助于我們在建立落點預測模型的過程中的整體把握,減少明顯理論錯誤發(fā)生的可能性,且有助于后期定性的驗證模型建立的準確性。 為充分說明典型火箭發(fā)動機殘骸的墜落過程,根據(jù)某次試射任務的觀測數(shù)據(jù),制作圖2~4所示。 圖2 V-t圖 圖3 H-t圖 圖4 sinθ-t圖 上述主要運動參數(shù)的分析根據(jù)分離點的高低、分離速度的大小和方向改變而稍有變化,但仍可在后續(xù)模型建立、修改和驗證的過程中起到有效指導作用。 殘骸再入大氣后的受力情況十分復雜,不僅受地球引力作用,還受空氣阻力和空氣動力矩的影響,在氣動加熱的情況下,還會使殘骸質(zhì)量、外形發(fā)生較大變化,這些都會在一定程度上影響殘骸落點。為使問題簡化,此處不考慮氣動加熱和殘骸外形的影響。殘骸受分離力的影響,還存在一定程度的翻滾,在10 km高度以上時由于空氣阻力微弱,其翻滾造成的影響可以忽略不計,在進入稠密大氣后殘骸的翻滾會很快受到阻力限制而趨向以某有效面積與速度方向保持穩(wěn)定角度。因此認為全程忽略空氣動力矩與翻滾的影響,將殘骸視為一個質(zhì)點仍能夠保證模型的有效精確。 殘骸無動力段的運動形式與火箭被動段的運動相似,在模型建立過程中可充分借鑒遠程火箭彈道學中關于被動段的成熟理論。在將殘骸看成一個質(zhì)點的情況下,殘骸主要受力為地心引力與空氣阻力?;鸺c殘骸分離點較高,飛行時間較長,不能忽視因地球轉(zhuǎn)動所產(chǎn)生的牽連慣性力和柯氏慣性力的影響。根據(jù)動力學基本理論,在發(fā)射坐標系中可得到火箭殘骸無動力段的矢量表達式為: (3) 為方便計算,在對彈道進行研究的時候一般使用發(fā)射坐標系,以發(fā)射點為原點,以射面與水平面的交線為x軸,以過原點并與發(fā)射點位置的水平面相垂直的軸為y軸,z軸與x軸和y軸成右手螺旋關系。式(3)在發(fā)射坐標系中的投影形式為: (4) 又根據(jù)運動學可得方程: (5) (6) 在抬升力、側向力和空氣動力矩忽略不計后,空氣給殘骸帶來的影響只??諝庾枇。該阻力方向恒定與速度方向相反,在速度坐標系中可表示為: X=CxdtqdtSmdt (7) 式中,Cxdt為殘骸阻力系數(shù),qdt=1/2ρV2為頭部速度頭,ρ為空氣密度,Smdt為殘骸的有效面積(根據(jù)經(jīng)驗可取最大橫截面積)。 即使通過上述簡化后,想要對X的值進行精確計算仍然是一個非常復雜的問題,難以僅通過理論推導和計算來完成。目前比較科學可靠的方法是使用實驗與模型相結合的方法,即通過空氣動力學理論盡可能詳細的建立接近現(xiàn)實變化的數(shù)學模型,再通過空氣動力學實驗的方法對模型中的關鍵系數(shù)與其他要素的變化規(guī)律進行觀察記錄。即使如此也僅能獲得當前實驗條件下的數(shù)據(jù),環(huán)境稍作變化就可能出現(xiàn)不同的情況。開發(fā)人員在火箭的研制過程中必然會根據(jù)其外形和材料進行空氣動力學研究,但由此獲得的圖表和數(shù)據(jù)只能在一定程度上代表火箭本身,無益于此處對發(fā)動機殘骸空氣動力系數(shù)的研究。 Cxdt在殘骸運動的過程中變化最為復雜,在以往案例中通常使用經(jīng)驗數(shù)值來表示。其影響因素很多,主要包括殘骸的外形和材料、殘骸在空中飛行的姿態(tài)、殘骸速度、飛行高度等。在掌握上述影響因素與Cxdt的關系后才可能使用數(shù)學表達式對Cxdt進行近似表示。 殘骸的外形和材料對空氣動力的影響很大,具體作用過程也最為復雜。簡單來說,長細比越大的物體,其阻力系數(shù)越小,表面材料越光滑且凸起越少的物體,其阻力系數(shù)越小。殘骸速度對Cxdt的影響主要通過改變殘骸周邊氣體的可壓縮程度來實現(xiàn),但要找到Cxdt與速度V之間的直接關系是極為復雜的,此處速度需要用馬赫數(shù)ma來表示,音速u與高度H之間本就存在復雜關系,且ma與氣體可壓縮程度的關系也并不是簡單線性,而是根據(jù)研究對象的形狀和材料不同表現(xiàn)出不同的曲線,因此在實際操作中往往需通過查閱實驗數(shù)據(jù)表格獲得。同時Cxdt與大氣的各項基本屬性例如密度、壓力、溫度等均密切相關,而大氣中這些屬性的最主要決定者就是高度H,因此Cxdt與高度H之間存在著復雜關系。綜上可得,在研究對象確定的情況下,對空氣動力系數(shù)Cxdt影響最大的因素便是飛行馬赫數(shù)ma和高度H。 同理,ρ受其他因素影響較小,可認為僅與高度H相關。 鑒于Cxdt的詳細估算十分復雜,因此在以往計算中常使用經(jīng)驗值代替,也有文獻顯示落點計算仍可達到一定精度。但如此則將其過分簡化,必然無法反映出飛行物或殘骸在大氣中的真實運行軌跡[3,6,11],尤其是在殘骸對大氣層的入射角度較小時,有阻力段在全射程中的占比增加,其誤差必然增大。 本研究根據(jù)國際標準大氣表(1976美國)的參數(shù)對空氣密度ρ和不同高度的音速u與幾何高度H的關系進行了探索,使用matlab2013b進行函數(shù)擬合。因空氣密度ρ的數(shù)值從H=0到H=80 km的范圍內(nèi)跨越6個數(shù)量級,導致用一個連續(xù)函數(shù)進行表達時必然會出現(xiàn)誤差平方和極低,但中高空擬合結果與實際數(shù)值出現(xiàn)1至2個數(shù)量級差距的問題。推測此處問題的本質(zhì)原因可能在于地球大氣的復雜分層導致無法使用一個連續(xù)函數(shù)進行合理表示。因此建議采用分段擬合的方式,可得到誤差平方和在10-5級別的表達式ρ=f(H)和u=h(H)。根據(jù)某火箭發(fā)動機殘骸Cxdt系數(shù)與馬赫數(shù)ma的關系表格: 表1 Cxdt與ma的關系 對Cxdt與ma的函數(shù)關系進行擬合,得到Cxdt=g(ma)。聯(lián)立三式可得Cxdt與H、V的函數(shù)關系式: Cxdt=I(H,V) (8) 將式(8)代入式(7)即可得到在速度坐標系中的空氣阻力X,需將其轉(zhuǎn)換至發(fā)射坐標系才能夠代入式(4)進行計算。根據(jù)速度坐標系與發(fā)射坐標系的轉(zhuǎn)換關系可得: (9) 其中:Xi(i=x,y,z)為空氣阻力X在發(fā)射坐標系各軸上的分量,θ是指殘骸速度矢量V與發(fā)射坐標系oxz平面即發(fā)射點所在水平面所成的角度,上為正,下為負,σ是指殘骸速度矢量V與發(fā)射坐標系oxy平面即射面所成的角度,從坐標原點沿發(fā)射坐標系x軸正方向看時,右為正,左為負。由定義可見,角θ變化范圍較大,角σ在主動段由于控制系統(tǒng)的作用,數(shù)值很小,在被動段其對射面的偏離主要來自分離時z方向的初始速度和橫向高空風的影響,在以往建模中一般忽略不計。但正如1.1.3中所作分析,由于殘骸的S與m比值相對火箭本身較大,更容易受偏離力的影響,所以此處不可忽略不計。由θ和σ的定義可推出: (10) 在將地球看做正常地球橢球體時,可通過下列公式計算得到地球引力在發(fā)射坐標系各軸上的分量[12]: (11) (12) (13) (14) B-φs=690.309″/3600/180*π*sin2B (15) (16) (17) (18) 式中,ω、fM、μ均為接近地球?qū)嶋H的總地球橢球體基本數(shù)據(jù),可通過查表獲得,r為殘骸到地心的矢徑,gr為引力加速度在地心坐標系矢徑r上的分量,gω為引力加速度在地球自轉(zhuǎn)角速度矢ω上的分量。R0i(i=x,y,z)為發(fā)射點地心矢徑R0在發(fā)射坐標系各軸上的投影。B為殘骸所處地理緯度,φs為殘骸所處的地心緯度,B0和φs0分別指發(fā)射點的地理緯度和地心緯度,北緯定義為正,南緯定義為負,C=B0-φs0為發(fā)射坐標系中oy軸與矢徑R0的夾角,即發(fā)射點地理緯度與發(fā)射點地心緯度的差值,可通過常用經(jīng)驗式(15)計算得。ωi(i=x,y,z)為地球自轉(zhuǎn)角速度矢ω在發(fā)射坐標系各坐標軸上的分量,B0、A0分別為發(fā)射點地理緯度和發(fā)射方位角[7-8]。 φs的計算方法除使用上述式(14)外,還可采用由發(fā)射坐標系向地心直角坐標系轉(zhuǎn)換,再向地心極坐標系轉(zhuǎn)換的方法。在需要得到地理緯度B時可采用式(15)進行計算。 根據(jù)定義可知,發(fā)射坐標系是固連在發(fā)射原點上的,因此在隨地球相對慣性坐標系做平移運動的前提下,還繞地球兩極做圓周運動。根據(jù)理論力學知識,導彈殘骸在運動過程中始終受柯氏慣性力的影響,其在發(fā)射坐標軸上的加速度分量可分別表示為[7]: (19) 因發(fā)射坐標系隨地球相對慣性參考系運動,因此必須考慮施加于導彈殘骸上的牽引慣性力,其在發(fā)射坐標軸上的加速度分量可分別表示為[7]: (20) 其中: (21) 聯(lián)立上述(4)至(21)式,在給定初始x、y、z和Vx、Vy、Vz的情況下通過對時間積分即可得到殘骸彈道數(shù)據(jù)。該結果數(shù)據(jù)為殘骸在發(fā)射坐標系中的坐標和速度,要得到常用地理經(jīng)、緯度坐標,還需要進行由發(fā)射坐標系至地心直角坐標系,再到地心極坐標系的轉(zhuǎn)換,然后使用式(15)得到地理經(jīng)、緯度。 本研究使用Matlab2013b進行仿真實驗,根據(jù)上述公式進行建模并以固定步長進行積分,直至滿足條件高度H小于等于設定高度H0時停止積分,輸出落點經(jīng)緯度及其它運動參數(shù)。 為評價模型精度,首先可采用實際落點(LMDs,Bs)與預測落點(LMDc,Bc)之間距離Δs來表示誤差的大小。因兩點距離較近,可忽略地球扁率影響得[13]: cos(Bs)cos(Bc)cos(LMDs-LMDc)] (22) 此處Δs的值即為預測值與實際值的絕對誤差,雖然能夠在一定程度上說明模型精度,但僅使用Δs無法全面比較不同落點計算模型之間的精度優(yōu)劣。不同模型所使用的數(shù)據(jù)均來自不同的任務或?qū)嶒灒l(fā)射點和目標點各異,殘骸飛行距離差別巨大,甚至達到上千公里。當殘骸射程相對增大時,落點計算模型所積累的誤差也隨之增大,Δs的絕對值也必然相應增大。因此為將不同射程的實測數(shù)據(jù)用來對比,模型精度的評價參數(shù)除Δs外還應當包括其與殘骸射程的比值Δs/L[14-15]。 使用某火箭位一固定地點向固定方向進行多次實射所采集的數(shù)據(jù)進行驗證,取步長dt=0.01 s,落點高度H0=0結果如表2和圖5所示。 根據(jù)對比計算的結果顯示,在殘骸射程由235.928 km增大至1 469.388 km的五次實驗中,仿真結果的絕對誤差分別為192.373 m、1 081.079 m、3 274.104 m、5 508.066 m、4 292.573 m,與3.1節(jié)中的推論相應和。隨著殘骸射程的增加絕對誤差Δs有明顯的增大,在射程達到1 469.388 km時,誤差最大達到5 508.066 m,誤差百分比達到了0.389 9%。綜合來看,五次仿真的平均相對精度a(Δs/L)=0.220 6%。 考慮到模型建立的主要目的,千分之二左右的誤差能夠滿足預定的落點預測功能,足以實現(xiàn)安全預警和便于回收等一系列目的,能夠為發(fā)射指揮人員的判斷提供可靠依據(jù)。 論證該落點計算模型的精度,不能僅依據(jù)現(xiàn)有數(shù)據(jù),還應當能夠在與同類型預測模型的對比中保持一定的優(yōu)勢地位。目前所閱相關文獻中均對所建模型的精度進行了一定的論證,但達到足夠精度并能夠給出足夠數(shù)據(jù)用以支撐模型之間橫向?qū)Ρ鹊奈墨I卻并不多見。決定采用文獻[3]和文獻[4]中的模型與本落點計算模型進行對比。 表2 精度分析結果 圖5 3號仿真軌跡圖示 文獻[3]中使用火箭殘骸實時定位信息為基礎,在殘骸落地前的40 km至20 km處測量得到殘骸具體運動參數(shù),用以計算殘骸在高度10 km處的位置。型號一火箭殘骸在高度相差10 km時平均偏差為0.63 km,高度相差20 km時平均偏差為0.72 km,高度相差30 km時平均偏差為1.36 km。型號二火箭殘骸在高度相差10 km時平均偏差為0.51 km,高度相差20 km時平均偏差為0.74 km,高度相差30 km時平均偏差為1.27 km。綜合來看該模型能夠?qū)⒔^對誤差控制在1.36 km至0.51 km之間,但相對誤差卻達到4.533%與5.1%之間。該方法雖能夠控制絕對精度,但其只對殘骸落地前的一小部分有空氣阻力彈道進行預測,且需要殘骸的實時定位信息為基礎數(shù)據(jù),要求在殘骸上內(nèi)置定位裝置,對觀測能力或者信息傳輸能力要求較高,適用性受到一定限制。 文獻[4]通過建立火箭不同仰角和側滑角下的氣動數(shù)據(jù)庫,對初始仰角、側滑角與射程、側位移的關系進行了詳細研究,使用非定常CFD方法能夠較為精準的預測火箭再入大氣層后的形變、解體和落點預測。進行了火箭發(fā)動機殘骸自30 km高度隕落的飛行試驗,繪圖結果顯示在射程約52 km的情況下誤差在1 km以內(nèi),但因并未給出具體數(shù)據(jù),所以無法進一步與本文方法比較。該方法預測的絕對誤差較小,但需要占用的計算資源較多,且不具有普適性,主要針對結構較為復雜或形態(tài)巨大的特定目標。 本研究所用方法從分離點的運動參數(shù)入手,能夠?qū)Πl(fā)動機殘骸運動的全過程進行預測,并且不依靠后續(xù)觀測數(shù)據(jù)的矯正,占用資源較少,計算速度較快,絕對誤差控制在任務單位的搜索能力范圍內(nèi),相對誤差較小,能夠適用于結構和材料較為簡單的所有發(fā)動機殘骸落點預測,符合執(zhí)行任務單位的能力要求。 本研究缺少殘骸的真實形狀、材料、剩余燃料等重要數(shù)據(jù),因此對空氣動力的模擬只局限于速度坐標系下沿x軸負方向的阻力,忽略了空氣升力和側向力及其力矩帶來的旋轉(zhuǎn),不能完全反映空氣對殘骸的影響,同時也無法估算出分離點的實際殘骸質(zhì)量和后續(xù)因氣動發(fā)熱造成的質(zhì)量損失[16-17]。限于信息獲取能力有限,代入的初始狀態(tài)為火箭本體與發(fā)動機分離時的運動數(shù)據(jù),忽略了分離過程對殘骸運動狀態(tài)的影響。也未能充分考慮高空風和擾動引力在殘骸近地飛行時對軌跡帶來的擾動[18-19]。后續(xù)還可使用神經(jīng)網(wǎng)絡等智能預測方法對同課題進行研究,嘗試在殘骸基本屬性未知的情況下對落點進行預測[20]。 根據(jù)某火箭實際使用和訓練任務的需求推動,在理論力學的基礎上建立動力學方程組,將地球大氣密度、空氣阻力系數(shù)和音速a與高度H的關系進行擬合,代入方程組得到滿足精度要求的殘骸落點計算模型,以分離點數(shù)據(jù)為殘骸初始狀態(tài),對時間積分即可得到殘骸全射程運動參數(shù),進而實現(xiàn)對落點的準確預測。 通過與任務實際數(shù)據(jù)進行對比,證明該方法誤差在合理范圍內(nèi),能夠滿足指揮員的輔助決策需求。通過與其他方法的橫向?qū)Ρ?,證明該方法在當前任務需求下具有計算速度快、信息需求少及經(jīng)濟性好的優(yōu)勢。后續(xù)仍可對殘骸受力的過程進一步細化,逐步加入分離力、氣動升力和側向力、擾動引力及質(zhì)量損失、形態(tài)變化和分離解體等因素的影響。1.2 殘骸動態(tài)分析
2 殘骸落點計算模型
2.1 殘骸所受空氣阻力X的計算
2.2 殘骸所受地球引力加速度g的計算
2.3 殘骸所受柯氏加速度的計算
2.4 殘骸所受牽連加速度的計算
3 模型精度驗證
3.1 精度分析方法
3.2 精度分析結果
3.3 橫向精度對比
3.4 模型問題分析
4 結束語