，， ，俞建，

(上海機(jī)電工程研究所， 上海 201109)

0 引 言

單側(cè)提拉發(fā)射技術(shù)具有反應(yīng)快、發(fā)射率高、結(jié)構(gòu)緊湊、占用空間小、載彈量大等優(yōu)點(diǎn)[1]，因此正逐漸應(yīng)用到國內(nèi)外戰(zhàn)術(shù)導(dǎo)彈型號中。單側(cè)提拉發(fā)射系統(tǒng)通過在活塞筒和導(dǎo)軌上開設(shè)排氣孔來對彈射所產(chǎn)生的燃?xì)膺M(jìn)行排導(dǎo)。導(dǎo)彈單側(cè)提拉發(fā)射過程中，燃?xì)鈺ㄟ^活塞筒排氣孔、導(dǎo)軌排氣孔以及導(dǎo)軌與活塞筒之間的間隙直接流入發(fā)射箱中，省去了專門的燃?xì)馀艑?dǎo)結(jié)構(gòu)。但是，采用這種燃?xì)馀艑?dǎo)方式時，燃?xì)馔ㄟ^導(dǎo)軌上的排氣孔直接作用到所彈射導(dǎo)彈彈體上，所產(chǎn)生的俯仰力矩對彈體的出箱姿態(tài)具有較大的影響，作用在導(dǎo)彈尾部的燃?xì)饬Πl(fā)動機(jī)噴口產(chǎn)生不利影響[2]。為了評估燃?xì)馍淞鲗楏w姿態(tài)的影響，需要通過流固耦合發(fā)射動力學(xué)來對單側(cè)提拉發(fā)射系統(tǒng)的工作過程進(jìn)行研究，以便為發(fā)射系統(tǒng)優(yōu)化設(shè)計(jì)提供依據(jù)。

單側(cè)提拉發(fā)射系統(tǒng)流固耦合發(fā)射動力學(xué)是指將計(jì)算流體力學(xué)和多體動力學(xué)應(yīng)用于單側(cè)提拉發(fā)射過程的研究。單側(cè)提拉發(fā)射氣體動力學(xué)與剛?cè)狁詈习l(fā)射動力學(xué)是發(fā)射穩(wěn)定性動力學(xué)研究中不可分割的兩大組成部分，兩者是相互耦合、相互影響的。當(dāng)導(dǎo)彈前滑塊脫離導(dǎo)軌后，在燃?xì)饬ψ饔孟聫楏w的運(yùn)動姿態(tài)必然會發(fā)生改變，運(yùn)動姿態(tài)的改變又必然會影響彈體周圍的流場，進(jìn)而影響作用在彈體上的燃?xì)饬Α楏w運(yùn)動姿態(tài)與燃?xì)饬鲌龅南嗷プ饔脤⑹箚蝹?cè)提拉發(fā)射系統(tǒng)的燃?xì)夂蛯?dǎo)彈處于特定的流固耦合振動形態(tài)[3]。因此本文將對單側(cè)提拉發(fā)射系統(tǒng)流固耦合動態(tài)特性進(jìn)行研究。

目前流固耦合的研究方法主要有兩種：一種是利用計(jì)算流體力學(xué)耦合彈體六自由度方程來進(jìn)行計(jì)算，該方法無法考慮導(dǎo)彈滑塊與導(dǎo)軌間的碰撞力對彈體姿態(tài)的影響[4-5]；另一種是把計(jì)算流體力學(xué)的計(jì)算結(jié)果導(dǎo)入多體動力學(xué)模型中實(shí)現(xiàn)弱流固耦合仿真[6]，該方法無法考慮彈體姿態(tài)變化對燃?xì)饬Φ挠绊?。單?cè)提拉發(fā)射過程中，接觸力和燃?xì)饬ο鋸椃蛛x過程中彈體的出箱姿態(tài)都具有重要影響[3]，因此上述兩種方法的仿真精度都不能滿足要求。

