減阻增程彈道的射程估算與特性分析

凌王輝，鮮 勇，郭瑋林，李 杰，張大巧

(火箭軍工程大學七系，西安 710025)

減阻增程彈道的射程估算與特性分析

凌王輝，鮮 勇，郭瑋林，李 杰，張大巧

(火箭軍工程大學七系，西安 710025)

為實現彈道導彈射程的快速估算和減阻設計后彈道增程的定量分析，提出一種基于導彈基本參數的射程快速精確估算方法。通過建立減阻模型，仿真計算得到氣動阻力系數表，并對速度計算公式進行逐項積分，利用高度近似值和氣動系數擬合結果得到關機點狀態量和射程的初步值，最終通過迭代計算得到滿足精度要求的射程結果。仿真結果表明，經過高度和氣動系數的迭代計算能有效快速得到滿足精度的計算結果，通過減阻設計，導彈射程提高了162.36 km，增程的效果明顯。

減阻設計；彈道增程；射程快速估算；氣動擬合；迭代計算

0 引 言

導彈的射程決定了導彈的火力攻擊范圍，是衡量導彈性能的重要參數之一。導彈擁有更遠的射程意味著在相同射程下，能投擲更大質量的戰斗部，具有更多的能量進行機動突防[1-2]。自彈道導彈誕生以來，推進劑的推陳出新、材料的更新換代和外形結構的改型升級使導彈的射程得到不斷提升，除此之外，通過外形減阻設計也可以達到提高射程的目的。美國的“三叉戟”Ⅰ型導彈就運用了激波桿設計，起到了減小大約52%阻力，增加550 km射程的效果。在此基礎之上，其改進的“三叉戟”Ⅱ型導彈進一步擴大了頭部的容積，增加了彈頭的投擲數量[3]。

仿真計算和風洞試驗表明，導彈的鈍頭體構型相比其余尖頭體構型，受到的氣動阻力最大[4]，加裝激波桿能改變頭部流場形態，形成斜激波，減小波后壓力與焓值，有效減小氣動阻力[5]。目前激波桿減阻的相關研究主要是實現減阻效果的仿真校驗與外形機構的優化設計[6-7]，但是對減阻設計后彈道增程效果的分析研究卻很少。

目前利用彈道基本參數實現射程估算的研究較為廣泛，如文獻[8]對模型和數據進行分離，設計了通用彈道仿真模型，實現了仿真多種彈道導彈彈道的功能，但此方法仍基于微分方程組的數值計算，計算量大，運算時間長；文獻[9]基于固定類型導彈的推力加速度模板對導彈主動段運動特性進行了仿真計算，估算性能良好，但構建的近常速運動模型忽略了氣動力，無法用于減阻與未減阻彈道的仿真對比。

而文獻[10-11]針對助推滑翔式高超聲速飛行器，根據其運動方程和飛行條件推導了估算公式，利用彈道參數對滑翔彈道進行了解析估算，并引入高度變化特性進一步提高了估算精度。基于此想法，本文以減阻設計的彈道導彈彈道為研究對象，根據飛行過程中的速度微分方程，推導了關機點狀態的計算公式，利用迭代計算，進一步提高了射程的計算精度，并針對具體減阻構型，實現了減阻增程效果的定量分析。

1 減阻原理

由圖1所示的流場示意圖，鈍頭錐體前為弓形激波，激波桿鈍頭體的流場結構則由斜激波、分離區、剪切層和再附激波組成。

鈍頭體結構前，氣流經過強烈壓縮，形成弓形激波，與之相比，減阻桿破壞了弓形激波，形成斜激波。斜激波后的氣流情況與弓形激波后的氣流情況差異較大，因此導致二者頭部的壓力分布不同。根據普朗特關系式，斜激波后的氣壓遠小于弓形激波后的氣壓，這使帶減阻桿外形的頭部壓力遠小于鈍頭錐體頭部壓力，因此減小了帶激波桿外形的氣動阻力。

2 仿真模型與網格

根據外形的具體設計參數，建立加裝激波桿導彈的實體模型。帶激波桿外形的模型如圖2所示。

