內(nèi)埋式彈艙與彈體相互影響的精細模擬

張群峰， 閆盼盼， 黎軍

(1.北京交通大學 土木工程學院， 北京 100044; 2.中國航空工業(yè)集團公司 沈陽飛機設(shè)計研究所， 遼寧 沈陽 110035)

內(nèi)埋式彈艙與彈體相互影響的精細模擬

張群峰1， 閆盼盼1， 黎軍2

(1.北京交通大學 土木工程學院， 北京 100044; 2.中國航空工業(yè)集團公司 沈陽飛機設(shè)計研究所， 遼寧 沈陽 110035)

為了研究內(nèi)埋式武器艙彈體投放過程中下落彈體與艙體之間強烈的相互作用，利用基于Menter SSTk-ω湍流模型的分離渦模擬方法，結(jié)合六自由度剛體動力學方程和重疊網(wǎng)格技術(shù)對某一簡化開式彈艙和彈體模型的三維流場進行了非定常計算，利用非平穩(wěn)信號處理方法平滑偽Wigner-Ville分布分析了艙內(nèi)測點壓力的時頻特性。研究結(jié)果表明：在彈體出艙過程中，剪切層被破壞，導致艙體內(nèi)流動狀態(tài)發(fā)生較大改變，不再呈現(xiàn)開式艙體自持振蕩特性，艙體內(nèi)噪聲的各階單調(diào)聲模態(tài)不明顯，受彈體的影響，艙體內(nèi)強度較高的渦結(jié)構(gòu)主要集中于艙體后緣，使艙體后緣噪聲水平升高；彈體出艙后剪切層迅速重建，自持振蕩回路再次形成，艙體流動狀態(tài)恢復為典型的純空艙流動狀態(tài)，但受彈體頭部加速氣流的影響，剪切層不穩(wěn)定性增強，導致艙體內(nèi)部噪聲水平升高。彈體在出艙過程中由于剪切層影響受到抬頭力矩的作用，使得彈體具有較大的迎角。彈體出艙后，在迎角的影響下，彈體開始受豎直向上的力和低頭力矩的作用，豎直向上的力會阻礙彈體的下落。同時，彈體下落過程受艙體強非定常流場的影響，彈體受力也存在強烈的波動。研究結(jié)果為內(nèi)埋武器艙優(yōu)化設(shè)計和制定武器安全投放控制規(guī)律提供參考。

流體力學； 內(nèi)埋艙； 武器分離； 重疊網(wǎng)格； SST分離渦模擬； 平滑偽Wigner-Ville分布

0 引言

武器裝載于內(nèi)埋艙有助于戰(zhàn)斗機降低阻力和減小雷達探測面積，其逐漸替代外掛式而被廣泛采用。然而當艙門開啟時，開式流動內(nèi)埋式武器艙內(nèi)部存在極其復雜的流動物理現(xiàn)象，彈艙內(nèi)會形成大幅度的壓力脈動，造成劇烈震蕩并產(chǎn)生刺耳的噪聲，對控制武器的電子系統(tǒng)及相關(guān)結(jié)構(gòu)造成疲勞破壞[1]。

為了研究內(nèi)埋式彈艙的復雜流場，研究人員通常把彈艙簡化為矩形空腔來研究流場的穩(wěn)態(tài)和動態(tài)特性，并已經(jīng)進行了大量的試驗及數(shù)值模擬研究，研究表明：腔體的長深比[2-3]、腔體形狀[4]、來流馬赫數(shù)[5]、來流邊界層厚度[6-7]及有無控制措施[8]等因素均會對腔體內(nèi)噪聲環(huán)境產(chǎn)生影響?；诤喕兛涨涣鲃拥难芯拷Y(jié)果，文獻[9-10]分別對艙體帶彈和彈體固定在艙體不同位置時的流場進行了研究，給出了彈體固定時艙體流場特性。值得提出的是，飛機在進行武器投放時，彈體下落過程中，彈艙剪切層存在破壞和重建過程。這種強烈的相互作用一方面使得空艙具有更強的非定常流動特性；另一方面，在彈體穿越剪切層的過程中，受艙體內(nèi)強非定常流場及剪切層的影響，彈體受力出現(xiàn)較大幅度的波動，彈體的下落姿態(tài)也發(fā)生較大的變化，這一階段很大程度上決定了整個發(fā)射過程的成敗。為了保護艙內(nèi)的機載設(shè)備和投放安全，同時保證彈體分離品質(zhì)，對彈體投放過程內(nèi)埋式艙體與彈體間相互影響進行深入研究顯得尤為重要。

現(xiàn)階段采用風洞試驗技術(shù)研究內(nèi)埋投放主要有可控軌跡試驗(CTS)和風洞自由飛試驗兩種，其中CTS基于準定常假設(shè)，對于強耦合條件下的分離問題預(yù)測，無法獲得強非定常非線性效應(yīng)。自由飛技術(shù)可以克服上述問題，但是難以精確控制投放分離的初始條件，同時無法獲得瞬態(tài)流場的精細化流動結(jié)構(gòu)。由此可見，由于受試驗設(shè)備和試驗安全的限制，采用試驗手段研究彈體投放過程中彈體和彈艙的強耦合流動特性具有一定的局限性。而數(shù)值模擬方法則有其優(yōu)越性，尤其是一些新型湍流模擬方法的出現(xiàn)，使得計算精度有較大幅度的提高。本文將通過求解三維非定常Navier-Stokes方程，采用基于SSTk-ω模型的SST-IDDES方法，并利用重疊網(wǎng)格技術(shù)[11-12]和求解剛體動力學方程實現(xiàn)彈體的運動，研究艙門打開后到彈體發(fā)射完成整個過程中，彈體與彈艙之間強烈的相互作用。

