柔性火箭彈氣動(dòng)力計(jì)算精度預(yù)測

李敏嬌， Laith K. Abbas，芮筱亭，王國平

(南京理工大學(xué) 發(fā)射動(dòng)力學(xué)研究所，南京 210094)

?

柔性火箭彈氣動(dòng)力計(jì)算精度預(yù)測

(南京理工大學(xué) 發(fā)射動(dòng)力學(xué)研究所，南京 210094)

為驗(yàn)證小變形小攻角下基于準(zhǔn)定常和線性理論進(jìn)行柔性火箭彈氣動(dòng)力計(jì)算的合理性，采用計(jì)算流體力學(xué)(CFD)方法預(yù)測了某火箭彈的升力及阻力系數(shù)，并在此基礎(chǔ)上與經(jīng)驗(yàn)公式相結(jié)合，獲得該火箭彈沿軸向的單位長度升力和阻力系數(shù)導(dǎo)數(shù)分布.應(yīng)用多體系統(tǒng)傳遞矩陣法計(jì)算了該火箭彈的固有頻率和振型，在準(zhǔn)定常及線性理論的假設(shè)條件下,結(jié)合已得到的氣動(dòng)力系數(shù)導(dǎo)數(shù),計(jì)算了基于火箭彈一階振型的兩個(gè)彎曲模型在不同攻角下的單位長度氣動(dòng)力系數(shù)分布.與CFD計(jì)算結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)該火箭彈在實(shí)際飛行條件下，該計(jì)算方法精度較好，在遠(yuǎn)超于實(shí)際的大變形時(shí)、攻角較大的情況下在頭部和彈身處仍具有一定精度，但尾部誤差較大，總體精度降低.

多體系統(tǒng)傳遞矩陣法；柔性火箭彈；振動(dòng)特性；CFD；氣動(dòng)力

柔性彈箭動(dòng)力學(xué)建模通常采用基于能量原理的拉格朗日法建立封閉的柔性彈箭氣動(dòng)彈性和飛行力學(xué)方程組，在小攻角和小變形條件下，對(duì)彈體彈性變形和大運(yùn)動(dòng)進(jìn)行耦合[1-6].其中，對(duì)彈箭振動(dòng)特性和氣動(dòng)力參數(shù)的準(zhǔn)確獲取是柔性彈道計(jì)算和穩(wěn)定性分析的關(guān)鍵.線性多體系統(tǒng)傳遞矩陣法特別適合于解決線性多體系統(tǒng)振動(dòng)特性、正交性、動(dòng)力響應(yīng)等問題[7]，兼具方便建模，計(jì)算速度快和精度高的優(yōu)點(diǎn).計(jì)算流體力學(xué)(CFD)方法可以很好的預(yù)測復(fù)雜幾何體的氣動(dòng)特性參數(shù)和流體流動(dòng)現(xiàn)象[8-11]，但是如果將其應(yīng)用到整個(gè)的封閉方程里則計(jì)算量太大，計(jì)算周期長，不利于快速計(jì)算和通過系統(tǒng)方程進(jìn)行理論分析.而基于氣動(dòng)力系數(shù)導(dǎo)數(shù)的氣動(dòng)力計(jì)算方法則在小變形和小攻角范圍內(nèi)可以兼顧精度高和計(jì)算簡單的特點(diǎn).

本文采用多體系統(tǒng)傳遞矩陣法和CFD與經(jīng)驗(yàn)公式相結(jié)合的方法，分別計(jì)算了火箭彈的振動(dòng)特性和氣動(dòng)力系數(shù)導(dǎo)數(shù)，在此基礎(chǔ)上計(jì)算彎曲火箭彈沿軸向單位長度氣動(dòng)力系數(shù)分布，并通過CFD方法對(duì)該方法的精度進(jìn)行了驗(yàn)證.

1 模型建立

1.1 火箭彈結(jié)構(gòu)建模

圖1 非均質(zhì)梁

火箭彈軸上的每一點(diǎn)所受的軸向力大小不同，等于該點(diǎn)至火箭彈頂端部分受到的慣性力.底部所受軸向力大小等于推力，頂端所受軸向力為零.Si取第i段質(zhì)心處所受的軸向力[3]為

圖2 第i段梁的狀態(tài)矢量

受均勻軸向壓力作用的均質(zhì)Euler-Bernoulli梁橫向自由振動(dòng)微分方程為[14]

推導(dǎo)可得長度為l的Euler-Bernoulli梁的傳遞矩陣為

其中：

