牛鈺森,姜 毅,李 靜,陳 苗,史少巖,于邵禎
(1.北京理工大學(xué) 宇航學(xué)院,北京 100081;2.海軍裝備研究院,北京 100161)
?
基于氣動耦合的導(dǎo)彈適配器初始分離彈道數(shù)值分析
牛鈺森1,姜毅1,李靜1,陳苗1,史少巖1,于邵禎2
(1.北京理工大學(xué) 宇航學(xué)院,北京 100081;2.海軍裝備研究院,北京 100161)
摘要:為了準(zhǔn)確預(yù)測與導(dǎo)彈分離后導(dǎo)彈適配器的初始分離彈道,使用基于局部重構(gòu)方法的動網(wǎng)格技術(shù),將適配器的運動與流場的變化耦合。通過CFD方法實時求解適配器飛行過程中的氣動載荷,并對適配器的六自由度運動微分方程進行求解。使用基于尺寸函數(shù)的局部重構(gòu)方法對適配器周圍網(wǎng)格進行更新。通過數(shù)值計算得到了適配器在初始分離彈道階段的氣動載荷系數(shù)曲線,以及適配器的分離速度、角速度曲線和分離飛行軌跡,并與試驗結(jié)果進行了對比,兩者高度一致。結(jié)果表明,在分離初速的作用下,適配器與導(dǎo)彈的距離逐漸增大,在氣動力作用下適配器發(fā)生偏轉(zhuǎn),直到初始分離彈道結(jié)束適配器未與導(dǎo)彈發(fā)生碰撞。
關(guān)鍵詞:適配器;CFD;六自由度;動網(wǎng)格;氣動耦合
安裝在導(dǎo)彈與發(fā)射箱之間的適配器能夠為導(dǎo)彈提供支撐,并且在發(fā)射過程中起到減震與導(dǎo)向的作用。對于安裝在導(dǎo)彈外表面的實心適配器,在導(dǎo)彈發(fā)射過程中隨著導(dǎo)彈一起運動,出箱后通過爆炸螺栓或者彈簧機構(gòu)與導(dǎo)彈解鎖分離。在初始分離彈道階段,適配器的飛行軌跡和姿態(tài)受氣動載荷的作用時刻發(fā)生變化,有可能與導(dǎo)彈發(fā)生碰撞造成發(fā)射事故。因此,在發(fā)射之前有必要對適配器的初始分離彈道進行計算分析。
傳統(tǒng)的非耦合方法計算數(shù)個攻角下適配器的氣動系數(shù),得到一系列離散的氣動載荷數(shù)據(jù),擬合出氣動載荷隨姿態(tài)的變化曲線,然后將曲線作為外力輸入,通過求解剛體的六自由度運動微分方程計算適配器的飛行姿態(tài)與軌跡。以上方法雖具有可操作性,但是為了得到精確的計算結(jié)果,需要計算大量的工況以得到更多的數(shù)據(jù)點,耗費時間和精力。并且,在真實的環(huán)境中適配器周圍的流場會隨著適配器與導(dǎo)彈間的距離發(fā)生變化,這種變化過程是傳統(tǒng)的非耦合計算方法不能實現(xiàn)的。
本文使用基于重構(gòu)方法的六自由度動網(wǎng)格技術(shù),將適配器的運動狀態(tài)與流場的變化相耦合。每一時刻適配器受到的氣動載荷由CFD程序直接求解,然后由剛體的六自由度運動微分方程計算出適配器的位移與角度變化。適配器的運動狀態(tài)更新后成為新的初始條件影響下一時刻氣動載荷的計算。如此循環(huán)迭代可以真正實時求解適配器在初始分離彈道階段的飛行姿態(tài)與軌跡。使用耦合方法計算剛體在氣動力驅(qū)動作用下的運動過程,已在國內(nèi)外關(guān)于飛機外掛物的分離彈道研究中得到了廣泛應(yīng)用[2-4],其正確性得到了驗證。
1適配器模型
適配器安裝在導(dǎo)彈中部靠近重心的位置處,總共4個,并以導(dǎo)彈中心線為軸線按90°陣列。發(fā)射筒與地面垂直,導(dǎo)彈被彈射出筒,整個發(fā)射系統(tǒng)模型如圖1所示。
圖1 發(fā)射系統(tǒng)示意圖
適配器整體呈長方形,靠近彈頭一端有向內(nèi)的楔形切面。內(nèi)外表面均呈圓弧形,以便與導(dǎo)彈和發(fā)射筒壁面貼合。適配器中心位置有放置彈簧結(jié)構(gòu)的盲孔,其外形如圖2所示。
圖2 導(dǎo)彈適配器幾何外形
發(fā)射前適配器被發(fā)射筒和彈體約束,彈簧處于壓縮狀態(tài)。導(dǎo)彈彈射出筒的過程中適配器與導(dǎo)彈一起運動,當(dāng)適配器離開發(fā)射筒,失去了筒壁的約束作用,彈簧開始伸展使適配器產(chǎn)生沿導(dǎo)彈徑向的速度。適配器的彈出距離x與彈簧工作時間t的關(guān)系為
(1)
式中:l為彈簧的壓縮量,k為彈簧的剛度,m為適配器的質(zhì)量。由此計算出彈簧完全伸展開需要的時間為10 ms,在此時間內(nèi)適配器的速度很小,因此氣動載荷對適配器產(chǎn)生的沖量很小,可以忽略不計。因為適配器彈出的時間極短,所以在此時間內(nèi)導(dǎo)彈與適配器沿軸線方向的速度也變化很小,可以認(rèn)為兩者的速度保持相同。以彈簧完全展開時為計算的初始時刻,此時導(dǎo)彈與適配器沿軸線方向的速度為20 m/s,適配器沿徑向的速度為5.8 m/s。將適配器下落到導(dǎo)彈尾部高度時指定為初始分離彈道的結(jié)束時刻。
計算適配器的平移運動選擇全局坐標(biāo)系較為方便,計算適配器的旋轉(zhuǎn)運動選擇隨體坐標(biāo)系較為方便。全局坐標(biāo)系以大地為基準(zhǔn),垂直于地面向上為X軸正方向,垂直于X軸且在流場模型的對稱平面內(nèi)為Y軸,Z軸由右手螺旋定則得到。適配器的隨體坐標(biāo)系以適配器中心孔為基準(zhǔn);經(jīng)過重心垂直于中心孔旋轉(zhuǎn)軸且在適配器對稱平面內(nèi)為X′軸,指向適配器楔形面為正方向;垂直于X′軸經(jīng)過重心與中心孔旋轉(zhuǎn)軸平行為Y′軸,指向中心孔開口側(cè)為正方向;Z′軸由右手螺旋定則確定。
設(shè)隨體坐標(biāo)系由全局坐標(biāo)系先繞Z軸旋轉(zhuǎn)φ角,再繞Y軸旋轉(zhuǎn)θ角,然后繞X軸旋轉(zhuǎn)ψ角得到。則從全局坐標(biāo)系到隨體坐標(biāo)系的轉(zhuǎn)換矩陣R為
(2)
r11=cosθcosψ, r12=cosθsinψ, r13=-sinθ
r21=sinφsinθcosψ-cosφsinψ
r22=sinφsinθcosψ+cosφcosψ, r23=sinφcosθ
r31=cosφsinθcosψ+sinφsinψ
r32=cosψsinθsinψ-sinφcosψ, r33=cosφcosθ
使用Gambit建模工具以四面體非結(jié)構(gòu)網(wǎng)格對整個流場計算域進行網(wǎng)格劃分,網(wǎng)格單元數(shù)量為200萬。為保證適配器氣動載荷計算的精度,對適配器周圍區(qū)域進行局部加密。并且在計算過程加密區(qū)域隨適配器一起運動,以確保適配器壁面附近的網(wǎng)格密度如圖3所示。
圖3 適配器與導(dǎo)彈周圍網(wǎng)格
2理論模型
2.1運動微分方程
適配器在空中飛行的過程中,除受到的氣動力和重力外沒有其他約束條件。因此,具有6個方向上的自由度。將適配器的運動分解為兩部分,即質(zhì)心的質(zhì)點運動以及剛體的轉(zhuǎn)動。
在全局坐標(biāo)系下建立質(zhì)心的質(zhì)點運動微分方程:
(3)
式中:FX、FY、FZ分別為X、Y、Z方向上適配器受到的作用力,vX、vY、vZ分別為這3個方向上的速度分量。在適配器的隨體坐標(biāo)系下建立剛體轉(zhuǎn)動的運動微分方程:
L·ω=∑M-ω×(L·ω)
(4)
式中:L為適配器的慣性張量,ω為角速度向量,M為作用在適配器上的力矩。Admas-Moulton多步積分法非常適合求解運動微分方程,因此使用Admas-Moulton隱式三步四階格式對以上微分方程進行數(shù)值求解。計算過程中式(3)與式(4)中的力與力矩通過積分適配器表面的壓強與切應(yīng)力得到,而壓強與切應(yīng)力由CFD程序解算流場控制方程得到。
2.2網(wǎng)格重構(gòu)方法
目前常用的動網(wǎng)格方法有:光滑方法、分層方法和重構(gòu)方法。光滑方法適用于網(wǎng)格變形較小的情況,分層方法適用于網(wǎng)格運動方向恒定的情況。而適配器在空中飛行的過程中,速度方向時刻發(fā)生變化,且有滾轉(zhuǎn)、俯仰和偏航運動,周圍網(wǎng)格的變形量較大,因此選用重構(gòu)動網(wǎng)格方法最為合適。
基于恒定尺寸閾值的局部網(wǎng)格重構(gòu)方法靈活性較差,有可能造成剛體運動方向前方的網(wǎng)格尺寸過小,而運動方向后方的網(wǎng)格尺寸過大,降低計算精度甚至引起迭代發(fā)散。并且隨著剛體的運動可能導(dǎo)致計算域內(nèi)網(wǎng)格單元不斷增多,使計算速度逐漸降低。而采用基于尺寸函數(shù)的網(wǎng)格重構(gòu)方法可以有效避免以上問題。其基本原理是,根據(jù)網(wǎng)格單元中心點與邊界的距離以及邊界網(wǎng)格單元的尺寸,計算此單元網(wǎng)格重構(gòu)時的尺寸基準(zhǔn)。重構(gòu)的網(wǎng)格尺寸為邊界網(wǎng)格單元尺寸的加權(quán)平均,其權(quán)值是重構(gòu)網(wǎng)格單元與邊界網(wǎng)格單元距離的反比,即:
sp=sb×γ
(5)
式中:
(6)
(7)
式中:FI與邊界上的單元格尺寸相關(guān);LI為網(wǎng)格單元中心點距離節(jié)點的距離;db是一個無量綱量,表征網(wǎng)格單元與邊界單元之間的距離;α為尺寸函數(shù)變化率,控制網(wǎng)格單元尺寸相對于邊界網(wǎng)格單元的大小;β為尺寸函數(shù)比例,控制網(wǎng)格單元尺寸隨db的變化程度。當(dāng)運動邊界靠近時,db減小,sp也減小,此時重構(gòu)單元網(wǎng)格加密,當(dāng)運動邊界遠(yuǎn)離時,db增大,sp也增大,此時重構(gòu)單元網(wǎng)格稀疏。
2.3流場控制方程
在基于歐拉描述的流場中,任意空間位置處的標(biāo)量φ在當(dāng)?shù)氐臅r間變化率是由:流入流出此控制體的對流流量,溫度梯度、濃度梯度等帶來的擴散流量以及其他現(xiàn)象產(chǎn)生的源項引起的。因此標(biāo)量φ的輸運方程可表示為
(8)
p=ρRT
(9)
h=∫cp(T)dT
(10)
式中:ρ為當(dāng)?shù)亓鲌龅拿芏?u為流場的速度矢量,為拉普拉斯算子,Γ為擴散系數(shù),Sφ為標(biāo)量φ的產(chǎn)生源項。當(dāng)φ=1時,可以得到流場的質(zhì)量守恒方程;當(dāng)φ分別取X、Y、Z方向上的速度標(biāo)量u、v、w時,可以得到動量守恒方程;當(dāng)φ取比焓h時,可以得到能量守恒方程。以上5個方程再加上描述氣體狀態(tài)的理想氣體狀態(tài)的式(9)以及描述比焓與溫度之間關(guān)系的式(10),就可以構(gòu)成封閉的控制方程組,對基本標(biāo)量ρ、u、v、w、p、T與h進行求解。
2.4Spalart-Allmaras湍流模型
在適配器的氣動流場中,適配器的楔形面與立面的交接過渡處存在曲率較大的流線以及較強的壓強梯度,當(dāng)適配器的攻角較大時還會存在貼體流動分離現(xiàn)象。對于這樣的氣動湍流問題,Spalart-Allmaras模型經(jīng)過一定的優(yōu)化,是專門為氣動流場設(shè)計的湍流模型。且在本文研究的適配器初始分離彈道問題中導(dǎo)彈發(fā)動機沒有點火,因此不存在燃?xì)馍淞鬟@樣的自由剪切流,選用此模型作為湍流粘性模型對于本文研究的問題十分合適。Spalart-Allmaras模型是一方程湍流粘性模型,通過輸運方程求解湍流粘性νT,其表達式為
(11)
式中:σν為常數(shù),Sν為湍流粘性源項。
3計算結(jié)果
初始時刻流場內(nèi)為標(biāo)準(zhǔn)大氣環(huán)境,壓強為101.325kPa,溫度為293.15K。流場中有沿Y軸正方向速度為3m/s的橫風(fēng)。按照初始時刻全局坐標(biāo)系下適配器質(zhì)心位置給適配器命名,質(zhì)心在Y軸正方向的為“SPQ-1”,質(zhì)心在Y軸負(fù)方向的為“SPQ-2”,質(zhì)心在Z軸正方向的為“SPQ-3”,質(zhì)心在Z軸負(fù)方向的為“SPQ-4”。
3.1適配器氣動載荷變化
初始時刻適配器“SPQ-1”、“SPQ-2”迎風(fēng)面的壓力分布云圖如圖4所示。
圖4 初始時刻適配器迎風(fēng)面壓力分布
由圖4可以看出,“SPQ-2”迎風(fēng)面平均壓強要高于“SPQ-1”。選取適配器楔形面在隨體坐標(biāo)系O′Y′Z′平面上的投影面積S為參考面積,選取適配器最長邊的一半為參考長度l,由公式(12)可得適配器無量綱參數(shù)氣動載荷系數(shù)。
(12)
式中:CF為作用在適配器上的氣動力系數(shù),CM為作用在適配器上的氣動力矩系數(shù),q∞為無窮遠(yuǎn)處自由流場的動壓。適配器氣動力系數(shù)CFX、CFY、CFZ在初始分離彈道階段隨時間變化曲線如圖5所示。
圖5 氣動力系數(shù)曲線
適配器有沿X軸正方向的速度,因此受到沿X軸負(fù)方向的氣動力作用??梢钥闯?個適配器X方向的氣動力變化趨勢一致。從初始時刻到0.15s這段時間內(nèi)適配器氣動力系數(shù)CFX幅值不斷增大,因為受圖5(b)、圖5(c)所示的氣動力矩作用,在此時間內(nèi)適配器繞O′Z′轉(zhuǎn)動,隨著旋轉(zhuǎn)角度的增大,適配器在OYZ平面內(nèi)的投影面積不斷增大,使X方向上的氣動力作用面積增大,使FX幅值上升。0.15s左右4個適配器繞O′Z′旋轉(zhuǎn)角度達到90°,適配器處于水平狀態(tài)如圖11(c)所示,之后適配器繼續(xù)旋轉(zhuǎn),OYZ平面內(nèi)的投影面積不斷減小,同時受FX作用vX不斷減小,2個因素共同作用下FX幅值下降如圖5(a)所示。在0.3s左右時適配器的旋轉(zhuǎn)角度達到180°,如圖11(d)所示,此時OYZ平面內(nèi)的投影面積達到最小,且vX小于初始時刻,因此FX幅值最小。在Y、Z方向上氣動力的變化過程與此類似。
不同時刻適配器表面的壓力分布如圖6所示。初始時刻,適配器處于豎直狀態(tài),此時X方向氣動力主要作用在適配器前端的楔形面上,如圖6(a)所示。楔形面外法線與O′X′軸有夾角,因此楔形面受到的氣動力有沿O′Y′軸正方向分量,在適配器質(zhì)心處產(chǎn)生的力矩沿O′Z′軸正向,這是致使適配器旋轉(zhuǎn)的主要因素。隨體坐標(biāo)系下適配器受到的氣動力矩系數(shù)CMZ′隨時間變化曲線如圖7所示。
圖6 適配器表面壓力云圖
圖7 氣動載荷系數(shù)CMZ′曲線
可以看出,在最初的一段時間內(nèi)適配器受到的氣動力矩有所上升,這是因為適配器繞O′Z′軸旋轉(zhuǎn)后增大了X方向的氣動力作用面積,適配器迎風(fēng)面位于質(zhì)心之前的區(qū)域提供了更多的沿O′Z′軸正向的氣動力矩,但是位于質(zhì)心之后的區(qū)域產(chǎn)生的是沿O′Z′軸負(fù)向的氣動力矩,兩者疊加氣動力矩幅值有較小的上升。當(dāng)適配器轉(zhuǎn)過一定的角度時,質(zhì)心后部區(qū)域產(chǎn)生的負(fù)向氣動力矩超過前部正向氣動力矩。并且隨著適配器轉(zhuǎn)角的增大,楔形面外法線與來流之間的夾角不斷增大,由此產(chǎn)生的正向氣動力矩逐漸減小,在2個因素的共同作用下致使適配器沿O′Z′軸正向的氣動力矩隨轉(zhuǎn)角的增大而減小。在0.15s作用時適配器旋轉(zhuǎn)至水平位置,此時CMZ′達到最小值,之后適配器繼續(xù)旋轉(zhuǎn)適配器背面和底部成為迎風(fēng)面,從而產(chǎn)生了正向氣動力矩使CMZ′有所上升,但是此時適配器的運動速度已經(jīng)減小很多,因此上升幅度不大。
3.2適配器分離彈道軌跡
由圖5、圖7可以看出,適配器“SPQ-3”與“SPQ-4”的氣動載荷系數(shù)體現(xiàn)出較好的對稱性與一致性,而適配器“SPQ-1”與“SPQ-2”的氣動載荷系數(shù)有較大的差別。這是因為適配器“SPQ-1”的分離速度沿Y軸正方向與風(fēng)的方向相同,而適配器“SPQ-2”的分離速度沿Y軸負(fù)方向與風(fēng)的方向相反,因此兩者相對于空氣的運動速度不同,造成兩者的氣動載荷也不相同。而氣動載荷的非對稱性又加劇了兩者Y方向上速度的差異。適配器質(zhì)心運動速度隨時間變化曲線如圖8所示。
可以看出,4個適配器在X方向上速度變化趨勢保持一致。適配器“SPQ-2”直接受到從Y軸負(fù)方向吹來的風(fēng)影響,vY的減小速度快于適配器“SPQ-1”。同時受風(fēng)力作用,本來在Y方向沒有速度的適配器“SPQ-3”與“SPQ-4”也逐漸產(chǎn)生了速度。在Z方向上沒有風(fēng)力作用,因此適配器“SPQ-3”與“SPQ-4”的速度vZ變化趨勢完全一致,且適配器“SPQ-1”與“SPQ-2”在Z方向上的速度始終為0。適配器質(zhì)心運動軌跡呈拋物線形如圖9所示。
圖8 適配器質(zhì)心速度曲線
圖9 適配器質(zhì)心軌跡
適配器與導(dǎo)彈分離0.575s后導(dǎo)彈底部在X方向上超過所有適配器,初始分離彈道結(jié)束。設(shè)此時適配器質(zhì)心與導(dǎo)彈軸線的垂直距離為s,結(jié)束時所有適配器的速度與位移見表1。
表1 初始分離彈道終點適配器運動參數(shù)
分離過程中繞隨體坐標(biāo)系O′Z′軸正方向的旋轉(zhuǎn)是主要的旋轉(zhuǎn)運動,致使適配器前端楔形面向彈體方向翻轉(zhuǎn),旋轉(zhuǎn)角度ψ與角速度ωZ隨時間的變化曲線如圖10所示。
初始分離彈道結(jié)束時所有適配器的旋轉(zhuǎn)角速度與角度見表1。數(shù)值計算得到的適配器運動姿態(tài)與實驗中高速攝影記錄的不同時刻適配器姿態(tài)變化對比如圖11所示。
由數(shù)值計算結(jié)果與試驗結(jié)果均可以看出,初始分離彈道結(jié)束時適配器未與導(dǎo)彈發(fā)生碰撞,適配器的分離是安全的。此外,數(shù)值計算預(yù)測的適配器飛行姿態(tài)和軌跡與高速攝影記錄的實驗結(jié)果保持一致。
圖10 適配器旋轉(zhuǎn)運動變化曲線
圖11 適配器飛行姿態(tài)
5結(jié)束語
本文建立了適配器與導(dǎo)彈的初始分離彈道模型,劃分了非結(jié)構(gòu)網(wǎng)格,使用基于六自由度動網(wǎng)格技術(shù)將適配器的分離運動與氣動流場相耦合,得到了初始分離彈道階段適配器的運動軌跡與姿態(tài),并與實驗結(jié)果進行了對比,可以得出以下結(jié)論:
①計算結(jié)果與實驗結(jié)果高度一致,表明此種研究方法具有較高的可信度和準(zhǔn)確性;
②在本文研究的工況中, 適配器在初始分離階
段沒有與導(dǎo)彈發(fā)生碰撞,分離過程是安全的;
③受氣動載荷的作用適配器產(chǎn)生旋轉(zhuǎn)運動,且適配器前端楔形面向彈體方向翻轉(zhuǎn);
④沿著風(fēng)速方向?qū)ΨQ分布的2個適配器,受到的氣動載荷不同,兩者的運動狀態(tài)也不同。
本文使用的研究方法對于研究適配器與導(dǎo)彈的分離過程,預(yù)防適配器與導(dǎo)彈的碰撞,以及地面設(shè)備和人員的防護有指導(dǎo)作用。
參考文獻
趙華,王敏杰,楊為,等.箱式發(fā)射導(dǎo)彈適配器.戰(zhàn)術(shù)導(dǎo)彈技術(shù),2007(4):42-50.
ZHAO Hua,WANG Min-jie,YANG Wei,et al.Adapters for canister-launched missile.Tactical Missile Technology,2007(4):42-50.(in Chinese)
MAHMOOD T,AIZUD N,ZAHIR S.CFD simulations of the store separation and its effect on the attitude of the parent vehicle//International Bhurban Conference on Applied Sciences.Piscataway,NJ,USA:IEEE,2012:198-202.
PANAGIOTOPOULOS E E,KYPARISSIS S D.CFD transonic store separation trajectory predictions with comparison to wind tunnel investigations.International Journal of Engineering,2010(3):538-553.
張培紅,王明,鄧有奇,等.武器分離及艙門開啟過程數(shù)值模擬研究.空氣動力學(xué)學(xué)報,2013,31(3):277-293.
ZHANG Pei-hong,WANG Ming,DENG You-qi,et al.Numerical simulation of store separation and door operation.Acta Aerodynamic Sincia,2013,31(3):277-293.(in Chinese)
BURDEN R L,FAIRES J D.Numerical analysis.Youngstown:Youngstown State University,2010.
POPE S B.Turbulent flows.UK:Cambridge University Press,2000.
ANDERSON J D.Fundamentals of aerodynamics.5th ed.New York:McGraw-Hill,2009.
Numerical Analysis of Initial Separation Trajectory of
Missile Adapter Based on Aerodynamic Coupling
NIU Yu-sen1,JIANG Yi1,LI Jing1,CHEN Miao1,SHI Shao-yan1,YU Shao-zhen2
(1.School of Aerospace Engineering,BIT,Beijing 100081,China;
2.Naval Acadme of Armament,Beijing 100161,China)
Abstract:In order to accurately predict the initial separation trajectory of the adapter after releasing from the missile,the motion of the adapter was coupled with the variation of the flow field by dynamic mesh technology based on local reconstruction.The aerodynamic loads of the adapter were computed in real time by CFD routine.The 6DOF motion differential equations of the adapter were solved.The models of missile and four adapters were built,and the computational domain was meshed by using tetrahedron cells.The grid cells around the adapter were updated by using local remeshing method based on size function.The initial separation trajectory of the adapter was calculated.The results accord with the experimental result.The curves of motion velocities and aerodynamic loads coefficients were drawn.Results show that the distance between adapter and missile increases because of the initial separation velocity,and the adapters rotate because of the aerodynamic force.The adapter doesn’t collide with missile untile the end of initial separation trajectory.
Key words:adapter;CFD;six degree of freedom;dynamic mesh;aerodynamic coupling.
中圖分類號:TJ303.4
文獻標(biāo)識碼:A
文章編號:1004-499X(2015)04-0052-07
作者簡介:牛鈺森(1987- ),男,博士研究生,研究方向為兵器發(fā)射理論與技術(shù)。E-mail:shenzhou1987@aliyun.com。
收稿日期:2015-06-11