1 數(shù)值計算方法

內(nèi)埋式艙體和彈體的內(nèi)外流場控制方程包括連續(xù)方程、動量方程、能量方程、氣體狀態(tài)方程、湍流模型方程和幾何守恒方程(SCL)[13-14]。限于篇幅，這里僅介紹本研究中所用到的SSTk-ω模型和SST DES模型、控制方程的離散方法、重疊網(wǎng)格法及計算方法驗證。

1.1 Menter SST k-ω模型

剪應(yīng)力修正的SSTk-ω模型混合了k-ω模型和k-ε模型的優(yōu)勢，在近壁面處使用k-ω模型，而在邊界層外部區(qū)域采用k-ε模型。SSTk-ω模型包含了修正的湍流黏性公式，考慮了湍流剪切應(yīng)力的效應(yīng)，能更精確地模擬逆壓梯度引起的分離和分離區(qū)大小，SSTk-ω模型公式[15-19]為

(1)

(2)

式中：ρ為密度；k為湍動能；t為時間；U為速度向量；Ug為網(wǎng)格移動速度，可以通過SCL方程進行求解得到；μ為分子黏性系數(shù)；μt為湍流黏性系數(shù)

(3)

(4)

Pk為湍動能生成項；ω為湍流耗散比，β=0.075；α=α1F1+α2(1-F1)，α1=5/9,α2=0.44；σk1=0.85；σω1=0.85；S為應(yīng)變率張量；

(5)

y為網(wǎng)格格心到最近壁面的距離，ν為運動黏度，

(6)

Pk=min(μtS2,10Cμρkω)；

(7)

Cμ=0.09；σω2=0.856.

1.2 SST DES模型

Spalart將大渦模擬(LES)方法和一方程 (S-A)湍流模型相結(jié)合，用S-A湍流模型和LES方法分別模擬以耗散和大渦輸運為主要流動特征的區(qū)域，并通過當?shù)鼐W(wǎng)格尺度與壁面距離的比較和判斷，實現(xiàn)這兩種模擬方法的自動轉(zhuǎn)換，構(gòu)建了分離渦模擬(DES)方法，被習慣性稱為DES97[20-21]。由于網(wǎng)格分布不當時，會產(chǎn)生因雷諾應(yīng)力不匹配而形成的模化應(yīng)力損耗，甚至出現(xiàn)非物理的分離。為了解決這一問題，2006年，Spalart等通過對邊界層參數(shù)進行控制，延遲了RANS湍流模型的作用，提出了Delayed-DES(DDES)模型[22]。

2008年Shur等[23]通過引入壁面?；鬁u模擬(WMLES)成功克服了對數(shù)律不匹配問題，并重新定義亞格子長度尺度為

ΔIDDES=min{max[cwdw,cwhmax,hwn],hmax},

(8)

提出了Improved-DDES。式中：hwn是垂直壁面方向的網(wǎng)格步長；dw為到壁面距離；cw為經(jīng)驗常數(shù)，取0.15；hmax為hwn的最大值。Strelets[24]提出以SST模型代替S-A模型的SST-DES模型，使得新的混合模型兼有了SST模型的諸多優(yōu)點。2012年Gritskevich[25]對SST DDES和IDDES模型的常數(shù)進行了進一步的精確校準，下述文中的SST DES均指的是SST-IDDES模型。

1.3 平滑偽Wigner-Ville分布時頻分析法

無論是彈體投放中彈體出艙穿越剪切層還是彈體出艙后頭部加速區(qū)與剪切層相互干擾，都會使得彈艙內(nèi)部壓力脈動發(fā)生較大變化而改變艙體聲學特性。由于彈體下落的影響，艙體內(nèi)部脈動壓力數(shù)據(jù)屬于非平穩(wěn)信號，因此不能采用傳統(tǒng)傅里葉變換來處理彈體下落過程中艙內(nèi)的壓力脈動，此時必須采用時頻分析技術(shù)進行艙內(nèi)監(jiān)測點壓力脈動的頻譜特性分析。

對于非平穩(wěn)信號，現(xiàn)在已有許多時頻分析技術(shù)，如短時傅里葉變換(STFT)、小波變換(WT)、Wigner-Ville分布(WVD)及平滑偽Wigner-Ville分布(SPWVD)。STFT方法存在時間分辨率和頻率分辨率的矛盾，導致時頻聚集性很差。WT時頻聚集性比STFT優(yōu)異，但是計算復雜耗時很長，不利于較大數(shù)據(jù)量的計算。WVD屬于二次型時頻描述法，由于二次型時頻方法固有特性，嚴重的交叉干擾對其應(yīng)用范圍有較大限制。SPWVD方法繼承了WVD的優(yōu)點同時對交叉干擾進行了平滑處理，因此在計算量和精度上都是適當?shù)摹１疚募催x用SPWVD方法對彈體下落過程中艙體內(nèi)部壓力非平穩(wěn)信號進行時頻分析，以得到艙體聲學特性變化。信號x(t)的SPWVD定義[26]為

SPWVDx(t,f)=

(9)

式中：g(u)h(τ)為實對稱窗函數(shù)，其長度可以獨立進行選擇。

1.4 重疊網(wǎng)格和離散格式