對(duì)于非均質(zhì)梁，連接點(diǎn)處在模態(tài)坐標(biāo)下的狀態(tài)矢量為Z0,1，Z1,2，Z2,3，…，Zi-1,i，…，Zn,n+1，每一段的傳遞方程為

(1)

因此總傳遞方程為

其中，Uall為總傳遞矩陣Uall=UnUn-1…Ui…U2U1.

(2)

1.2 氣動(dòng)力模型

火箭彈彈身細(xì)長，翼面較小，并且翼面在弦長方向的尺寸占彈身總長很少，對(duì)于這樣的柔性彈，在線性理論范圍內(nèi)可以將彈體沿軸向劃分為若干微段，各段所受到的氣動(dòng)力為作用于該段壓力中心的合力.在準(zhǔn)定常理論和線性理論的假設(shè)下，不計(jì)氣流的尾流影響，各段的氣動(dòng)力只和該段的當(dāng)?shù)毓ソ怯嘘P(guān)，認(rèn)為變形后的彈體各段當(dāng)?shù)厝嵝怨ソ堑扔谠摱蔚暮铣伤俣认蛄颗c該段弦線之間的瞬時(shí)夾角，采用數(shù)值模擬方法計(jì)算出不同Ma數(shù)下的氣動(dòng)力系數(shù)導(dǎo)數(shù)的分布，可直接將氣動(dòng)力表述成當(dāng)?shù)厝嵝怨ソ堑暮瘮?shù)[5-6,15-16].

基于線性理論和準(zhǔn)定常理論，彈體所受整體的阻力、側(cè)向力和升力分別為

(3)

(4)

2 計(jì)算結(jié)果

2.1 火箭彈振動(dòng)特性

利用已推導(dǎo)的多體系統(tǒng)傳遞矩陣法公式計(jì)算某火箭彈被動(dòng)段時(shí)的頻率和振型，被動(dòng)段時(shí)質(zhì)量不變，推力為零.如圖3所示，當(dāng)火箭彈所分段數(shù)超過50時(shí)1階頻率的計(jì)算結(jié)果沒有明顯的變化，說明將火箭彈簡化為50段均質(zhì)梁即可以保證計(jì)算的精度.圖4為該火箭彈前4階振型，x/L=0為火箭彈尾部.

圖3 1階頻率隨單元數(shù)的變化

2.2 火箭彈氣動(dòng)力系數(shù)計(jì)算

由于火箭彈為軸對(duì)稱，其法向和側(cè)向具有相同的氣動(dòng)特性，因此本文中對(duì)側(cè)向氣動(dòng)力不予考慮.如圖5所示，火箭彈為剛體模型，沿著軸向?qū)⒒鸺龔椙谐?3段，每一個(gè)小段的中間點(diǎn)作為氣動(dòng)力的作用點(diǎn)，采用CFD方法計(jì)算沿軸向的氣動(dòng)力系數(shù)分布[17-18].計(jì)算域來流條件為靜壓101 325 Pa，靜溫288 K，計(jì)算域外邊界設(shè)為壓力遠(yuǎn)場，即基于Rieman不變量的無反射邊界條件，火箭彈物面設(shè)為無滑移壁面.采用SSTk-ω湍流模型，基于密度求解器進(jìn)行求解，離散格式采用AUSM格式和隱式求解格式，對(duì)流通量采用二階迎風(fēng)格式.本文采用多Block生成流場拓?fù)浣Y(jié)構(gòu)，外O Block生成彈體邊界層的方法，生成了高質(zhì)量的結(jié)構(gòu)化網(wǎng)格，在邊界層內(nèi)對(duì)網(wǎng)格進(jìn)行加密，保證y+<1，網(wǎng)格數(shù)為310萬.如圖7所示，獲得模型在α=0°和α=2°情況下的氣動(dòng)參數(shù)并與以往的實(shí)驗(yàn)數(shù)據(jù)比較，CL和CD的誤差均在10%以內(nèi)，驗(yàn)證了該數(shù)值計(jì)算方法的可行性與準(zhǔn)確性.圖6和圖8分別是Ma=2，α=2°工況下的速度云圖和火箭彈的升力、阻力系數(shù)分布.通過式(3)～(4)可得該火箭彈的氣動(dòng)力系數(shù)導(dǎo)數(shù)分布.

圖4 火箭彈前4階振型

圖5 火箭彈面網(wǎng)格

