基于氣動耦合的導彈適配器初始分離彈道數值分析

牛鈺森,姜　毅,李　靜,陳　苗,史少巖,于邵禎

(1.北京理工大學 宇航學院,北京 100081;2.海軍裝備研究院,北京 100161)

?

基于氣動耦合的導彈適配器初始分離彈道數值分析

牛鈺森1,姜毅1,李靜1,陳苗1,史少巖1,于邵禎2

(1.北京理工大學 宇航學院,北京 100081;2.海軍裝備研究院,北京 100161)

摘要:為了準確預測與導彈分離后導彈適配器的初始分離彈道,使用基于局部重構方法的動網格技術,將適配器的運動與流場的變化耦合。通過CFD方法實時求解適配器飛行過程中的氣動載荷,并對適配器的六自由度運動微分方程進行求解。使用基于尺寸函數的局部重構方法對適配器周圍網格進行更新。通過數值計算得到了適配器在初始分離彈道階段的氣動載荷系數曲線,以及適配器的分離速度、角速度曲線和分離飛行軌跡,并與試驗結果進行了對比,兩者高度一致。結果表明,在分離初速的作用下，適配器與導彈的距離逐漸增大，在氣動力作用下適配器發生偏轉，直到初始分離彈道結束適配器未與導彈發生碰撞。

關鍵詞:適配器;CFD;六自由度;動網格;氣動耦合

安裝在導彈與發射箱之間的適配器能夠為導彈提供支撐,并且在發射過程中起到減震與導向的作用。對于安裝在導彈外表面的實心適配器,在導彈發射過程中隨著導彈一起運動,出箱后通過爆炸螺栓或者彈簧機構與導彈解鎖分離。在初始分離彈道階段,適配器的飛行軌跡和姿態受氣動載荷的作用時刻發生變化,有可能與導彈發生碰撞造成發射事故。因此,在發射之前有必要對適配器的初始分離彈道進行計算分析。

傳統的非耦合方法計算數個攻角下適配器的氣動系數,得到一系列離散的氣動載荷數據,擬合出氣動載荷隨姿態的變化曲線,然后將曲線作為外力輸入,通過求解剛體的六自由度運動微分方程計算適配器的飛行姿態與軌跡。以上方法雖具有可操作性,但是為了得到精確的計算結果,需要計算大量的工況以得到更多的數據點,耗費時間和精力。并且,在真實的環境中適配器周圍的流場會隨著適配器與導彈間的距離發生變化,這種變化過程是傳統的非耦合計算方法不能實現的。

本文使用基于重構方法的六自由度動網格技術,將適配器的運動狀態與流場的變化相耦合。每一時刻適配器受到的氣動載荷由CFD程序直接求解,然后由剛體的六自由度運動微分方程計算出適配器的位移與角度變化。適配器的運動狀態更新后成為新的初始條件影響下一時刻氣動載荷的計算。如此循環迭代可以真正實時求解適配器在初始分離彈道階段的飛行姿態與軌跡。使用耦合方法計算剛體在氣動力驅動作用下的運動過程,已在國內外關于飛機外掛物的分離彈道研究中得到了廣泛應用[2-4],其正確性得到了驗證。

1適配器模型

適配器安裝在導彈中部靠近重心的位置處,總共4個,并以導彈中心線為軸線按90°陣列。發射筒與地面垂直,導彈被彈射出筒,整個發射系統模型如圖1所示。

圖1　發射系統示意圖

適配器整體呈長方形,靠近彈頭一端有向內的楔形切面。內外表面均呈圓弧形,以便與導彈和發射筒壁面貼合。適配器中心位置有放置彈簧結構的盲孔,其外形如圖2所示。

圖2　導彈適配器幾何外形

發射前適配器被發射筒和彈體約束,彈簧處于壓縮狀態。導彈彈射出筒的過程中適配器與導彈一起運動,當適配器離開發射筒,失去了筒壁的約束作用,彈簧開始伸展使適配器產生沿導彈徑向的速度。適配器的彈出距離x與彈簧工作時間t的關系為

(1)

式中:l為彈簧的壓縮量,k為彈簧的剛度,m為適配器的質量。由此計算出彈簧完全伸展開需要的時間為10 ms,在此時間內適配器的速度很小,因此氣動載荷對適配器產生的沖量很小,可以忽略不計。因為適配器彈出的時間極短,所以在此時間內導彈與適配器沿軸線方向的速度也變化很小,可以認為兩者的速度保持相同。以彈簧完全展開時為計算的初始時刻,此時導彈與適配器沿軸線方向的速度為20 m/s,適配器沿徑向的速度為5.8 m/s。將適配器下落到導彈尾部高度時指定為初始分離彈道的結束時刻。