重疊網(wǎng)格法中存在兩套或多套相互交疊的網(wǎng)格，分別稱為背景網(wǎng)格區(qū)域和重疊網(wǎng)格區(qū)域。在計算中通過將背景網(wǎng)格中被重疊區(qū)遮擋的網(wǎng)格隱藏，通常稱之為“挖洞”，形成一套計算網(wǎng)格。如果重疊網(wǎng)格區(qū)域發(fā)生移動，則按一定規(guī)律重新進行“挖洞”，來得到新的計算網(wǎng)格，即可實現(xiàn)重疊網(wǎng)格區(qū)域的運動。背景網(wǎng)格與重疊網(wǎng)格相交的邊界區(qū)域通過一定的方法進行插值來實現(xiàn)變量值的傳遞，并保證數(shù)值解在交界處光滑過渡。

控制方程對流項的空間離散采用具有2階精度的Roe格式，擴散項采用中心差分格式。選取修正的Venkatakrishnan[27]限制器保證2階精度插值且具有TVD性質(zhì)，同時又具有較小的數(shù)值耗散。非定常計算的時間離散涉及控制方程中每一個與Δt有關(guān)的項。物理量φ隨時間變化描述為

F(φ)包含所有的空間離散項。計算中引入雙重時間步法，即在控制方程中引入虛擬時間項，利用物理時間步求解真實解，而每一物理時間步通過虛擬時間迭代達到收斂。

1.5 數(shù)值方法驗證

1.5.1 對空腔噪聲計算方法的驗證

開式流動腔體會在腔體唇口處形成剪切層，剪切層由于自身不穩(wěn)定性或外界擾動形成脫落渦，脫落渦向后傳播并與腔后壁碰撞，產(chǎn)生強烈的噪聲形成聲波，聲波通過空腔內(nèi)部向上游傳遞，與腔體前壁相撞后激發(fā)了剪切層內(nèi)新的渦生成、發(fā)展與脫落，再次與后壁相撞，形成二次聲波，在滿足一定的空氣動力條件和相位條件時，腔內(nèi)流動形成自激振蕩。自激振蕩產(chǎn)生后會形成聲波反饋回路，誘發(fā)腔內(nèi)產(chǎn)生強烈噪聲，出現(xiàn)多個聲壓級峰值的激振頻率。

為了研究空腔氣動聲學特性，歐洲QinetiQ組織對M219空腔進行了多輪風洞試驗得到了可信的試驗數(shù)據(jù)。M219腔體為典型開式空腔[28]，尺寸：長L為0.508 m，寬W和深D均為0.101 6 m，長深比L/D=5∶1. 模型幾何尺寸如圖1所示。

圖1 M219模型幾何形狀及監(jiān)測點分布Fig.1 Sketch of M219 cavity anddistribution of monitoring points

模型的計算網(wǎng)格數(shù)為1 300萬，腔體內(nèi)網(wǎng)格尺寸為1 mm，第一層邊界層網(wǎng)格取2×10-3mm,以保證y+值在1左右。來流馬赫數(shù)Ma=0.85，來流靜溫T=266.53 K，來流靜壓p=63 000 Pa. 與試驗相同，在腔體底面等間距的布置了10個壓力監(jiān)測點。計算時間步選取1×10-5s，總的計算步為20 000步，總采樣時長為0.2 s. 對空腔流動，可用Rossiter[29]、Heller等[30-31]改進的Rossiter半經(jīng)驗公式來估算各階振蕩頻率：

(10)

式中：m=1,2,3，…；r、k為常數(shù)分別取0.25、0.57；Ma∞為遠方來流馬赫數(shù)；U∞為遠方來流速度。由(10)式求得的各階頻率為148 Hz、363 Hz、578 Hz.

對得到的壓力時域分布曲線進行傅里葉變換，得到各測點的聲壓級頻譜圖，如圖2所示。從圖2中可以看出求得的聲壓級頻譜曲線與試驗測得的結(jié)果及公式計算得到的結(jié)果三者吻合得很好，證明本文所采用的數(shù)值方法和網(wǎng)格分布可以較為準確地求解空腔噪聲問題。

圖2 腔底聲壓級分布Fig.2 SPL along the center line of cavity bottom

利用SPWVD時頻分析技術(shù)對驗證算例中計算結(jié)果進行時頻分析。取后緣監(jiān)測點10的壓力時域曲線進行處理，得到其能量密度時頻分布如圖3所示。圖3中顏色條即表示能量密度SPWVDx(t,f)，其單位是Pa2/(s·Hz)，較高的能量密度意味著較強的噪聲水平。由圖3可知能量密度主要集中分布在3個頻率附近，分別對應(yīng)著M219腔體的1～3階模態(tài)，其中1階模態(tài)能量密度較弱，2階、3階模態(tài)能量密度較強，說明腔體內(nèi)噪聲以2階、3階模態(tài)為主導，1階模態(tài)噪聲水平較弱，這與驗證算例得到的聲壓級頻譜圖相一致。同時可以看出各階模態(tài)的能量密度隨著時間增加并不保持為常數(shù)：在0～0.04 s時間內(nèi)能量密度主要分布在1階、2階模態(tài)，而0.04 s之后，1階、2階模態(tài)能量密度減弱，轉(zhuǎn)為集中在第3階模態(tài)。在0.14 s，2階模態(tài)能量密度又有所增加、3階模態(tài)能量密度有所減弱。這一現(xiàn)象即文獻[32]中所提到的腔體噪聲主導模態(tài)隨時間的轉(zhuǎn)換。通過以上分析可以看出，腔體內(nèi)的壓力脈動信號嚴格意義上是非平穩(wěn)信號，因此通常采用的傅里葉變換是一種近似處理，而SPWVD方法可以精確地分析腔體內(nèi)氣動聲學時頻特性。