注：彩圖見電子版(http://hit.alljournals.cn)(2016年第10期)

圖6 速度云圖

注：彩圖見電子版(http://hit.alljournals.cn)(2016年第10期)

圖7 某火箭彈氣動(dòng)力系數(shù)隨Ma變化曲線及與實(shí)驗(yàn)數(shù)據(jù)對(duì)比

圖8 Ma=2, α=2°時(shí)火箭彈氣動(dòng)力系數(shù)分布

Fig.8 Rocket’s aerodynamic coefficient distribution withMa=2,α=2°

2.3 彎曲模型氣動(dòng)力計(jì)算方法驗(yàn)證

由于彈體變形基本基于1階振型，因此本文基于1階振型建立了彎曲的彈體模型.考慮到建模時(shí)可能存在尺寸誤差及CFD計(jì)算誤差，同時(shí)為了分析本文中氣動(dòng)力計(jì)算方法的適用范圍，分別選擇廣義坐標(biāo)ζ1=0.025和ζ1=0.100，對(duì)應(yīng)的模型示意，如圖9所示.圖10為兩個(gè)模型的彈體斜率，ζ1=0.025和ζ1=0.100遠(yuǎn)大于火箭彈實(shí)際正常飛行時(shí)的變形量，但是當(dāng)彈體發(fā)生顫振或者彎曲發(fā)散時(shí)變形劇烈，很有可能達(dá)到這樣的變形程度.

分別基于氣動(dòng)力系數(shù)導(dǎo)數(shù)和采用CFD方法計(jì)算了火箭彈在這兩種彎曲變形下，Ma=2、不同攻角的氣動(dòng)力系數(shù)分布，并對(duì)計(jì)算結(jié)果進(jìn)行對(duì)比.為了捕捉到微小變形對(duì)氣流產(chǎn)生的影響，第1層網(wǎng)格高度需要足夠小，此處設(shè)置y+<1，網(wǎng)格數(shù)為380萬.對(duì)于ζ1=0.025，即變形相對(duì)較小的情況下，攻角取0°、2°、4°時(shí)的單位長度阻力系數(shù)、單位長度升力系數(shù)分布的對(duì)比結(jié)果如圖11所示，可以發(fā)現(xiàn)兩種方法結(jié)果吻合較好，阻力系數(shù)偏差在5%以內(nèi)，升力系數(shù)偏差在10%以內(nèi).對(duì)于ζ1=0.100，即變形相對(duì)較大的的情況下，攻角取0°、2°、4°時(shí)的單位長度阻力系數(shù)、單位長度升力系數(shù)分布的對(duì)比結(jié)果如圖12所示，在0°攻角時(shí)，兩種方法的結(jié)果吻合較好；在2°攻角時(shí)彈的頭部和尾部等局部阻力系數(shù)有一定的誤差，但是誤差較小；在4°攻角時(shí)，頭部和尾部等局部誤差增大.ζ1=0.100的變形下在4°攻角對(duì)稱面速度云圖如圖13所示，由于變形增大，攻角增大，頭部及彈身的尾流對(duì)火箭彈尾部的流動(dòng)產(chǎn)生了復(fù)雜影響，此時(shí)，線性理論已經(jīng)不再適用.

圖9 彎曲彈體模型

圖10 彎曲彈體斜率

圖11 ζ1=0.025的彎曲火箭彈沿彈軸單位長度氣動(dòng)力系數(shù)分布

圖12 ζ1=0.100的彎曲火箭彈沿彈軸單位長度氣動(dòng)力系數(shù)分布

圖13 α=4°時(shí)ζ1=0.100的彎曲火箭彈速度云圖

注：彩圖見電子版(http://hit.alljournals.cn)(2016年第10期)

3 結(jié) 論

1) 推導(dǎo)了考慮軸向力作用的Euler-Bernoulli梁的傳遞矩陣，將某真實(shí)火箭彈簡化成自由-自由的非均質(zhì)Euler-Bernoulli梁，應(yīng)用多體系統(tǒng)傳遞矩陣法計(jì)算了其頻率和振型.

2) 采用基于SSTk-ω湍流模型的CFD方法計(jì)算了該火箭彈在不同Ma數(shù)、不同攻角下的氣動(dòng)力系數(shù)，與實(shí)驗(yàn)數(shù)據(jù)相比最大誤差不超過10%，證明了采用數(shù)值方法計(jì)算火箭彈氣動(dòng)力系數(shù)分布可以保證足夠的精度以及采用數(shù)值方法驗(yàn)證彎曲火箭彈氣動(dòng)力建模的可行性.

3) 基于該火箭彈1階振型建立了兩種不同變形程度的彎曲彈體模型，在準(zhǔn)定常理論和線性理論的假設(shè)條件下，結(jié)合氣動(dòng)力系數(shù)導(dǎo)數(shù)進(jìn)行兩個(gè)彎曲模型的氣動(dòng)力計(jì)算，得到了Ma=2，不同攻角兩個(gè)模型的單位長度升力系數(shù)分布和阻力系數(shù)分布，并采用CFD方法對(duì)計(jì)算結(jié)果進(jìn)行驗(yàn)證.通過驗(yàn)證發(fā)現(xiàn)，該火箭彈在實(shí)際飛行可能的變形范圍內(nèi)，基于準(zhǔn)定常理論和線性理論的柔性彈箭氣動(dòng)力計(jì)算方法結(jié)果精度較好，在遠(yuǎn)超于實(shí)際的大變形時(shí)，該方法在彈頭和彈身處仍具有一定精度，但尾部誤差較大，計(jì)算精度降低.