進行計算區域離散化，將連續的計算區域劃分為多個子區域，確定每個區域中的節點，從而生成網格。本文采用O型的非結構網格，生成的網格如圖3～4所示，其中圖3為生成的總體網格，圖4為物面處的網格，放大了非結構網格在物面處的細節，貼體性較好，密集均勻，能滿足捕捉流場細節的需要。

3 估算方法

假設導彈在射面內飛行，如圖5所示，則導彈的主動段質心運動方程為：

(1)

式中：V為導彈的速度，P為推力，m為飛行過程中的實時質量，g為重力加速度，θ為彈道傾角，X為氣動阻力。

推力P和發動機有效噴氣速度ue為：

(2)

式中：u為發動機燃氣相對彈體的噴射速度，Sa為噴管截面積，pa為燃氣靜壓，p為大氣靜壓。

將式(2)代入式(1)，可得：

(3)

式(3)的兩端關于時間τ進行積分可得：

(4)

當積分上限為tk時，即可求得導彈在關機點的速度Vk：

(5)

式中：Vu k為真空環境、不計重力條件下導彈主動段關機點的速度，ΔV1k為到關機點重力造成的速度損失，ΔV2k為到關機點氣動阻力造成的速度損失，ΔV3k為到關機點大氣壓力造成的速度損失。

(6)

導彈飛行過程中質量與初始質量的比值μ=m/m0，則：

(7)

將式(7)代入式(5)，可得：

(8)

式中：Sm為最大橫截面積，Cx為阻力系數，ρ為大氣密度。

表1 各個變量的數量級(國際單位制)Table 1 The order of magnitude of the variable (International system of units)

由表1可知，Vu k的數量級為103，而ΔV1k數量級為102且接近103。在整個飛行過程中，同一高度下，相對密度ρ/ρ0與相對壓強p/p0的數量級基本一致，可得：

(9)

除與ΔV3k數量級相同的積分項，ΔV2k的其余積分項CxV2變化范圍為0～106，在積分區域的面積數量級為105，Sm與Sa數量級一致，而ΔV3k的常數項p0大小為101325，所以ΔV2k與ΔV3k數量級一致，二者近似相等。ΔV3k的常數項數量級為102，p/p0成指數下降趨勢，積分項明顯小于1，由此可知，Vu k>ΔV1k>ΔV3k≈ΔV2k。

3.1Vu k的估算

直接積分可得Vu k：

(10)

式(10)即為齊奧爾科夫斯基公式，用于計算導彈主動段關機點的理想速度。根據現有發動機工藝水平確定ue和μk，進而估算Vu k的大小。

3.2ΔV1k的估算

將重力加速度g近似為地面重力加速度g0，對sinθ進行積分即可得到ΔV1k的近似值。因為導彈在主動段以小攻角飛行，根據θ=φ-α，將彈道傾角近似為飛行程序角φ。導彈飛行程序分為垂直起飛段、程序轉彎段、瞄準段三段，根據典型的彈道飛行程序，取θ(μ)為：

(11)

(12)

其中,θk<π/2且為定值，則積分下限隨μ0增大而增大，滿足0≤b(μ0-0.45)≤0.5b。被積表達式為超越方程，無法得到其積分表達式，需對三角函數級數展開，可得：

(13)

根據交錯級數的收斂定理可知級數必然收斂，被積函數級數展開滿足一致收斂性，當θk<π/2時，令μ0=0.45，利用式(13)計算可得：

(14)

3.3ΔV2k的估算

根據Vu k和ΔV1k的估算公式，可得到主動段飛行過程中任意質量比μ下的理想速度Vu和重力造成的速度損失ΔV1。利用Vu和ΔV1近似表示飛行過程中的速度V，可得到導彈在發射坐標系下y軸的位移大小。根據圖5所示，x軸位移遠遠小于地球半徑，則tanβk≈sinβk，可利用位移y近似表示高度：

(15)

因為當地聲速根據高度確定，計算得到h即可得到導彈飛行馬赫數和當地空氣密度。高度和速度是時間的函數，通過式(7)可知,高度和速度可看作是質量比μ的函數。根據高度和馬赫數插值求出整個飛行過程中的氣動阻力系數，通過擬合近似也可將阻力系數看作質量比的函數，即Cx(μ)。取步長Δμ，利用梯形積分公式即可得到ΔV2k：