圖3 M219腔體底部能量密度時頻分布圖Fig.3 Time-frequency distribution of energy density on the bottom of M219 cavity

1.5.2 彈體投放問題的驗證

機翼/掛架/帶舵外掛物(WPFS)模型是美國發(fā)起的第一次分離投放計算流體力學計算驗證時使用的試驗?zāi)Ｐ?，CTS試驗測試由阿諾德工程發(fā)展中心(AEDC)完成，具體的幾何尺寸及結(jié)果數(shù)據(jù)詳見參考文獻[33]。本文以該模型試驗來驗證采用本文所選用數(shù)值方法并應(yīng)用重疊網(wǎng)格技術(shù)和剛體動力學方程模擬彈體投放過程的準確性。

WPFS模型幾何外形如圖4所示，模型網(wǎng)格劃分為兩個區(qū)域，包含機翼及掛架的背景網(wǎng)格區(qū)域和包含彈體的重疊網(wǎng)格區(qū)域。背景網(wǎng)格和重疊網(wǎng)格區(qū)域交界處網(wǎng)格尺寸之比保持在1.0～1.2，以保證插值的精度。背景網(wǎng)格區(qū)域網(wǎng)格總數(shù)為1 150萬，重疊網(wǎng)格區(qū)域網(wǎng)格總數(shù)為320萬。壁面第一層網(wǎng)格尺度為2×10-3mm. 計算時間步取1×10-3s.

圖4 WPFS模型幾何外形Fig.4 Wing-pylon-store geometry

彈體發(fā)射時受到兩個彈射力的作用，大小分別為10 679.4 N和42 717.5 N，彈射力作用距離為0.1 m，作用時間約為0.05 s，詳細參數(shù)見表1. 總的計算時間為0.3 s.

表1 彈體參數(shù)及發(fā)射條件Tab.1 Missile parameters and launching conditions

圖5為3個方向彈體質(zhì)心位移隨時間變化曲線。從圖5中可以看出，計算得到的彈體位移變化與試驗符合較好，由于彈射力作用，彈體下落方向位移要遠大于另外兩個方向的位移。圖6為彈體滾轉(zhuǎn)、偏航、俯仰姿態(tài)角隨時間變化曲線，同樣由于彈體彈射力作用，在下落初期，彈體俯仰角度增長最快，但撤去彈射力后，在重力和氣動力共同作用下俯仰角又開始減小，其余兩個方向姿態(tài)角保持持續(xù)增加的趨勢。可以看出本文采取的數(shù)值方法可以準確地模擬彈體下落運動。

圖5 彈體質(zhì)心位移隨時間變化Fig.5 Displacement of missile Cg versus time

圖6 彈體姿態(tài)角隨時間變化Fig.6 Missile attitude angle versus time

2 計算模型和網(wǎng)格劃分

2.1 計算模型

為了更接近真實地模擬內(nèi)埋艙彈體投放過程，本文選取的艙體計算模型為等比放大的M219腔體模型。艙體尺寸為4.57 m×0.914 m× 0.914 m，彈艙長深比保持5∶1. 彈體模型即為WPFS中彈體，其長度為3.38 m，彈徑0.508 m，質(zhì)心位于距彈頭1.42 m處，彈體質(zhì)量為907.2 kg，x、y、z軸3個方向的轉(zhuǎn)動慣量分別為27 kg·m2、488 kg·m2、488 kg·m2.

選擇的坐標系為x軸為逆航向，z軸向上，y軸向右，則按此坐標系定義，彈射力為負，彈體抬頭為正、左偏航為正。

2.2 計算條件

來流馬赫數(shù)Ma=0.85，靜壓p=6.3×104Pa，靜溫T=266.5 K，基于每米的雷諾數(shù)Re=1.35×107. 來流迎角為α=0°，入口條件設(shè)置為遠場自由來流條件，艙體壁面均采用無滑移壁面條件。彈體投放采用彈射投放方式，彈體所受的彈射力作用規(guī)律與WPFS中彈射力作用規(guī)律一致，詳見表1.

2.3 艙體網(wǎng)格尺寸選擇

為了能準確地模擬大尺寸艙體內(nèi)噪聲問題，需要適當?shù)臄?shù)值方法和恰當?shù)木W(wǎng)格密度。對于本文所采用的數(shù)值方法已經(jīng)得到了驗證，但是對于大尺寸艙體內(nèi)網(wǎng)格尺度的選擇，還需進一步確定。本文將通過試算的方式來確定恰當?shù)木W(wǎng)格尺度。在1.5.1節(jié)已經(jīng)看到，通過Rossiter公式可以準確地預(yù)測艙體的各階模態(tài)頻率，因此本文將以Rossiter公式求得的頻率作為依據(jù)來檢驗網(wǎng)格尺寸是否滿足要求，根據(jù)Rossiter公式求得各階模態(tài)的頻率分別為16 Hz、40 Hz、64 Hz. 本文分別選取了艙體內(nèi)4種網(wǎng)格尺寸，分別為16 mm、12 mm、8 mm和6 mm. 壁面第一層網(wǎng)格尺度保持與驗證算例相一致。最終4種情況總網(wǎng)格數(shù)目分別為350萬、700萬、2 400萬和4 200萬。

計算時間步長為5×10-5s，總計算步為20 000步，總物理時間為1 s. 計算中按照相同的方式在艙體底面均勻布置10個壓力監(jiān)測點，來流條件保持與驗證算例相同。取監(jiān)測點10將得到的結(jié)果進行傅里葉變換，得到聲壓級頻譜如圖7所示。