計算適配器的平移運動選擇全局坐標系較為方便,計算適配器的旋轉運動選擇隨體坐標系較為方便。全局坐標系以大地為基準,垂直于地面向上為X軸正方向,垂直于X軸且在流場模型的對稱平面內為Y軸,Z軸由右手螺旋定則得到。適配器的隨體坐標系以適配器中心孔為基準;經過重心垂直于中心孔旋轉軸且在適配器對稱平面內為X′軸,指向適配器楔形面為正方向;垂直于X′軸經過重心與中心孔旋轉軸平行為Y′軸,指向中心孔開口側為正方向;Z′軸由右手螺旋定則確定。

設隨體坐標系由全局坐標系先繞Z軸旋轉φ角,再繞Y軸旋轉θ角,然后繞X軸旋轉ψ角得到。則從全局坐標系到隨體坐標系的轉換矩陣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建模工具以四面體非結構網格對整個流場計算域進行網格劃分,網格單元數量為200萬。為保證適配器氣動載荷計算的精度,對適配器周圍區域進行局部加密。并且在計算過程加密區域隨適配器一起運動,以確保適配器壁面附近的網格密度如圖3所示。

圖3　適配器與導彈周圍網格

2理論模型

2.1運動微分方程

適配器在空中飛行的過程中,除受到的氣動力和重力外沒有其他約束條件。因此,具有6個方向上的自由度。將適配器的運動分解為兩部分,即質心的質點運動以及剛體的轉動。

在全局坐標系下建立質心的質點運動微分方程:

(3)

式中:FX、FY、FZ分別為X、Y、Z方向上適配器受到的作用力,vX、vY、vZ分別為這3個方向上的速度分量。在適配器的隨體坐標系下建立剛體轉動的運動微分方程:

L·ω=∑M-ω×(L·ω)

(4)

式中:L為適配器的慣性張量,ω為角速度向量,M為作用在適配器上的力矩。Admas-Moulton多步積分法非常適合求解運動微分方程,因此使用Admas-Moulton隱式三步四階格式對以上微分方程進行數值求解。計算過程中式(3)與式(4)中的力與力矩通過積分適配器表面的壓強與切應力得到,而壓強與切應力由CFD程序解算流場控制方程得到。

2.2網格重構方法

目前常用的動網格方法有:光滑方法、分層方法和重構方法。光滑方法適用于網格變形較小的情況,分層方法適用于網格運動方向恒定的情況。而適配器在空中飛行的過程中,速度方向時刻發生變化,且有滾轉、俯仰和偏航運動,周圍網格的變形量較大,因此選用重構動網格方法最為合適。

基于恒定尺寸閾值的局部網格重構方法靈活性較差,有可能造成剛體運動方向前方的網格尺寸過小,而運動方向后方的網格尺寸過大,降低計算精度甚至引起迭代發散。并且隨著剛體的運動可能導致計算域內網格單元不斷增多,使計算速度逐漸降低。而采用基于尺寸函數的網格重構方法可以有效避免以上問題。其基本原理是,根據網格單元中心點與邊界的距離以及邊界網格單元的尺寸,計算此單元網格重構時的尺寸基準。重構的網格尺寸為邊界網格單元尺寸的加權平均,其權值是重構網格單元與邊界網格單元距離的反比,即:

sp=sb×γ

(5)

式中:

(6)

(7)

式中:FI與邊界上的單元格尺寸相關;LI為網格單元中心點距離節點的距離;db是一個無量綱量,表征網格單元與邊界單元之間的距離;α為尺寸函數變化率,控制網格單元尺寸相對于邊界網格單元的大小;β為尺寸函數比例,控制網格單元尺寸隨db的變化程度。當運動邊界靠近時,db減小,sp也減小,此時重構單元網格加密,當運動邊界遠離時,db增大,sp也增大,此時重構單元網格稀疏。

2.3流場控制方程

在基于歐拉描述的流場中,任意空間位置處的標量φ在當地的時間變化率是由:流入流出此控制體的對流流量,溫度梯度、濃度梯度等帶來的擴散流量以及其他現象產生的源項引起的。因此標量φ的輸運方程可表示為

(8)

p=ρRT

(9)

h=∫cp(T)dT

(10)