(16)

3.4ΔV3k的估算

與估算ΔV2k的方法相似，將彈道傾角看作μ的函數，根據Vu k、ΔV1k和ΔV2k的估算值，近似表示飛行過程中的速度V，利用速度計算h：

ΔV1(1-nΔμ)-ΔV2(1-nΔμ))sin(θ(1-nΔμ))

(17)

根據高度得到當地壓強計算ΔV3k：

(18)

根據主動段關機點參數估算導彈全射程L的公式[12]，可利用計算得到的Vu k、ΔV1k、ΔV2k、ΔV3k和h估算全射程，即可實現利用導彈的有效噴氣速度、最大橫截面積、氣動參數等易于估算的基本參數對減阻增程結果的定量分析。其中有效噴氣速度、關機點質量比等參數由目前的工藝水平所確定，最大橫截面積、氣動系數等參數由導彈的外形所確定。

(19)

式中：R為地球半徑,βk為主動段射程角,βc為被動段射程角,Θk為關機點速度傾角,υk為能量參數。

3.5射程的迭代計算

在對高度和ΔV2k進行計算時，只利用了部分速度損失值，忽略了ΔV2k、ΔV3k的影響，得到的計算結果將影響全射程的計算精度，同時直接影響因氣動阻力造成的速度損失ΔV2k的計算結果，進而影響減阻增程的定量分析結果。

在得到Vu k、ΔV1k、ΔV2k、ΔV3k的基礎上，將結果再次代入進行迭代計算，得到滿足精度要求的射程計算結果。計算流程如圖6所示。

4 仿真結果與分析

仿真計算未減阻構型和加裝激波桿構型的氣動阻力系數如圖7～8所示。兩種構型的氣動阻力系數曲線在跨聲速時變化劇烈，在Ma2～Ma25范圍內曲線平滑。整個過程飛行的高度和馬赫數隨時間增加而單調遞增，可以根據估算得到部分節點的高度和馬赫數插值計算氣動阻力系數，并利用擬合得到的氣動阻力系數表達式計算ΔV2k。

擬合得到的氣動阻力系數如圖11所示，與其他擬合結果相比，分段二階擬合效果最好，利用質量比分段二次函數表示阻力系數的公式為：

(20)

同理，分段擬合減阻構型阻力系數的效果如圖12所示，得到的近似公式為：

(21)

最終計算得到減阻、未減阻構型在關機點的高度、速度和二者的射程如表2所示。對結果進行迭代計算，在射程誤差滿足小于0.01 km的條件下，減阻和未減阻構型的全射程計算分別迭代3和5次得到收斂值，二者精確的關機點狀態量和射程如表2所示。通過迭代結果可以看到將ΔV2k、ΔV3k代入，對高度h、ΔV2k、ΔV3k進行再次計算能有效提高計算精度。

與迭代計算的結果相比，忽略ΔV2k、ΔV3k對高度的影響進行計算，即忽略了氣動力、大氣壓力的影響，根據式(15)可知，計算得到的未減阻和減阻彈道關機點高度一致。二者相比，減阻構型使關機點速度增加52.53 m/s，全射程增加134.02 km。迭代計算結果顯示，減阻構型使關機點的高度增加2.32 km，速度增加61.69 m/s，全射程增加162.36 km。計算結果的比較說明：迭代計算ΔV2k、ΔV3k進一步提高了主動段射程的計算精度，有利于減阻與未減阻彈道全射程的比較與分析。

表2 一次計算、迭代計算得到的關機點狀態量和射程Table 2 State of the shutdown point and range by first calculation and iteration

未減阻構型和減阻構型各個速度迭代結果隨質量比的變化如圖14～15所示。在關機點因高度、氣動阻力、壓強導致速度減少量如表3所示，其中因高度增加而導致速度減小的量占速度減小總量ΔVs的75%左右，采用減阻構型使ΔV2k在ΔVs中的占比由原先的16.29%減小為11.46%。數據表明,ΔV1k= 800.5 m/s、ΔV2k變化范圍為100 ～ 180 m/s、ΔV3k≈ 120 m/s，遠小于速度V≈ 5800 m/s，驗證了估算方法中各個速度數量級的推算結果，說明在對ΔV1k、ΔV2k、ΔV3k進行計算時，可以忽略量與量之間的耦合作用進行近似計算。