圖7 不同網(wǎng)格密度下艙底聲壓級分布Fig.7 SPLs along the center line of cavity bottom at different mesh density

通過圖7中的對比可以看出，當網(wǎng)格尺度從16 mm，減小到6 mm的過程中，計算得到的聲壓級均與理論公式求得的聲壓級基本相符。然而對于低頻區(qū)域的寬頻噪聲，16 mm和12 mm網(wǎng)格得到的寬頻噪聲偏大，8 mm與6 mm網(wǎng)格的寬頻噪聲結(jié)果較小且相差不大。綜合考慮計算精度、計算資源和計算時間，選取 8 mm的網(wǎng)格作為基準網(wǎng)格尺度。本文后續(xù)艙體網(wǎng)格劃分及計算時間步長選取均采用本小節(jié)選取的尺度。

彈體網(wǎng)格劃分時保證交界處重疊網(wǎng)格尺寸與背景網(wǎng)格尺寸之比小于1.2. 彈體壁面第一層網(wǎng)格厚度同樣保持為2×10-3mm，總網(wǎng)格數(shù)目為450萬。

3 結(jié)果分析

3.1 彈體投放對艙體聲學特性影響

彈體下落算例中，在建立起流場之后，又對艙體帶彈狀態(tài)進行了0.3 s的非定常計算，之后彈體開始下落，下落時長為0.45 s，總計算時間為0.75 s. 取艙體內(nèi)部前緣監(jiān)測點2和后緣監(jiān)測點8的壓力數(shù)據(jù)進行SPWVD時頻分析，得到能量密度時頻分布圖，如圖8所示。

圖8 艙體底部監(jiān)測點能量密度時頻分布圖Fig.8 Energy density distributions on the bottom of cavity

通過能量密度時頻分布圖可以看出，在彈體下落之前的0.3 s內(nèi)，能量密度集中在各階模態(tài)對應(yīng)的3個頻率附近，由于艙體帶彈的影響，頻率值與Rossiter公式求得的3個模態(tài)的頻率16 Hz、40 Hz、64 Hz有小幅偏差。

在0.3～0.6 s這一階段，整個艙體內(nèi)氣動聲學特性完全改變，從前后緣兩個測點都可以看出，能量密度不再分布于各階模態(tài)對應(yīng)的頻率上，而是集中于低頻區(qū)域，后緣監(jiān)測點的能量密度幅值較彈體下落之前存在著大幅度的增加，前緣點能量密度幅值雖然有所增加但幅度較小。這說明由于彈體的下落，艙體內(nèi)噪聲水平有所增高，后緣點噪聲增加幅度更大，噪聲環(huán)境更惡劣。通過圖9的截面馬赫數(shù)分布云圖可以看出，這段時間為彈體穿越艙體唇口剪切層出艙過程，在此階段艙體剪切層被破壞，開式空艙自持振蕩的回路被阻斷。因此艙體內(nèi)氣動聲學特性發(fā)生顯著改變，不再呈現(xiàn)典型的開式艙體氣動特性。同時由于彈體對剪切層的干擾以及彈體下落后艙體容積改變產(chǎn)生的容積效應(yīng)，使得艙體內(nèi)壓力脈動更加劇烈，噪聲強度增大。圖10為艙體內(nèi)部Q值等值面分布圖。從圖10中可以看出，由于彈體下落的影響，艙體內(nèi)渦結(jié)構(gòu)主要分布在艙體后緣區(qū)域，與之相比艙體前部流動較為穩(wěn)定，因此與艙體后緣相比前緣噪聲水平較低。彈體下落后，艙體后緣噪聲水平升高，增強的噪聲前傳使得前緣噪聲水平高于彈體下落前噪聲水平。

圖10 Q準則等值面分布圖Fig.10 Iso-surface contribution of Q criterion

在0.60 s之后，能量密度又再次集中在3個模態(tài)對應(yīng)的頻率附近。但是艙體后緣點各階頻率上的能量密度較彈體下落前都有所增大。當彈體頭部穿過剪切層之后，彈體完成出艙過程，艙體剪切層迅速重建，自持振蕩回路再次形成，艙體流動狀態(tài)恢復為典型的純空艙流動狀態(tài)。從圖9(c)的馬赫數(shù)分布云圖中可以看出，由于彈體前緣頭部的加速氣流與剪切層相互作用，使得剪切層不穩(wěn)定性增強，同時會使得剪切層內(nèi)渦傳播速度增加，脫落渦具有更大的動能與后緣撞擊，導致艙體后緣壓力脈動更加劇烈，因此艙體后部噪聲水平升高。

3.2 艙體對彈體投放的影響

內(nèi)埋彈體投放過程中存在彈體穿越剪切層過程，此時彈體受到的氣動力呈現(xiàn)出強烈的非定常特性，對彈體的姿態(tài)產(chǎn)生較大的影響。