式中:ρ為當地流場的密度,u為流場的速度矢量,為拉普拉斯算子,Γ為擴散系數,Sφ為標量φ的產生源項。當φ=1時,可以得到流場的質量守恒方程;當φ分別取X、Y、Z方向上的速度標量u、v、w時,可以得到動量守恒方程;當φ取比焓h時,可以得到能量守恒方程。以上5個方程再加上描述氣體狀態的理想氣體狀態的式(9)以及描述比焓與溫度之間關系的式(10),就可以構成封閉的控制方程組,對基本標量ρ、u、v、w、p、T與h進行求解。

2.4Spalart-Allmaras湍流模型

在適配器的氣動流場中,適配器的楔形面與立面的交接過渡處存在曲率較大的流線以及較強的壓強梯度,當適配器的攻角較大時還會存在貼體流動分離現象。對于這樣的氣動湍流問題,Spalart-Allmaras模型經過一定的優化,是專門為氣動流場設計的湍流模型。且在本文研究的適配器初始分離彈道問題中導彈發動機沒有點火,因此不存在燃氣射流這樣的自由剪切流,選用此模型作為湍流粘性模型對于本文研究的問題十分合適。Spalart-Allmaras模型是一方程湍流粘性模型,通過輸運方程求解湍流粘性νT,其表達式為

(11)

式中:σν為常數,Sν為湍流粘性源項。

3計算結果

初始時刻流場內為標準大氣環境,壓強為101.325kPa,溫度為293.15K。流場中有沿Y軸正方向速度為3m/s的橫風。按照初始時刻全局坐標系下適配器質心位置給適配器命名,質心在Y軸正方向的為“SPQ-1”,質心在Y軸負方向的為“SPQ-2”,質心在Z軸正方向的為“SPQ-3”,質心在Z軸負方向的為“SPQ-4”。

3.1適配器氣動載荷變化

初始時刻適配器“SPQ-1”、“SPQ-2”迎風面的壓力分布云圖如圖4所示。

圖4　初始時刻適配器迎風面壓力分布

由圖4可以看出，“SPQ-2”迎風面平均壓強要高于“SPQ-1”。選取適配器楔形面在隨體坐標系O′Y′Z′平面上的投影面積S為參考面積,選取適配器最長邊的一半為參考長度l,由公式(12)可得適配器無量綱參數氣動載荷系數。

(12)

式中:CF為作用在適配器上的氣動力系數,CM為作用在適配器上的氣動力矩系數,q∞為無窮遠處自由流場的動壓。適配器氣動力系數CFX、CFY、CFZ在初始分離彈道階段隨時間變化曲線如圖5所示。

圖5　氣動力系數曲線

適配器有沿X軸正方向的速度,因此受到沿X軸負方向的氣動力作用。可以看出4個適配器X方向的氣動力變化趨勢一致。從初始時刻到0.15s這段時間內適配器氣動力系數CFX幅值不斷增大,因為受圖5(b)、圖5(c)所示的氣動力矩作用,在此時間內適配器繞O′Z′轉動,隨著旋轉角度的增大,適配器在OYZ平面內的投影面積不斷增大,使X方向上的氣動力作用面積增大,使FX幅值上升。0.15s左右4個適配器繞O′Z′旋轉角度達到90°,適配器處于水平狀態如圖11(c)所示,之后適配器繼續旋轉,OYZ平面內的投影面積不斷減小,同時受FX作用vX不斷減小,2個因素共同作用下FX幅值下降如圖5(a)所示。在0.3s左右時適配器的旋轉角度達到180°,如圖11(d)所示,此時OYZ平面內的投影面積達到最小,且vX小于初始時刻,因此FX幅值最小。在Y、Z方向上氣動力的變化過程與此類似。

不同時刻適配器表面的壓力分布如圖6所示。初始時刻,適配器處于豎直狀態,此時X方向氣動力主要作用在適配器前端的楔形面上,如圖6(a)所示。楔形面外法線與O′X′軸有夾角,因此楔形面受到的氣動力有沿O′Y′軸正方向分量,在適配器質心處產生的力矩沿O′Z′軸正向,這是致使適配器旋轉的主要因素。隨體坐標系下適配器受到的氣動力矩系數CMZ′隨時間變化曲線如圖7所示。

圖6　適配器表面壓力云圖

圖7　氣動載荷系數CMZ′曲線