整個飛行過程中，兩種構型的飛行速度差隨質量比的變化曲線如圖16所示。在飛行初始階段，速度較低，氣動阻力較小，二者速度差別不大，隨著質量的減小，二者的氣動阻力差變大，速度差也隨之增加，但飛行到高空后，空氣密度減小，二者的氣動阻力急劇減小，速度差也隨之減小。

表3 速度減少量Table 3 Speed reduction

計算出的關機點狀態量和射程如表4所示，得到減阻和未減阻導彈的仿真彈道如圖18所示。其中坐標原點為地心，發射點、目標點和地心所在平面截取地球得到的半圓為圖中所示的地球表面。

表4 射程和關機點狀態量計算結果比較Table 4 Comparison of the range and the state at shutdown point

與數值仿真計算結果相比，本文計算得到的關機點狀態量、射程的偏差在數值仿真結果中所占百分比如表5所示，其中關機點高度的計算結果偏差較大，關機點的狀態量中速度對全射程的影響更大，結果表明,本文方法能獲得較高精度的關機點狀態值和全射程量。

表5 射程和關機點狀態量的計算偏差Table 5 Computing deviation of the range and the state at shutdown point

在相同的運行環境下，利用數值積分算法進行計算，當推力和秒耗量為恒定值時，數值積分需用時26.55 s，當推力和秒耗量為真實值時，增加了推力和秒耗量的差值計算，用時45.42 s。數值積分計算過程中，每經過一個時間步長，需計算當前時間下的推力和秒耗量，就需要進行一次插值計算，大大增加了計算量。在多次迭代計算的情況下，本文的方法只用了1.46 s，說明相比數值積分算法，本文的方法能以較高的計算速度獲得較高精度的關機點狀態和射程估算值。

5 結 論

1)提出適用于未減阻彈道和減阻設計彈道射程快速、高精度的計算方法。根據彈道導彈飛行過程中的速度計算公式，推導了基于導彈基本參數進行射程估算的計算公式，可在彈道設計過程中根據設計的基本參數論證導彈是否滿足射程指標，為參數的優化提供依據，或根據現有工藝水平和具體外形估算外軍導彈的最大射程。

2)利用函數擬合近似計算阻力系數。對原有的氣動阻力系數表通過函數分段擬合的方法計算氣動阻力，保證了因阻力造成的損失速度的計算精度，提高了計算速度。

3)建立ΔV2k和ΔV3k的迭代模型。將近似計算結果作為初值代入迭代模型進行計算，提高了關機點速度的計算精度。經過高度和氣動系數的迭代計算能快速有效得到滿足精度的計算結果。通過減阻設計，整個飛行過程關機點的速度增加61.69 m/s，導彈的射程增加了162.36 km，增程的效果明顯。

[1] 鮮勇, 李少朋, 李振華, 等. 基于梯度粒子群算法的縱橫向機動跳躍彈道設計及優化[J]. 彈道學報, 2015, 27(3):1-6. [Xian Yong, Li Shao-peng, Li Zhen-hua, et al. Design and optimization of vertical wavy and crosswise maneuvering trajectory based on grads particle swarm algorithm [J]. Journal of Ballistics, 2015, 27(3): 1-6.]

[2] 劉炳琪, 鮮勇, 李振華, 等. 基于非連續點火的變射面空間“M”形彈道設計及優化[J]. 固體火箭技術, 2015, 38(5):608-613. [Liu Bing-qi, Xian Yong, Li Zhen-hua, et al. Design and optimization of changeable launching plane ‘M’ maneuvering trajectory based on discontinuous ignition [J]. Journal of Solid Rocket Technology, 2015, 38(5): 608-613.]

[3] Gauer M, Paull A. Numerical investigation of a spiked blunt nose cone at hypersonic speeds [J]. Journal of Spacecraft and Rockets, 2008, 45 (3): 459-471.