圖11、圖12分別為彈體豎直方向受力及速度隨時間變化曲線，圖13、圖14、圖15為彈體下落過程俯仰力矩、俯仰角速度及俯仰角隨時間變化曲線。從0.30 s～0.35 s彈體受彈射力作用，此時彈體所受合力較大，具有很大的加速度，短時間內(nèi)速度由0增加到-3.5 m/s，由于彈射力的作用，使得彈體在這一階段內(nèi)受抬頭力矩的作用，抬頭角速度迅速增加，但由于彈射力作用時間短，彈體迎角約1.6°. 0.35 s時彈射力消失，彈體在重力和氣動力共同作用下下落出艙。在0.35～0.42 s時，彈體受力較平穩(wěn)，向下合力維持在10 kN左右，彈體下落的加速度減小，速度增幅下降但繼續(xù)保持穩(wěn)定增長。圖16為0.38 s艙體中心截面馬赫數(shù)分布云圖及彈體表面壓力分布云圖，圖16 (a)中截面為馬赫數(shù)分布，彈體表面為壓力分布，圖16 (b)為彈體下表面壓力分布云圖，圖16 (c)為彈體上表面壓力分布云圖。圖17為彈體上下表面中心線壓力系數(shù)分布曲線。由圖17可知：彈體此時具有一定的迎角，彈體頭部所在區(qū)域流動速度較低，壓力較高；彈體尾部區(qū)域更靠近剪切層，所以流動速度較高，壓力較低。因此彈體此階段受到較小抬頭力矩的作用，彈體俯仰角速度小幅度增加，彈體迎角持續(xù)增大，到0.42 s時，彈體迎角達到7°左右。

圖11 彈體豎直方向受力隨時間變化Fig.11 Vertical force of missile versus time

圖13 彈體俯仰力矩隨時間變化Fig.13 Pitch moment of missile versus time

圖14 彈體俯仰角速度隨時間變化Fig.14 Pitch angular velocity of missile versus time

圖15 彈體俯仰角隨時間變化Fig.15 Pitch angle of missile versus time

圖16 t=0.38 s時艙體截面馬赫數(shù)及彈體表面壓力分布云圖Fig.16 Mach number contours of cavity section and pressure distribution of missile for t=0.38 s

圖17 t=0.38 s時彈體上下表面對稱線壓力系數(shù)分布Fig.17 Pressure coefficients of symmetric lines on upper and lower surfaces for t=0.38 s

圖18 t=0.49 s時截面馬赫數(shù)及彈體表面壓力分布云圖Fig.18 Mach number contours of cavity section and pressure distribution of missile for t=0.49 s

圖19 t=0.49 s時彈體上下表面對稱線壓力系數(shù)分布Fig.19 Pressure coefficient distribution of symmetric lines on upper and lower surfaces for t=0.49 s

從0.42 s起，彈體開始穿越剪切層，圖18 (a)為0.49 s時艙體中心截面馬赫數(shù)分布云圖及彈體表面壓力分布云圖。圖18(b)、圖18(c)分別為彈體下表面和上表面壓力分布云圖。圖19為彈體上下表面中心線壓力系數(shù)分布曲線。從圖18 (a)中可知，彈體頭部位置剪切層受到下落彈體的阻礙，撞擊在彈體下表面，一部分動能轉(zhuǎn)變?yōu)閴耗?。從圖18 (b)中可以看出，彈體頭部下表面受到剪切層沖擊的區(qū)域局部壓力顯著增高。從圖18(b)、圖18(c)及圖19中可以看出彈體上表面壓力明顯低于下表面。彈體頭部下表面的局部高壓區(qū)使得彈體俯仰力矩保持為抬頭力矩，上下表面較大的壓力差導致氣動力合力方向為豎直向上，且隨著彈體出艙過程不斷增大。因此在彈體穿越剪切層階段，彈體俯仰角速度及抬頭角度持續(xù)增加，且彈體下落受到較大的阻礙，下落速度增速減緩。

隨著彈體進一步下落，彈體頭部穿過剪切層，作用在彈體上的抬頭力矩減小；同時彈體受抬頭角速度的影響，彈體迎角持續(xù)增加，0.52 s時，彈體迎角增加到14°，彈體后部下表面壓力繼續(xù)增加，尤其是彈翼下表面壓力及高壓區(qū)面積也明顯增大。受此影響，最終使彈體合力矩變?yōu)榈皖^力矩，且隨著彈體的迎角增加、低頭力矩持續(xù)增大。當俯仰力矩變?yōu)榈皖^力矩后，彈體抬頭角速度開始減小，但迎角持續(xù)增加。到0.53 s時向下的合力減小為0，隨后變?yōu)橄蛏系暮狭Γ瑥楏w下落速度不增反降。

直到0.67 s之后，彈體俯仰角速度減小到0并進一步變?yōu)榈皖^角速度，彈體迎角達到最大值后開始減小，彈體的豎直方向合力開始下降，迎角的減小也使得彈體受到的低頭力矩逐漸降低。

4 結(jié)論

1) 彈體投放過程中內(nèi)埋艙體與彈體存在強烈的相互作用。彈體下落出艙穿越剪切層會使艙體唇口剪切層被破壞，導致開式艙體流動狀態(tài)發(fā)生改變。艙體自持振蕩回路被阻斷，艙內(nèi)噪聲不再呈現(xiàn)明顯的各階單調(diào)聲模態(tài)。受彈體的影響，艙體內(nèi)渦結(jié)構(gòu)主要集中于艙體后緣，且強度較高使艙體后緣噪聲水平升高。

2)彈體出艙后，艙體唇口處剪切層迅速重建，艙體恢復典型空腔流動特性。但是受彈體頭部加速氣流的影響，剪切層不穩(wěn)定性增強，脫落渦強度增高，與艙體后壁撞擊產(chǎn)生的噪聲更大，使得艙體內(nèi)部聲壓級水平總體上升。