為了更準(zhǔn)確地研究復(fù)雜接觸力和燃?xì)饬l件下單側(cè)提拉發(fā)射系統(tǒng)箱彈的動態(tài)分離規(guī)律，為外彈道設(shè)計(jì)提供依據(jù)，有必要引入新的方法來模擬強(qiáng)流固耦合行為，以提高仿真精度。發(fā)射氣體動力學(xué)直接耦合多體動力學(xué)提供了解決方案。該方法既考慮了導(dǎo)彈滑塊與導(dǎo)軌之間的接觸碰撞對彈體姿態(tài)的影響，也實(shí)現(xiàn)了燃?xì)饬εc彈體姿態(tài)的直接流固耦合仿真，極大地提高了單側(cè)提拉發(fā)射仿真的精度。流固耦合協(xié)同仿真環(huán)境的構(gòu)建基于MATLAB/Simulink軟件，發(fā)射氣體動力學(xué)流場及燃?xì)饬Φ那蠼饣贑FD軟件Fluent，單側(cè)提拉發(fā)射系統(tǒng)剛?cè)狁詈蟿恿W(xué)的求解基于多體動力學(xué)仿真軟件ADAMS。

1 單側(cè)提拉發(fā)射系統(tǒng)組成

如圖1所示，本文所研究的單側(cè)提拉發(fā)射系統(tǒng)主要由活塞筒、燃?xì)獍l(fā)生器(兼活塞)、制動錐、活塞桿、托架、鎖彈掛鉤、導(dǎo)軌、導(dǎo)彈、導(dǎo)彈滑塊以及發(fā)射箱等部件組成[1,3]。單側(cè)提拉發(fā)射過程中各部件之間的運(yùn)動關(guān)系為：隨著燃?xì)獍l(fā)生器高壓室內(nèi)的燃?xì)庠丛床粩嗟亓魅氲蛪菏?，活塞下腔壓?qiáng)逐漸增大，在完全克服重力和靜摩擦力的作用后，活塞桿提拉托架，托架通過滑塊帶動導(dǎo)彈沿導(dǎo)軌向上運(yùn)動；導(dǎo)彈前滑塊脫離導(dǎo)軌后，托架與制動錐發(fā)生撞擊，同時燃?xì)庑遁d孔開始卸載，燃?xì)饬魅雽?dǎo)軌內(nèi)腔，并通過導(dǎo)軌排氣孔進(jìn)行排氣，所排燃?xì)鈺苯幼饔玫綇楏w上，直到導(dǎo)彈完全出箱。

圖1 單側(cè)提拉發(fā)射系統(tǒng)組成Fig. 1 Components of one-side pulling launching system

2 發(fā)射氣體動力學(xué)模型

根據(jù)原始模型，在滿足計(jì)算要求的基礎(chǔ)上對發(fā)射箱、導(dǎo)軌、導(dǎo)彈和活塞筒模型進(jìn)行簡化處理。使用前處理工具，將整個模型關(guān)于導(dǎo)彈周面對稱建立幾何模型并生成網(wǎng)格，所生成的網(wǎng)格導(dǎo)入計(jì)算流體力學(xué)程序中，建立計(jì)算區(qū)域。模型中運(yùn)動邊界的移動涉及動網(wǎng)格技術(shù)，因本文所研究的單側(cè)提拉發(fā)射系統(tǒng)計(jì)算模型中的運(yùn)動邊界外形非常規(guī)整，所以可采用動態(tài)分層法來進(jìn)行動網(wǎng)格控制[7-9]。湍流模型選用Renormalization Group(RNG)格式帶增強(qiáng)壁面處理的k-ε湍流模型[9-11]。所有模塊加載到計(jì)算流體力學(xué)程序中形成的發(fā)射氣體動力學(xué)模型，如圖2所示。

圖2 發(fā)射氣體動力學(xué)模型Fig. 2 Launching gas dynamic model

單側(cè)提拉發(fā)射系統(tǒng)發(fā)射氣體動力學(xué)模型的流體控制方程主要包括連續(xù)性方程、動量守恒方程和能量守恒方程，其通用形式為

(1)

式中：t為時間；ρ為空氣密度；u為速度矢量；φ為流場通量；Γ為擴(kuò)散系數(shù)；S為源項(xiàng)。

3 多體動力學(xué)模型