[1] 徐明友. 彈箭飛行動(dòng)力學(xué)[D]. 北京: 國防工業(yè)出版社, 2003.

XU Mingyou. Flight dynamics of missiles[D]. Beijing: National Defence of Industry Press,2003.

[2] 何斌, 芮筱亭, 陸毓琪. 柔性彈箭飛行力學(xué)建模研究[J]. 彈道學(xué)報(bào), 2006, 18(1): 22-24. DOI：10.3969/j.issn.1004-499X.2006.01.006.

HE Bin, RUI Xiaoting, LU Yuqi. A study on flight dynamic modeling of flexible shell rocket[J]. Journal of Ballistics, 2006, 18(1): 22-24. DOI：10.3969/j.issn.1004-499X.2006.01.006.

[3] POURTAKDOUST S H, ASSADIAN N. Investigation of thrust effect on the vibrational characteristics of flexible guided missiles[J]. Journal of Sound and Vibration, 2004, 272 (1/2): 287-299. DOI:10.1016/S0022-460X(03)00779-X.

[4] JOSHI A. Free vibration characteristics of variable mass rocket having large axial thrust/acceleration [J]. Journal of Sound and Vibration, 1995, 187 (4): 727-736. doi:10.1006/jsvi.1995.0559.

[5] PLATUS D H. Aeroelastic stability of slender, spinning missiles[J]. Journal of Guidance, Control and Dynamics, 1992, 15(1): 144-151. DOI: 10.2514/3.20812.

[6] 黃曉鵬, 樊紅祥, 蔣厚洸. 細(xì)長旋轉(zhuǎn)飛行器的氣動(dòng)彈性研究[J]. 北京理工大學(xué)學(xué)報(bào), 1998, 18(4): 431-435. DOI : 10. 15918/j.tbit1001-0645. 1998. 04. 008.

HUANG Xiaopeng, FAN Hongxiang, JIANG Houguang. Aeroelastic analysis of slender spinning vehicle[J]. Journal of Beijing Institute of Technology, 1998, 18(4): 431-435. DOI : 10. 15918/j.tbit1001-0645. 1998. 04. 008.

[7] 芮筱亭, 贠來峰, 陸毓琪, 等. 多體系統(tǒng)傳遞矩陣法及其應(yīng)用[M]. 北京: 科學(xué)出版社, 2008.

RUI Xiaoting, YUN Laifeng, LU Yuqi, et al. Transfer matrix method of multibody system and its applications[D]. Beijing: Science Press, 2008.

[8] MENTER F R. Two-equation eddy-viscosity models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605. DOI: 10.2514/3.12149.

[9] KLEB W L, WOOD W A, GNOFFO P A, et al. Computational aeroheating predictions for X-34[J]. Journal of Spacecraft and Rockets, 1999, 36(2): 179-188. DOI: 10.2514/2.3448.

[10]SAHU J, EDGE H L, HEAVEY K R, et al. Computational fluid dynamics modeling of multi-body missile aerodynamic interference[R]. U.S. Army Research Laboratory, ARL-TR-1765, Aberdeen Proving Ground, MD, 1998.

[11]DESPIRITO J, PLOSTINS P. CFD Prediction of M910 projectile aerodynamics: unsteady wake effect on magnus moment[C]//AIAA Atmospheric Flight Mechanics Conference and Exhibit 20-23 August 2007. Hilton Head, South Carolina: AIAA-2007-6580, 2007.