[4] 封貝貝, 陳大融, 汪家道, 等. 頭部形狀對超聲速飛行器力學性能影響分析[J]. 飛行力學, 2012, 30(6):537-540. [Feng Bei-bei, Chen Da-rong, Wang Jia-dao, et al. Affect of head shape on flight dynamics of supersonic vehicles [J]. Flight Dynamics, 2012, 30(6): 537-540.]

[5] Mansour K, Khorsandi M. The drag reduction in spherical spiked blunt body [J]. Acta Astronautica, 2014, 99: 92-98.

[6] Tahani M, Karimi M S, Motlagh A M, et al. Numerical investigation of drag and heat reduction in hypersonic spiked blunt bodies [J]. Heat and Mass Transfer, 2013, 49 (10): 1369-1384.

[7] 李永紅, 高川, 唐新武. 激波針氣動特性及外形參數優化研究[J]. 兵工學報, 2016, 37(8):1415-1420. [Li Yong-hong, Gao Chuan, Tang Xin-wu. Drag reduction characteristics and design optimization of spikes [J]. Acta Armamentarii, 2016, 37(8): 1415 -1420.]

[8] 王海峰, 趙久奮, 王偉. 面向任務級彈道導彈的通用彈道設計仿真研究[J]. 計算機仿真, 2016, 33(6):64-68. [Wang Hai-feng, Zhao Jiu-fen, Wang Wei. Research on general trajectory design and simulation for mission level ballistic missile [J]. Computer Simulation, 2016, 33(6): 64-68.]

[9] 張濤, 安瑋, 周一宇. 基于推力加速度模板的主動段彈道跟蹤方法[J]. 宇航學報, 2006, 27(3):385-389. [Zhang Tao, An Wei, Zhou Yi-yu. The trajectory tracking method in boost phase using thrust acceleration profile [J]. Journal of Astronautics, 2006, 27(3): 385-389.]

[10] 王潔瑤, 江涌, 鐘世勇, 等. 高超聲速遠程導彈彈道解析估算與特性分析[J]. 宇航學報，2016, 37(4):403-410. [Wang Jie-yao, Jiang Yong, Zhong Shi-yong, et al. Analytical estimation and analysis of trajectory performance for hypersonic long-range missiles [J]. Journal of Astronautics, 2016, 37(4): 403-410.]

[11] 王潔瑤, 江涌, 鐘世勇. 助推-滑翔彈道高精度滑翔射程解析估算方法[J]. 宇航學報, 2016, 37(5): 519-525. [Wang Jie-yao, Jiang Yong, Zhong Shi-yong. A high accuracy analytical estimation method for glide range of the boost-glide trajectories [J]. Journal of Astronautics, 2016, 37(5): 519-525.]

[12] 張毅, 肖龍旭, 王順宏. 彈道導彈彈道學[M]. 長沙:國防科技大學出版社,2005:183-189.

RangeEstimationandAnalysisforExtendedTrajectoryBasedonDragReductionDesign

LING Wang-hui, XIAN Yong, GUO Wei-lin, LI Jie, ZHANG Da-qiao

(7th Department, Rocket Force University of Engineering, Xi’an 710025, China)

In order to realize the quantitative analysis of the trajectory and the evaluation of the extend range by using the drag reduction design, a fast and accurate estimation method is proposed based on the basic parameters of a missile. By establishing a drag reduction model, an aerodynamic drag coefficient table is calculated. Then the speed calculation formula is integrated item by item. The initial value of the state at the shutdown point and the range are obtained by the height approximation and the fitting aerodynamic coefficient. After the iteration, the result satisfies the accuracy requirements of the range. The simulation results show that the heuristic calculation of the height and the aerodynamic coefficients can effectively and quickly get the satisfying results. Under the drag reduction design, the final missile range is increased by 162.36 km and the extended-range effect is obvious.

Drag reduction design; Trajectory extension; Range rapid estimation; Aerodynamic coefficient fit; Iteration

TP 731

A

1000-1328(2017)10- 1048- 09

10.3873/j.issn.1000-1328.2017.10.005

2017- 05- 24;

2017- 07- 19

凌王輝(1992-)，男，碩士生，主要從事飛行器設計、制導方面的研究。

通信地址:陜西省西安市灞橋區同心路2號(710025)

電話:15399420435

E-mail:ling_wanghui@163.com