在建立彈射系統(tǒng)多剛體模型的基礎(chǔ)上，結(jié)合關(guān)鍵部件柔性體模型和制動錐簡化力學(xué)模型，依據(jù)其拓?fù)浣Y(jié)構(gòu)關(guān)系在ADAMS仿真平臺中建立單側(cè)提拉發(fā)射系統(tǒng)的多體動力學(xué)模型[1,3]，如圖3所示。

圖3 單側(cè)提拉系統(tǒng)多體動力學(xué)模型Fig. 3 Muti-body dynamic model of one-side pulling launching system

模型中導(dǎo)彈滑塊與導(dǎo)軌間、托架與導(dǎo)軌間、托架與導(dǎo)彈滑塊間、鎖彈掛鉤與導(dǎo)彈滑塊間的約束關(guān)系采用ADAMS軟件中Contact約束來實(shí)現(xiàn)。通過該約束，ADAMS能夠采用非線性彈簧阻尼模型實(shí)現(xiàn)柔性體與剛體間的接觸力計(jì)算。

模型中的燃?xì)鈴椛淞腿細(xì)鈹_動力(力矩)由ADAMS中的單向力元模擬，通過定義狀態(tài)變量由ADAMS/Controls從Simulink中實(shí)時獲取。

單側(cè)提拉發(fā)射系統(tǒng)多體動力學(xué)方程為

(2)

4 流固耦合協(xié)同仿真方法

基于Simulink平臺,利用Fluent和ADAMS協(xié)同仿真方法處理單側(cè)提拉發(fā)射系統(tǒng)發(fā)射氣體動力學(xué)和發(fā)射多體動力學(xué)之間流固耦合問題的基本思路是[12-13]：利用接口程序聯(lián)接發(fā)射氣體動力學(xué)求解器和單側(cè)提拉發(fā)射系統(tǒng)多體動力學(xué)求解器，將Fluent作為接口程序中計(jì)算發(fā)射氣體動力學(xué)的計(jì)算引擎，實(shí)時提供彈體在當(dāng)前姿態(tài)下的燃?xì)廨d荷；接口程序再傳遞當(dāng)前燃?xì)廨d荷給ADAMS，ADAMS計(jì)算燃?xì)廨d荷作用下彈體的速度和姿態(tài)等數(shù)據(jù)；接口程序再將彈體速度和姿態(tài)等數(shù)據(jù)傳遞給Fluent,用于發(fā)射氣體動力學(xué)計(jì)算。發(fā)射氣體動力學(xué)和多體動力學(xué)間實(shí)時耦合計(jì)算過程如圖4所示。

Fluent與ADAMS間協(xié)同仿真過程中共享數(shù)據(jù)的定義如表1所示。彈體的運(yùn)動參數(shù)包括速度、姿態(tài)角、姿態(tài)角速度等；彈體所受燃?xì)廨d荷參數(shù)包括俯仰力(力矩)、偏航力(力矩)和滾轉(zhuǎn)力(力矩)。Fluent所需的彈體速度和姿態(tài)由ADAMS計(jì)算獲得，ADAMS所需的燃?xì)饬?力矩)由Fluent更新。布爾型數(shù)據(jù)Aflag為彈體速度和姿態(tài)是否寫入成功的標(biāo)志，F(xiàn)flag為燃?xì)饬κ欠駥懭氤晒Φ臉?biāo)志。Aflag和Fflag等于1時，表示成功寫入了所需數(shù)據(jù)。

圖5為單側(cè)提拉發(fā)射系統(tǒng)基于Simulink平臺的Fluent和ADAMS流固耦合協(xié)同仿真流程圖。協(xié)同仿真主要由3個部分組成：Fluent仿真、共享數(shù)據(jù)和ADAMS仿真。Fluent仿真部分包括 Define_execute_at_end和Define_cg_motion兩個接口代碼，其中：Define_execute_at_end用于實(shí)現(xiàn)存儲并更新導(dǎo)彈彈體燃?xì)饬?shù)據(jù)；Define_cg_motion用于實(shí)現(xiàn)等待彈體運(yùn)動數(shù)據(jù)更新和讀取彈體運(yùn)動數(shù)據(jù)，同時為Fluent傳遞動網(wǎng)格參數(shù)。協(xié)同仿真過程中通過Fflag和Aflag來判斷讀取數(shù)據(jù)或是否等待，當(dāng)讀取數(shù)據(jù)結(jié)束后更新Fflag和Aflag。