可以看出，在最初的一段時間內適配器受到的氣動力矩有所上升,這是因為適配器繞O′Z′軸旋轉后增大了X方向的氣動力作用面積,適配器迎風面位于質心之前的區域提供了更多的沿O′Z′軸正向的氣動力矩,但是位于質心之后的區域產生的是沿O′Z′軸負向的氣動力矩,兩者疊加氣動力矩幅值有較小的上升。當適配器轉過一定的角度時,質心后部區域產生的負向氣動力矩超過前部正向氣動力矩。并且隨著適配器轉角的增大，楔形面外法線與來流之間的夾角不斷增大,由此產生的正向氣動力矩逐漸減小,在2個因素的共同作用下致使適配器沿O′Z′軸正向的氣動力矩隨轉角的增大而減小。在0.15s作用時適配器旋轉至水平位置,此時CMZ′達到最小值,之后適配器繼續旋轉適配器背面和底部成為迎風面,從而產生了正向氣動力矩使CMZ′有所上升,但是此時適配器的運動速度已經減小很多,因此上升幅度不大。

3.2適配器分離彈道軌跡

由圖5、圖7可以看出，適配器“SPQ-3”與“SPQ-4”的氣動載荷系數體現出較好的對稱性與一致性,而適配器“SPQ-1”與“SPQ-2”的氣動載荷系數有較大的差別。這是因為適配器“SPQ-1”的分離速度沿Y軸正方向與風的方向相同,而適配器“SPQ-2”的分離速度沿Y軸負方向與風的方向相反,因此兩者相對于空氣的運動速度不同,造成兩者的氣動載荷也不相同。而氣動載荷的非對稱性又加劇了兩者Y方向上速度的差異。適配器質心運動速度隨時間變化曲線如圖8所示。

可以看出,4個適配器在X方向上速度變化趨勢保持一致。適配器“SPQ-2”直接受到從Y軸負方向吹來的風影響,vY的減小速度快于適配器“SPQ-1”。同時受風力作用,本來在Y方向沒有速度的適配器“SPQ-3”與“SPQ-4”也逐漸產生了速度。在Z方向上沒有風力作用,因此適配器“SPQ-3”與“SPQ-4”的速度vZ變化趨勢完全一致,且適配器“SPQ-1”與“SPQ-2”在Z方向上的速度始終為0。適配器質心運動軌跡呈拋物線形如圖9所示。

圖8　適配器質心速度曲線

圖9　適配器質心軌跡

適配器與導彈分離0.575s后導彈底部在X方向上超過所有適配器,初始分離彈道結束。設此時適配器質心與導彈軸線的垂直距離為s,結束時所有適配器的速度與位移見表1。

表1　初始分離彈道終點適配器運動參數

分離過程中繞隨體坐標系O′Z′軸正方向的旋轉是主要的旋轉運動,致使適配器前端楔形面向彈體方向翻轉,旋轉角度ψ與角速度ωZ隨時間的變化曲線如圖10所示。

初始分離彈道結束時所有適配器的旋轉角速度與角度見表1。數值計算得到的適配器運動姿態與實驗中高速攝影記錄的不同時刻適配器姿態變化對比如圖11所示。

由數值計算結果與試驗結果均可以看出,初始分離彈道結束時適配器未與導彈發生碰撞,適配器的分離是安全的。此外,數值計算預測的適配器飛行姿態和軌跡與高速攝影記錄的實驗結果保持一致。

圖10　適配器旋轉運動變化曲線

圖11　適配器飛行姿態

5結束語

本文建立了適配器與導彈的初始分離彈道模型,劃分了非結構網格,使用基于六自由度動網格技術將適配器的分離運動與氣動流場相耦合,得到了初始分離彈道階段適配器的運動軌跡與姿態,并與實驗結果進行了對比,可以得出以下結論:

①計算結果與實驗結果高度一致,表明此種研究方法具有較高的可信度和準確性;

②在本文研究的工況中, 適配器在初始分離階

段沒有與導彈發生碰撞,分離過程是安全的;

③受氣動載荷的作用適配器產生旋轉運動,且適配器前端楔形面向彈體方向翻轉;

④沿著風速方向對稱分布的2個適配器,受到的氣動載荷不同,兩者的運動狀態也不同。

本文使用的研究方法對于研究適配器與導彈的分離過程,預防適配器與導彈的碰撞,以及地面設備和人員的防護有指導作用。

參考文獻

趙華,王敏杰,楊為,等.箱式發射導彈適配器.戰術導彈技術,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.

張培紅,王明,鄧有奇,等.武器分離及艙門開啟過程數值模擬研究.空氣動力學學報,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

文獻標識碼:A

文章編號:1004-499X(2015)04-0052-07

作者簡介:牛鈺森(1987- ),男,博士研究生,研究方向為兵器發射理論與技術。E-mail:shenzhou1987@aliyun.com。

收稿日期:2015-06-11