[12]臧濤成. 大長細(xì)比模擬彈模態(tài)分析與計(jì)算[J]. 彈箭與制導(dǎo)學(xué)報(bào), 2003, 23(2): 40-42. DOI：10.3969/j.issn.1673-9728.2003.02.013.

ZANG Taocheng. Great slenderness ratio imitative projectile mode analysis and theories calculation[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2003, 23(2): 40-42. DOI：10.3969/j.issn.1673-9728.2003.02.013.

[13]王良明. 柔性彈丸飛行動(dòng)力學(xué)研究[J]. 兵工學(xué)報(bào), 1998, 19(3): 208-213.

WANG Liangming. Flight dynamics of flexible projectiles[J]. Acta Armamentarii, 1998, 19(3): 208-213.

[14]ABBAS L K, LI Minjiao, RUI Xiaoting. Transfer matrix method for the determination of the natural vibration characteristics of realistic thrusting launch vehicle—Part I[J/OL]. Mathematical Problems in Engineering, 2013: 764673. http://www.hindawi.com/journals/mpe/2013/764673/ .DOI:10.1155/2013/764673.

[15]王華畢, 黃曉鵬, 吳甲生. 大長徑比旋轉(zhuǎn)火箭彈氣動(dòng)彈性數(shù)值計(jì)算[J]. 彈箭與制導(dǎo)學(xué)報(bào), 2006, 26(2): 242-245.DOI: 10.3969/j.issn.1673-9728.2006.02.081.

WANG Huabi, HUANG Xiaopeng, WU Jiasheng. The aeroelastic numerical calculation of large f ineness ratio spinning rocket[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2006, 26(2): 242-245. DOI: 10.3969/j.issn.1673-9728.2006.02.081.

[16]JEGARKANDI M F, NOBARI A S, SABZEHPRAVAR M, et al. Aeroelastic stability consideration of supersonic flight vehicle using nonlinear aerodynamic response surfaces[J]. Journal of Fluids and Structures, 2009, 25(6): 1079-1101. DOI:10.1016/j.jfluidstructs.2009.04.006.

[17]ABBAS L K, CHEN Dongyang, Rui Xiaoting. Normal force computation for axisymmetric multistage launch vehicle[J]. Applied Mechanics and Materials, 2013, 419: 23-29.DOI: 10.4028/www.scientific.net/AMM.419.23.

[18]ABBAS L K, CHEN Dongyang, Rui Xiaoting. Numerical calculation of effect of elastic deformation on aerodynamic characteristics of a rocket[J/OL]. International Journal of Aerospace Engineering, 2014: 478534. http://www.hindawi.com/journals/ijae/2014/478534. DOI:10.1155/2014/478534.

(編輯 張 紅)

The prediction of flexible rocket’s aerodynamic force computational accuracy

LI Minjiao, Laith K. Abbas, RUI Xiaoting, WANG Guoping

(Institute of Launch Dynamics, Nanjing University of Science and Technology, Nanjing 210094, China)

To verify the reasonability of calculating the flexible rocket’s aerodynamic forces based on the quasi steady theory and the linear theory when both the deformation and the attack angle are small enough, the computational fluid dynamics (CFD) method is used to calculate the rocket’s lift and drag coefficients. Using this method and the empirical formula, the rocket’s aerodynamic coefficient derivatives distributions are obtained. The transfer matrix method of the linear multibody system dynamics (MSTMM) is applied to get the rocket’s frequencies and mode shapes. Under the assumption of quasi-steady theory and linear theory, using the obtained aerodynamic coefficient derivatives distribution, the aerodynamic forces of two bended rochets’ the deformations of which are based on the first order mode shape are calculated. To judge the reasonability of the linear theory, the CFD method is also used to simulate the two models. Compared with the CFD method, the proposed method has high precision under the practical flight conditions. In addition, when the bending degree is much higher than the possible deformation, the precision of the rocket’s nose is still fine but very low in the tails position for the proposed method.

transverse matrix method of multi-body system; flexible rocket; vibration characteristic; CFD; aerodynamic force

10.11918/j.issn.0367-6234.2016.10.013

2015-03-09

國家自然科學(xué)基金(11472135, 61304137)；新世紀(jì)優(yōu)秀人才支持計(jì)劃(NCET-10-0075)；高校博士點(diǎn)基金(20113219110025, 20133219110037)

李敏嬌(1988—)，女，博士生;Laith K. Abbas(1965—)，男，教授，博士生導(dǎo)師

Laith K. Abbas，laithabbass@yahoo.com

V211

A

0367-6234(2016)10-0091-06