圖4 發(fā)射氣體動力學(xué)和多體動力學(xué)間實(shí)時耦合計(jì)算過程Fig. 4 Real-time coupling simulation process between launching gas dynamics and muti-body dynamics

數(shù)據(jù)域FdataAdataFflagAflag數(shù)據(jù)類型結(jié)構(gòu)體結(jié)構(gòu)體布爾型布爾型數(shù)據(jù)含義彈體氣動參數(shù)彈體運(yùn)動參數(shù)Fdata數(shù)據(jù)寫入標(biāo)志Adata數(shù)據(jù)寫入標(biāo)志數(shù)據(jù)更新源FluentADAMSFluent接口程序

圖5 Fluent和ADAMS流固耦合協(xié)同仿真流程圖Fig. 5 Co-simulation flow of fluid-structure coupling between Fluent and ADAMS

5 流固耦合仿真計(jì)算結(jié)果及其分析

單側(cè)提拉發(fā)射系統(tǒng)中導(dǎo)軌排氣孔的大小、數(shù)量、位置、形狀等都可能對燃?xì)饬鞯呐艑?dǎo)產(chǎn)生一定的影響，全面考慮所有因素將會使燃?xì)饬鲌龅挠?jì)算變得復(fù)雜而困難[2]。因此，本文在確定導(dǎo)軌排氣孔的大小、位置、形狀以及排列方式的前提下，僅就導(dǎo)軌排氣孔數(shù)量對發(fā)射箱內(nèi)燃?xì)饬鲌龊蛷楏w姿態(tài)的影響進(jìn)行研究。這里設(shè)計(jì)了5種研究工況，如表2 所示。

表2 計(jì)算工況

導(dǎo)彈在運(yùn)動過程中會受到X、Y、Z3個方向的力和力矩的影響，其坐標(biāo)系遵循右手定則，定義如圖6所示。因?yàn)檎麄€系統(tǒng)模型是關(guān)于XOY平面對稱的，所以可以認(rèn)為Z方向上的合力為零，在計(jì)算過程中只考慮X、Y方向上的作用力以及由此產(chǎn)生的在Z方向上的力矩。

圖6 坐標(biāo)系定義Fig. 6 Reference frame

燃?xì)馔ㄟ^活塞筒上的卸載孔進(jìn)入導(dǎo)軌內(nèi)腔，在向?qū)к壪露肆鲃拥倪^程中，依次通過導(dǎo)軌上的多個燃?xì)馀艢饪走M(jìn)入發(fā)射箱內(nèi)，噴射在導(dǎo)彈彈體表面形成多個高壓力區(qū)。這些高壓力區(qū)就是導(dǎo)彈彈體所受干擾力的主要來源。圖7給出了某時刻不同工況條件下彈體上的壓力分布云圖。

工況1 工況2 工況3 工況4 工況5 圖7 某時刻不同工況條件下彈體上的壓力分布云圖Fig. 7 Pressure nephograms of missile

如圖8所示，隨著導(dǎo)軌上的排氣孔數(shù)量逐漸減少，導(dǎo)彈所受的指向Z正向的力矩呈現(xiàn)先增大后減小的趨勢。這是因?yàn)榕艢饪诇p少后燃?xì)庾饔迷趯?dǎo)彈上的位置向后移動，X方向作用力的力臂增長，并且由于排氣孔數(shù)量的減少，燃?xì)鈱?dǎo)彈的沖擊作用增大，體現(xiàn)在導(dǎo)彈X方向作用力的變化曲線上。但是，燃?xì)庾饔迷趯?dǎo)彈上的受力面積逐漸減小，這會減小導(dǎo)彈X方向上的受力從而減小導(dǎo)彈質(zhì)心的力矩。

圖8 導(dǎo)彈質(zhì)心所受力矩Fig. 8 Moment of center of mass