3)彈體下落過程受艙體強非定常流場的影響，彈體受力及力矩均出現(xiàn)較大幅的波動。彈體穿越剪切層時頭部受剪切層撞擊產(chǎn)生局部高壓區(qū)，使彈體產(chǎn)生抬頭力矩導致彈體姿態(tài)角發(fā)生改變出現(xiàn)抬頭。彈體姿態(tài)角的抬頭又使彈體受力發(fā)生改變，氣動力阻礙彈體下落，影響飛機與彈體安全分離。

References)

[1] 馮必鳴, 聶萬勝, 車學科. 超聲速條件下內(nèi)埋式武器分離特性的數(shù)值分析[J]. 飛機設(shè)計, 2009, 29(4): 1-5. FENG Bi-ming, NIE Wan-sheng, CHE Xue-ke. Simulation of the store separation from a cavity at supersonic speed[J]. Aircraft Design, 2009, 29(4)：1-5. (in Chinese)

[2] 楊黨國, 范召林, 李建強, 等. 超聲速空腔流激振蕩與聲學特性研究[J]. 航空動力學報, 2010, 25(7): 1567-1572. YANG Dang-guo, FAN Zhao-lin, LI Jian-qiang. Cavity flow oscillations and aeroacoustic characteristics at supersonic speeds [J]. Journal of Aerospace Power, 2010, 25(7): 1567-1572. (in Chinese)

[3] 謝露, 艾俊強, 李權(quán), 等. 長深比對空腔流動與聲學特性的影響分析[J]. 航空工程進展, 2014, 5(1): 18-24. XIE Lu, AI Jun-qiang, LI Quan, et al. Effect of length-to-depth ratio on flow and aero acoustic characteristics of cavity [J]. Advances in Aeronautical Science and Engineering, 2014, 5(1): 18-24. (in Chinese)

[4] 徐路, 桑為民, 雷熙薇. 三維內(nèi)埋式彈艙流動特性及形狀影響數(shù)值分析[J]. 應(yīng)用力學學報, 2011, 28(1): 85-89. XU Lu, SANG Wei-min, LEI Xi-wei. Numerical analysis of flow characteristics and shape effect of three dimensional internal weapons bay [J]. Chinese Journal of Applied Mechanics, 2011, 28(1): 85-89. (in Chinese)

[5] 王一丁, 陳濱琦, 郭亮, 等. 空腔噪聲非線性數(shù)值模擬[J]. 國防科技大學學報, 2015, 37(4): 151-157. WANG Yi-ding, CHEN Bin-qi, GUO Liang, et al. Nonlinear numerical simulation of cavity noise [J]. Journal of National University of Defense Technology, 2015, 37(4): 151-157. (in Chinese)

[6] 張群峰,閆盼盼,黎軍.邊界層厚度對腔體氣動聲學特性影響數(shù)值模擬[J].航空動力學報, 2016, 31(3): 717-725. ZHANG Qun-feng, YAN Pan-pan, LI Jun. Numerical simulation on influence of boundary-layer thickness to cavity aeroacoustic characteristics [J]. Journal of Aerospace Power, 2016, 31(3):717-725. (in Chinese)

[7] 楊黨國,羅新福,李建強.來流邊界層厚度對開式空腔氣動聲學特性的影響分析[J]. 空氣動力學學報, 2011, 29(4): 486-490. YANG Dang-guo, LUO Xin-fu, LI Jian-qiang. Analysis of aeroacoustic characteristics in open cavities influenced by boundary-layer thickness[J]. Acta Aerodynamica Sinica, 2011, 29(4): 486-490. (in Chinese)

[8] 黎軍, 李天, 張群峰. 開式流動腔體的流動機理與控制[J]. 實驗流體力學,2008,22(1):80-83. LI Jun, LI Tian, ZHANG Qun-feng. The mechanism and control of open cavity flow[J].Journal of Experiments in Fluid Mechanics, 2008, 22(1): 80-83. (in Chinese)

[9] 黃長強, 國海峰, 唐上欽, 等. 超聲速帶彈武器艙氣動特性數(shù)值研究[J]. 兵工學報,2013, 34(8): 975-980. HUANG Chang-qiang, GUO Hai-feng, TANG Shang-qin, et al. Numerical research on aerodynamic characteristics of cavity with store at supersonic speeds[J].Acta Armamentarii, 2013, 34(8): 975-980. (in Chinese)

[10] 司海青, 王同光. 數(shù)值模擬有外掛物的空腔流動[J].空氣動力學學報,2007, 25(3): 404-409. SI Hai-qing, WANG Tong-guang. Numerical simulations of the cavity with a store [J]. Acta Aerodynamica Sinica,2007, 25(3):404-409. (in Chinese)

[11] 朱自強. 應(yīng)用計算流體力學[M]. 北京：北京航空航天大學出版社, 1998: 173-174. ZHU Zi-qiang. The application of computational fluid dynamics [M]. Beijing: Beijing University of Aeronautics and Astronautics Press, 1998: 173-174. (in Chinese)

[12] 閻超. 計算流體力學方法及應(yīng)用[M].北京：北京航空航天大學出版社，2006: 197-217. YAN Chao. The computational fluid dynamics method and its application [M]. Beijing, Beijing University of Aeronautics and Astronautics Press, 2006: 197-217. (in Chinese)

[13] Tamura Y, Fujii K. Conservation law for moving and transformed grids[C]∥Proceedings of 6th National Symposium on Computational Fluid Dynamics. Tokyo, Japan: Japan Society of Computational Fluid Dynamics, 1993: 519-522.