如圖9所示，從后滑塊處所受力矩變化趨勢上可以看出，從工況1到工況5，后滑塊所受的力矩隨著導(dǎo)軌上排氣孔數(shù)量的減少迅速減小。這是因?yàn)殡S著排氣孔數(shù)量的減少，一方面作用在導(dǎo)彈上的受力面積逐漸減小，另一方面X方向作用力到后滑塊的力臂也逐漸縮短，兩方面的因素綜合到一起導(dǎo)致后滑塊所受的力矩減少。

圖9 導(dǎo)彈后滑塊所受力矩Fig. 9 Moment of back slider of missile

如圖10所示，與工況1相比，工況2的導(dǎo)彈X方向受力峰值有所增加，之后的工況下導(dǎo)彈X方向受力峰值依次下降。這是因?yàn)槎律弦粋€排氣孔之后，原來從該排氣孔排出的氣體從剩余的排氣孔排出，所以質(zhì)量流量增加，燃?xì)獾膭訅荷仙?，對?dǎo)彈的沖擊作用加大。但是，因?yàn)樽饔迷趯?dǎo)彈上的受力面積減小，所以導(dǎo)彈整體X方向受力沒有明顯的增加，并且隨著排氣孔數(shù)量的減少迅速下降。同時從導(dǎo)彈X方向的受力曲線上可以觀察到其幅值有幾次震蕩，這是由于燃?xì)庠诎l(fā)射箱內(nèi)來回反射造成的。隨著時間的推移，燃?xì)庵饾u排出發(fā)射箱外，發(fā)射箱內(nèi)的流場趨于平穩(wěn)，導(dǎo)彈X方向上的受力也趨于平穩(wěn)。

圖10 導(dǎo)彈質(zhì)心X方向受力對比Fig. 10 Force of center of mass for X direction

如圖11所示，從工況1 到工況5，導(dǎo)彈Y方向上的受力逐漸增大。這是因?yàn)樵居幸徊糠秩細(xì)鈴膶?dǎo)軌上的排氣孔排出后會沿著導(dǎo)彈與導(dǎo)軌之間的間隙向發(fā)射箱口運(yùn)動并排出箱外，并且隨著導(dǎo)軌上排氣孔的數(shù)量逐漸減少，這部分燃?xì)馑嫉谋壤絹碓缴?，更多的燃?xì)馔ㄟ^活塞筒與導(dǎo)軌后部的縫隙排出(燃?xì)獾乃俣确较蛑赶虬l(fā)射箱的底部，燃?xì)鉀_擊到發(fā)射箱的底部后會向前反射并產(chǎn)生壓縮波，燃?xì)庾矒舻綄?dǎo)彈底部時將會產(chǎn)生導(dǎo)彈Y方向上受力的峰值)。這個過程也可以從導(dǎo)彈X方向上的受力峰值和導(dǎo)彈Y方向上的受力峰值出現(xiàn)的時間差上得到印證。

圖11 導(dǎo)彈Y 方向受力對比Fig. 11 Force for Y direction

如圖12～13所示，除了工況3，其它工況下導(dǎo)彈底部中心點(diǎn)的壓強(qiáng)相差不大，同時峰值時間很短，超過3個標(biāo)準(zhǔn)大氣壓的作用時間范圍在1～3 ms 之間，這與導(dǎo)彈Y方向受力變化趨勢相符合。從溫度的變化曲線上可以看出，溫度的變化有很大隨機(jī)性，這是因?yàn)閺膶?dǎo)軌后部與活塞筒之間的間隙排出的燃?xì)猓瑫纬赏牧鳒u。當(dāng)這部分燃?xì)饩淼綄?dǎo)彈尾部時就會導(dǎo)致溫度升高，但是隨著燃?xì)庵饾u排出發(fā)射箱，彈底部的溫度整體上呈現(xiàn)下降的趨勢。

圖12 導(dǎo)彈底部中心點(diǎn)壓強(qiáng)對比Fig. 12 Pressure of center of missile bottom

圖13 導(dǎo)彈底部中心點(diǎn)溫度對比Fig. 13 Temperature of center of missile bottom

圖14～25給出了單側(cè)提拉發(fā)射系統(tǒng)流固耦合仿真計(jì)算結(jié)果中彈體位移、速度和姿態(tài)隨時間變化的情況。從彈體的位移、速度和姿態(tài)計(jì)算結(jié)果可以看出：①導(dǎo)彈前滑塊離開導(dǎo)軌前，燃?xì)膺€沒有開始卸載，各工況的位移、速度、角速度等曲線重合；②當(dāng)導(dǎo)彈前滑塊剛脫離導(dǎo)軌時，燃?xì)庑遁d孔開始卸載，隨著導(dǎo)軌排氣孔數(shù)量的減少，在同一時刻點(diǎn)，導(dǎo)彈的橫向(X方向)位移、橫向速度、俯仰角、俯仰角速度逐漸減??；③導(dǎo)軌上有5個排氣孔和有4個排氣孔條件下，箱彈分離過程中彈體各姿態(tài)參數(shù)變化較接近，說明5個排氣孔和4個排氣孔的作用效果相當(dāng)；④為了為外彈道設(shè)計(jì)提供準(zhǔn)確的姿態(tài)參數(shù)，需要合理選取導(dǎo)軌排氣孔的個數(shù)。

圖14 X方向位移對比Fig. 14 Displacement for X direction

圖15 Y方向位移對比Fig. 15 Displacement for Y direction

圖16 Z方向位移對比Fig. 16 Displacement for Z direction

圖17 X方向速度對比Fig. 17 Velocity for X direction

圖18 Y方向速度對比Fig. 18 Velocity for Y direction

圖19 Z方向速度對比Fig. 19 Velocity for Z direction

圖20 偏航角對比Fig.20 Yaw angle

圖21 俯仰角對比Fig.21 Pitching angle

圖22 滾轉(zhuǎn)角對比圖Fig. 22 Roll angle

圖23 偏航角速度對比Fig.23 Yaw angular velocity

圖24 俯仰角速度對比Fig.24 Pitching angular velocity

圖25 滾轉(zhuǎn)角速度對比Fig.25 Roll angular velocity

6 結(jié) 論

在提出單側(cè)提拉發(fā)射系統(tǒng)設(shè)計(jì)方案的基礎(chǔ)上，借助發(fā)射氣體動力學(xué)和多體動力學(xué)，基于Simulink軟件構(gòu)建了ADAMS和Fluent相結(jié)合的協(xié)同仿真環(huán)境。采用交替求解剛?cè)狁詈蟿恿W(xué)和發(fā)射氣體動力學(xué)的方式實(shí)現(xiàn)了單側(cè)提拉發(fā)射系統(tǒng)流固耦合協(xié)同仿真計(jì)算，評估了導(dǎo)彈單側(cè)提拉發(fā)射系統(tǒng)流固耦合動力學(xué)行為。主要研究結(jié)論如下：

1) 發(fā)射氣體動力學(xué)直接耦合多體動力學(xué)進(jìn)行聯(lián)合仿真是開展導(dǎo)彈發(fā)射動力學(xué)分析的一條有效途徑，研究結(jié)果對發(fā)射系統(tǒng)優(yōu)化設(shè)計(jì)具有重要的工程價值。

2) 隨著導(dǎo)軌上的排氣孔數(shù)量逐漸減少，導(dǎo)彈所受的指向Z正向的力矩呈現(xiàn)先增大后減小的趨勢；導(dǎo)彈后滑塊所受的力矩隨著導(dǎo)軌上排氣孔數(shù)量的減少迅速減小。

3) 導(dǎo)彈前滑塊離開導(dǎo)軌前，燃?xì)膺€沒有開始卸載，各工況的位移、速度、角速度等曲線重合；當(dāng)導(dǎo)彈前滑塊剛脫離導(dǎo)軌時，燃?xì)庑遁d孔開始卸載，隨著導(dǎo)軌排氣孔的減少，在同一時刻點(diǎn)，彈的橫向(X方向)位移、橫向速度、俯仰角、俯仰角速度逐漸減小。

4) 導(dǎo)軌上有5個排氣孔和有4個排氣孔條件下，箱彈分離過程中彈體各姿態(tài)參數(shù)變化較接近，說明5個排氣孔和4個排氣孔的作用效果相當(dāng)；為了為外彈道設(shè)計(jì)提供準(zhǔn)確的姿態(tài)參數(shù)，需要合理選取導(dǎo)軌排氣孔的個數(shù)。