[14] 劉偉. 細長機翼搖滾機理的非線性動力學分析及數(shù)值模擬方法研究[D]. 長沙: 國防科學技術(shù)大學, 2004. LIU Wei. Nonlinear dynamics analysis for mechanism of slenderwing rock and study of numerical simulation method [D]. Changsha： National University of Defense Technology, 2004. (in Chinese)

[15] Durbin P A, Pettersson B A. Statistical theory and modeling for turbulent flows[M]. 2nd ed. New York: A John Wiley Sons Ltd, 2011.

[16] Wilcox D C. Turbulence modeling for CFD[M]. La Canada, CA: DCW Industries, 1998: 142-148.

[17] Menter F R. Two-equation eddy-viscosity turbulence modeling for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[18] Menter F R, Kuntz M. A daptation of eddy viscosity turbulence models to unsteady separated flow behind vehicles[M]∥McCallen R, Browand F, Ross J. The aerodynamics of heavy vehicles: trucks, buses and trains. Berlin: Springer Berlin Heidelberg, 2004:339-352.

[19] Menter F R, Kuntz M, Langtry R. Ten years of industrial experience with the SST turbulence model [C]∥ Proceedings of the Fourth International Symposium on Turbulence, Heat and Mass Transfer. Danbury, CT: Begell House, 2003: 625-632.

[20] Spalart P R. Detached eddy simulation [J]. Annual Review of Fluid Mechanics, 2009, 41(1):203-229.

[21] Spalart P R, Jou W H, Strelets M, et al. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[C]∥1st AFOSR international Conference on DNS/LES. Columbus, OH: AFOSR, 1997.

[22] Spalart P R, Deck S, Shur M L, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densities[J]. Theoretical and Computational Fluid Dynamics, 2006, 20(3): 181-195.

[23] Shur M L, Spalart P R, Strelets M K, et al. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities [J]. International Journal of Heat and Fluid Flow, 2008, 29(6): 1638-1649.

[24] Strelets M. Detached eddy simulation of massively separated flows [C]∥39th AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV US: AIAA, 2001.

[25] Gritskevich M. Development of DDES and IDDES formulations for thek-ωshear stress transport model [J]. Flow, Turbulence and Combustion, 2012, 88(3): 431-449.

[26] Auger F, Flandrin P. Improving the readability of time-frequency and time-scale representations by the reassignment method [J]. IEEE Transactions on Signal Processing, 1995, 43(5): 1068-1089.

[27] Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions[C]∥ 31st Aerospace Sciences Meeting and Exhibit. Reno, NV, US: AIAA, 1994:152-164.

[28] Henshaw M J C. M219 cavity case: verification and validation data for computational unsteady aerodynamics, RTO-TR-26[R]. London: QinetiQ, 2002.

[29] Rossiter J E. Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds [R]. Farnborough: Aeronautical Research Council, 1964.

[30] Heller H H, Holmes D G, Covert E E. Flow induced pressure oscillations in shallow cavities [J]. Journal of Sound and Vibration, 1971, 18(4):545-546.

[31] Heller H H, Bliss D B. The physical mechanism of flow-induced pressure fluctuations in cavities and concepts for their suppression [C]∥2nd AIAA Aeroacoustics Conference. Hampton,VA,US: AIAA, 1975.

[32] Zhuang N. Experimental investigation of supersonic cavity flow and their control [D]. Tallahassee, Florida, US: Florida State University, 2007.

[33] Heim E R. CFD wing/pylon/finned store mutual interference wind tunnel experiment, AEDC-TSR-91-P4 [R]. Tennessee: Arnold Engineering Development Center, 1991.

Elaborate Simulation of Interaction Effect between Internal Weapon Bay and Missile

ZHANG Qun-feng1，YAN Pan-pan1，LI Jun2

(1.School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China;2.Shenyang Aircraft Design and Research Institute,Aviation Industry Corporation of China，Shenyang 110035, Liaoning, China)

To study the strong interaction effect between internal weapon bay and missile during the weapon release, the unsteady flow around a simplified weapon bay and missile model is simulated using SST k-ω improved delayed detached eddy simulation (IDDES) method, six-degrees-of-freedom rigid body dynamics equations and overset mesh method. The fluctuation of pressure in the weapon bay is analyzed using the smooth pseudo Winger-Vile distribution (SPWVD) method, and the time-frequency characteristics of pressure are obtained. The research results show that, when a missile is delivered from the weapon bay, the shear layer is destroyed to lead to the change in the flow structure, the self-sustained oscillation disappears in the weapon bay, and there is no obvious cavity tone. Stronger vortexes concentrate in the trailing edge of cavity, and thus the sound pressure level (SPL) in the trailing edge of the weapon bay increases. After the missile is delivered from the weapon bay, the shear layer is rapidly reestablished and the self-sustained oscillation shows up again. The impact of the accelerated flow near the head of the missile leads to enhance the instability in shear layer. Accordingly, the sound pressure level in the weapon bay is enhanced. The missile flies at a large angle of attack since it is affected by upward force moment during passing through the shear layer. After the missile is delivered from the weapon bay, it begins to be affected by vertical upward force and pitch down moment. The vertical upward force can hinder the missile from falling. Meanwhile, the force and moment acting on the missile fluctuate strongly under the effect of the unsteady flow in the weapon bay.

fluid mechanics; internal weapon bay; weapon release; overset mesh; SST DES model; smooth pseudo Winger-Vile distribution

2016-05-03

國家自然科學基金項目(11172283)

張群峰(1972—)，男，講師，博士。E-mail: zhangqunfeng@263.net

V211.3

A

1000-1093(2016)12-2366-11

10.3969/j.issn.1000-1093.2016.12.024