李曉東,姜蔚婉,2,屈 崑,*,蔡晉生
(1. 西北工業大學 航空學院,西安 710072;2. 中國航空工業集團 成都飛機設計研究所,成都 610091)
內埋式武器在投放分離過程中,彈艙將演變為(帶艙門)大尺度空腔。大尺度空腔是航空航天武器裝備必不可少且普遍采用的一種典型結構形式,主要用于作戰武器裝載、動力燃油存儲和起落架存放等,具有尺度大、結構復雜等特點。在動載荷激勵下,空腔薄壁結構容易產生嚴重的振動和噪聲,不僅影響武器裝備的作戰效能,甚至危及武器裝備安全。因此,空腔繞流是空氣動力學領域的一個重要且典型問題。
針對空腔中流動和噪聲的復雜現象與物理機理,國內外科研工作者通過風洞試驗和計算手段對干凈彈艙(空腔)模型[1-4]進行了廣泛深入的研究。具體包括分析流致振蕩現象的成因[5-6],對空腔類型[7]、來流馬赫數、邊界層厚度等[3,8]參數進行詳細研究,并探討了多種不同的控制措施[2,9],對空腔設計具有重要的指導意義。Rossiter[10]提出了預測空腔噪聲模態頻率的半經驗模型;Heller 等[11]改進了此模型,將其推廣到高馬赫數(Ma= 2、3)的可壓縮流動中。
空腔是彈艙的一種理想模型,實際中艙門的存在會使得空腔內非定常流動更加復雜,因而研究學者們對帶艙門的空腔流動開展了研究和分析。首先是對帶有靜態艙門的空腔流動的研究:Blair 和Stallings[12]采用風洞試驗研究了不同艙門開度和艙門高度對導彈類物體在超聲速空腔中投放的影響,結果表明艙門開度對彈體的氣動載荷有顯著影響,并且降低艙門高度會減少彈體的氣動載荷;Lawson 和Barakos[13]采用脫體渦模擬(DES)方法對艙門開啟90°的空腔標模M219 和1303 無人機模型進行了模擬,并分析了艙門打開對彈體拋投的影響。Casper 等[14]采用試驗的方式模擬了帶有艙門的復雜空腔中彈體拋投的流固耦合問題,發現增加空腔復雜度會引發彈體展向共振。其次在動態艙門的研究上,Sheta 等[15]采用RANS/LES混合方法,模擬了超聲速(Ma= 1.44)空腔艙門從5°打開到35°過程中的流動,發現在艙門運動到30°時,剪切層和艙內聲波產生較強的相互作用,腔內壁面出現強的壓力脈動;Loupy 等[16]對比分析了固定開度艙門和動態艙門對馬赫數為0.85 的空腔流動、載荷和氣動噪聲的影響,發現艙門打開過程中的過渡階段,艙門氣動載荷增大,流動的脈動和總聲壓級增強,對空腔結構產生顯著影響。
國內多位學者采用RANS 方法模擬分析靜態艙門[17]和艙門開啟過程中[18]空腔動態氣動特性和載荷的變化。吳繼飛等[19-20]對Ma= 0.6 空腔艙門的開閉問題開展了研究,設計了不同的運動方式,發現艙門運動時氣動載荷發生明顯變化,艙門打開30°時氣動載荷最大,并分析了腔內和艙門的總聲壓級和聲壓譜特性。由于帶艙門的空腔流動非定常效應顯著,且流動的湍流結構復雜,所以傳統的RANS方法捕捉復雜流動結構比較困難。閆盼盼等[21]采用基于k?ω的分離渦模擬方法(IDDES)計算了艙門不同開啟姿態對內埋武器投放的影響,發現艙門會改變超聲速流場中的波系結構,對彈體偏航影響較大。
考慮運動艙門時,空腔中氣動與噪聲特性會變得更加復雜,預測難度增加,而國內對帶運動艙門的復雜空腔流動與噪聲的研究較少。為了研究艙門運動情況下的空腔復雜流動現象,深入分析靜、動態艙門對空腔復雜流動的影響,本文采用嵌套重疊網格方法對艙門和空腔分別生成網格, 采用改進的脫體渦模擬方法對亞聲速(Ma= 0.6)空腔流場進行數值模擬,從流場特性、湍流結構、脈動特性上分析艙門對湍流流場和腔內噪聲的影響。


針對空腔流動強脈動、強湍流的特點,采用RANS/LES 混合方法進行計算,本文采用改進的脫體渦湍流模擬方法(Improved Delayed DES,IDDES)[22],解析空腔內分離區域的湍流結構和脈動特性?;谝陨系臄抵捣椒?,本文發展了一套基于有限體積方法的多塊結構網格求解器,并通過多種復雜流動問題驗證其準確性和可靠性[23-26]。

式中,

為驗證高精度數值格式在空腔問題中的有效性,采用本文提出的格式對馬赫數0.85 空腔標模M219進行模擬分析,并與參考文獻[32]的LES 計算結果進行對比。網格單元總數量約800 萬,空腔內部與附近的網格數量約380 萬。對于馬赫數為0.85 的流動,時間步長為 Δt=2×10?5s, 約為Tref=L/U∞的百分之一,足以捕捉第一階Rossiter 頻率。圖1 給出了來流馬赫數為0.85 時,中截面上時均流向速度、縱向速度和湍動能的分布,結果與參考文獻符合得較好。由此可見,本文采用的計算方法可以比較準確地捕捉空腔復雜湍流的時均統計特性。

圖1 空腔標模M219 在Ma = 0.85 時均速度與湍動能分布Fig. 1 Wall-normal profiles of the mean velocity and turbulence kinetic energy in canonical cavity M219 at Ma = 0.85
針對運動物體的網格生成問題,發展了快速的動態嵌套重疊網格方法[33]。首先通過采用分層嵌套重疊網格策略(如圖2)對網格進行劃分,然后采用自適應樹結構完成網格切割和貢獻單元搜索。引入壁面距離概念,設計了統一的分層嵌套重疊策略,同時實現了重疊網格切割和嵌套網格切割兩個功能,完成重疊最小化的功能,使得重疊網格組裝質量得到了改善。

圖2 分層嵌套重疊網格策略Fig. 2 Hierarchically embedded strategy of overset grids
根據運動規律的不同,可將所有運動物體分為若干動態組,而所有靜態物體歸為唯一的靜態組。各組分別構造獨立的二叉樹和八叉樹輔助網格,并且每個組內的流場網格和二/八叉樹網格的坐標均保持不變。即,各個動態分組仍然用各自原有的局部坐標系,但在每個時間步上記錄下各動態組運動的位移。而在更新切割信息時,將所有被切割的坐標點變換到切割物體所在的局部坐標內,從而完成運動問題的網格切割。對于同一組內的切割,由于相對位置不變,只需進行一次切割,后續不需再改變。此方法使得二叉樹和八叉樹在運動過程中不必重新構造,只需進行坐標變換操作,大大節省了計算時間。表1 給出了艙門運動時網格組裝時間的對比。對于總計2 100 多萬的網格單元來說,兩種動態嵌套重疊網格組裝方法均可滿足該類非定常流動計算的效率要求。

表1 考慮運動艙門的空腔動態重疊網格組裝的耗時統計Table 1 Assembly time of overset grids for a cavity with moving doors
空腔標模C201 由中國空氣動力研究與發展中心設計,空腔底部、前后壁分別布置了靜壓測點和脈動壓力測點(見圖3)。本文采用的計算模型如圖4(a)所示,與試驗的模型外形一致。空腔長度(L)為200 mm,深度(D)為33.3 mm,寬度(W)為66.7 mm,空腔長深比L/D為6,寬深比W/D為2;與空腔相連壁面的長度和寬度分別是Wa和La。

圖3 空腔標模C201 中截面的壓力測點Fig. 3 Distribution of pressure probes in the central section of canonical cavity C201
參照文獻[1]的數據選取艙門,艙門厚度1.5 mm,艙門四周采用弧形物面進行過渡。此構型與空腔標模M219 的艙門外形有一定相似性。在艙門與空腔、艙門與艙門之間保留1 mm(約為艙門深度的3%)縫隙,這樣可以保證縫隙處保留足夠的網格用于數值模擬。艙門三維模型如圖4(b)所示。

圖4 空腔標模C201 與艙門計算模型Fig. 4 Geometry of canonical cavity C201 and doors
整個模型網格采用一個網格層,包含三個網格簇:空腔網格簇(其中包含背景網格)和兩個艙門網格簇。圖5 給出了空腔的網格,由于網格在展向(z方向)是對稱的,所以圖中顯示了一半網格。邊界條件設置如下:在圖5 給出的模型中,紅色為無滑移物面邊界條件;綠色表示的是空腔前后的自由端,采用對稱邊界條件;展向的兩個邊界均采用對稱邊界條件;最上面的外場邊界設置為入口/出口邊界;前后外場分別為入口/出口邊界。空腔網格采用結構網格,共計約1 600 萬網格單元。

圖5 空腔標模C201 計算模型與網格(每4 個點顯示一個網格點)Fig. 5 Computational domain and grids for canonical cavity C201
艙門采用沿著物面法向將物面網格向外推進的方法生成貼體網格,然后在外圍采用矩形外場,使得外場網格的周向尺寸得到控制。單個艙門的網格單元共計約454 萬,圖6(a)給出了艙門的物面網格,艙y+=1門的前向剖視圖如圖6(b)所示。艙門壁面的第一層網格厚度由 給出。艙門開度的定義如圖6(c)所示。

圖6 艙門計算網格與艙門打開角度(開度)定義Fig. 6 Computational grids for cavity doors and the definition of the door-opening angle
剪切層流動是空腔復雜流場的重要因素之一,所以需要保留高質量的空腔剪切層網格。通過改進空腔壁面距離的計算方法,使得在四種不同艙門打開角度下均能保留空腔剪切層網格(見圖7),從而保證剪切層流動計算的精細和準確。

圖7 四種不同艙門開度在 x/L=0.5處yOz 視圖上網格的組裝結果Fig. 7 The grid slice of assembled overset grids in yOz planes for cavities with different door-opening angles
對無艙門的空腔標模C201 亞聲速流動(Ma=0.6)進行模擬,來流條件為T∞=286.66 K,Re=12.14×106。 時 間步 長 采用 Δt=1.0×10?6s,約為 參 考 時 間Tref=L/U∞的千分之一,足以分辨流動的細節。采用雙時間步推進求解非定常流場,每步偽時間推進采用多重網格加速收斂,并采用MPI 并行技術提高計算效率。非定常計算采用IDDES 模擬方法,在廣州天河二號超級計算機上進行。采用SA 湍流模型計算的穩態流場作為初場,在非定常流場充分發展之后(約15Tref~20Tref),開始非定常流場的數據采集和時均場的累計。時均場的累計大約為 50Tref,約為0.05 s,足以分辨第一階Rossiter 頻率。
圖13(a)給出了中截面時均流場的壓強系數分布和流線形狀,在空腔后部存在一個較小的第二環流。在較低速來流條件下,剪切層厚度增長較大,大范圍沖擊到后壁、角落和底部。沖擊后壁與角落的氣流較為靠近剪切層核心區,流速較高,形成局部高壓環流;而沖擊底部的氣流較為遠離剪切層核心區,流速較低,被角落的高壓環流所阻礙,提前回流形成了主環流。
圖8 對比了三種不同網格尺度下腔內底部壓力系數的分布與試驗值的比較。隨著網格的加密,計算結果趨于穩定。盡管其值高于試驗結果,但是計算結果可以較好地預測壓力系數的變化。從圖中可以看出,底部壓強分布表現為前段平緩下降、中段增長和后段快速上升,最后氣流沖擊后壁形成高壓區域。

圖8 三種網格對腔內底部壓力系數時均分布的影響Fig. 8 Grid convergence of the mean pressure coefficients at cavity floor compared to experimental distribution
總聲壓級(overall sound pressure level,OASPL)可以通過湍流模擬得到的壓力脈動來獲得,其定義如下:


圖9 空腔標模C201 腔內底部總聲壓級Fig. 9 OASPL distribution at the floor of canonical cavity C201 without doors

圖10 空腔標模C201 底部聲壓級成分合成示意圖Fig. 10 Sketch of the development of OASPL in canonical cavity C201
通過脈動壓力功率譜分析空腔底部的頻譜特性,對脈動壓力p′采用快速傅里葉分析得到功率譜密度(PSD), 然后將PSD 用SPL[13]呈現,進行功率譜的分析。圖11 給出了空腔標模C201 腔內底部第3 和第11 監測點脈動壓強功率譜,圖中灰色條帶是根據Heller 公式計算的峰值頻率帶。計算結果與試驗結果吻合,主導峰值頻率都是第二峰值,因此在工程上需要特別注意第二峰值頻率的影響。

圖11 空腔標模C201 腔內底部p3 和p11 測點脈動壓強功率譜Fig. 11 Comparison of SPL of fluctuating pressure at probes p3 and p11 between numerical simulation and experiment
在滿足準定常假設時,可以對靜態艙門直接模擬分析。采用發展的嵌套重疊網格方法和IDDES 湍流模擬方法,對艙門開度分別為30°、60°、90°和120°的四種工況進行準定常流動模擬,分析其流動時均特性和流場動態特性。從時均流場(中截面時均流場和底部壓力分布)、流場的動態特性(總聲壓級、脈動壓強功率譜)和瞬時場湍流結構三個方面分析艙門對空腔流場產生的影響。本節采用的數值方法、計算方式與干凈空腔模擬方法一致。
將空腔底部的靜壓測點時均壓力與干凈空腔的結果進行對比分析,圖12 給出了不同狀態下底部壓力系數的分布。從腔內底部壓力分布來看,艙門打開30°時,底部壓力系數的分布與干凈空腔的差別較大,前部壓強有所升高,后部壓強明顯降低。在艙門打開60°、90°、120°時,腔內前半部分底部壓力系數基本與干凈空腔的計算結果一致,后半部分壓力系數隨著艙門打開角度的增加逐漸趨于干凈空腔。開啟角度為90°和120°時,時均壓強系數分布與干凈空腔基本一致。

圖12 四種艙門開度(30°、60°、90°、120°)對空腔底部壓力系數的影響Fig. 12 Impacts of the door-opening angle on mean pressure coefficients at the cavity floor
圖13 給出了干凈空腔和不同艙門開度的中截面時均流場壓力系數云圖和流線分布情況。從圖中可以看出在艙門打開角度為30°時,空腔前緣附近的平板上,壓力系數會升高,這是由于艙門的阻礙使得流動減速造成的。這種壓強升高作用也使得空腔前部的負壓區域的壓強系數相對干凈空腔要偏高。同時,對于30°開啟角度,外流與空腔附近的流動受到艙門的阻擋,只能通過中截面附近的開口相互作用。特別是阻礙外流對剪切層的能量供給,使得后壁附近的壓強相對干凈空腔明顯降低,導致30°時,空腔底面后部壓強系數水平整體降低。此外,后壁附近的高壓也受到艙門阻礙,在縱向方向只能在中截面開口處泄壓,后部流線有更加明顯地向上彎曲。

圖13 干凈空腔和四種艙門開度下中截面平均流場壓力系數與流線分布Fig. 13 Impacts of door-opening angles on the streamlines and mean pressure contours in the middle section of cavity C201
艙門開度為60°、90°和120°時,空腔前緣平板同一位置的壓力系數逐漸降低,且后壁高壓區和后壁拐角處低壓區的范圍較艙門開度為30°時逐漸變大,并與干凈空腔的結果趨于一致。這與“艙門打開角度越大,對流動的限制明顯越小”的物理原因是一致的。
另外,相比于干凈空腔,隨著艙門打開角的不斷增大,主渦位置沿流向不斷地后移并靠近干凈空腔主渦位置。艙門打開120°時,湍流結構基本上與干凈空腔的一致,但是在前壁附近的流場結構與干凈空腔有所區別。由此可知,隨著空腔打開角度的增加,腔內時均流場的形態、渦的結構與干凈空腔的結果逐漸趨于一致。
圖14 給出了四種艙門開度下空腔底部總聲壓級的分布。從圖中可以看出,艙門開度30°時,底部各個位置的總聲壓級較干凈空腔均明顯下降,但變化趨勢與干凈空腔的結果基本一致。艙門打開90°時,總聲壓級與120°的變化趨勢一致,前者的值略高1~2 dB。艙門開度為60°時,總聲壓級腔內前半部分的分布特點與90°、120°的結果基本一致;而在腔內后半部分,其總聲壓級明顯有所改變。
本文認為導致噪聲升高的原因在于艙門阻礙了腔內噪聲向外傳播。從圖14 可以看出,艙門分別開啟60°、90°、120°時,三種條件下的聲壓級分布幾乎一致,因此可知噪聲升高與艙門開啟角度無關,只與艙門的存在有關。艙門的存在使得空腔內部的脈動能量無法以聲波的形式沿展向傳播,只能沿著兩艙門開口方向傳播,第二階模態的能量無法得到有效釋放,這樣就導致腔內噪聲強度升高。本文中艙門開度為90°時,總聲壓級的計算結果與空腔標模M219 在馬赫數0.85 下的試驗結果[1]分布特點相似,呈現出“W”形態的分布。

圖14 干凈彈艙與四種艙門開度下的彈艙底部總聲壓級分布Fig. 14 Impacts of door-opening angles on OASPL at the cavity floor
空腔中的噪聲主要由寬頻白噪聲和窄頻噪聲(Rossiter 模態)兩部分組成[13]。白噪聲是由流動結構產生,其來源為自由來流的脈動、剪切層和湍流等。第二部分也稱為Rossiter 模態,是不同流場結構直接相互作用產生的,其特點為幅值大且集中在某幾個頻率。圖15 給出了四種不同艙門開啟角度下,腔內9 個脈動測點的脈動壓強功率譜曲線,柱狀條帶是Heller 公式計算的峰值頻率。從圖中可以看出:

圖15 不同艙門開度下腔內9 個脈動測點脈動壓強功率譜Fig. 15 SPL of fluctuating pressure at nine probes on cavity floors for cases with different door-opening angles
1)艙門開度為30°時,二階Rossiter 所占優勢顯著減弱,并且該模態頻率偏高,其他高頻(第三、第四)模態已經不再明顯。這是由于艙門打開較小時,空腔流動中剪切層與回傳壓縮波的相互作用受到限制,所以流場結構相互作用的強度減弱,致使聲壓幅值降低。
2)60°、90°和120°艙門開度下,各脈動壓力測點的Rossiter 模態明顯,且模態頻率與經驗公式基本一致。這說明此時流場中的結構與干凈空腔一致。二階Rossiter 模態占主導,幅值明顯比其他模態高10 dB以上。
3)圖中虛線、實線分別表示腔內前、后半部分采樣點的結果。從圖中可以看出,四種不同艙門開啟角度下,后半部分的脈動強度明顯高于前半部分。
為研究艙門對腔內脈動特性的影響,進一步選取后壁(p11)的脈動壓強測點上的脈動壓強功率譜進行比較分析,如圖16 所示。各圖中黑色曲線為干凈空腔試驗曲線,褐色曲線為干凈空腔計算曲線,其他彩色曲線為各角度情況下的計算曲線。

圖16 不同艙門開度腔內后壁測點p11 脈動壓強功率譜Fig. 16 SPL of fluctuating pressure at probe p11 on cavity floors for cases with different door-opening angles
對于后壁測點p11 來說,當開啟角為30°時,功率水平明顯低于干凈空腔,其主峰頻率比干凈空腔的主峰頻率偏高。這可能是由于小角度開啟情況下,空腔構型更接近封閉腔體,其展向或者縱向的共振成分加強。而在其他開啟角度情況下,變化基本一致,峰值頻率大致和干凈空腔一致,但是第二模態峰值都顯著升高,三、四模態峰值沒有升高甚至幅值降低。這也印證了第二階模態被強化的觀點。同時這也說明,當打開角度較大時,在各處,艙門對主模態的影響是一致的,是通過改變流場的主結構來影響噪聲的產生。在前壁和腔內底部測點的脈動壓強功率譜中,出現了相似的變化和影響,不在此贅述。
整體來看,出現不同的變化的原因是,在艙門打開小角度時,空腔內的法向(y向)和展向(z向)流動均受到限制,流場脈動會受到抑制,噪聲也會減小。艙門打開角度較大(60°以后)時,艙門的限制主要在展向,對法向的限制已經不明顯;但是此時艙門的存在會增加空腔流動的復雜性,阻斷空腔近壁面噪聲沿展向向外傳播,從而會增強腔內的噪聲。
為了進一步分析不同艙門開度下空腔流動的特點,通過Q等值面比較瞬時場中湍流結構的異同。正的Q值區域表示渦張量比應變率張量大的流動區域,是與旋渦有關的部分,是常被用作判斷旋渦區域的標準[34]。圖17 給出了干凈空腔和四種不同艙門開度下,瞬時流場中的Q= 3 000 等值面,采用來流方向的無量綱速度u進行染色。

圖17 瞬時流場的Q = 3 000 等值面(用u 進 行染色)Fig. 17 Instantaneous iso-surfaces (extracted by Q = 3 000) colored by u for clean cavity and the cavity with stationary doors
從瞬時圖中可以看出,由于艙門的存在,空腔內的流動會受到不同程度的影響,具體表現為:
1)30°的湍流結構明顯受到艙門的限制,從而使得流場的湍流結構較干凈空腔要小,艙內渦結構很稀疏,尤其是空腔前部區域。
2)在艙門打開60°時,渦結構變得密集,可以明顯看到艙內的流動結構沿法向(y向)向外運動,但是空腔前部仍然比較稀疏。
3)到90°和120°時,渦結構的豐富程度變得相似。在艙門開度為120°時,艙門仍然會限制空腔展向的流動。這從側面反映了艙門的存在會阻礙空腔周圍噪聲的傳播,從而使得腔內噪聲一定程度地增大。
4)相比干凈空腔的瞬時Q準則等值面圖,可以看出,有艙門時空腔上方的剪切層渦結構運動高度更高,這仍然是剪切層在展向受艙門限制的結果。而隨著艙門開度的不斷增大,展向的空腔流動會受到不同程度的限制。
5)流動經過前緣平板后形成的剪切層,在流場結構中會出現大尺度的剪切層渦結構。艙門對腔內前半部分湍流結構的相互作用較弱;對后腔內半部分不同流動結構之間相互作用的影響較大,從而影響空腔流動的脈動和噪聲特性。
為進一步模擬真實艙門打開過程對空腔流動和聲場的影響,本節采用前文的空腔和艙門,在考慮艙門運動的情況下,對空腔的湍流流動進行模擬。對考慮運動的非穩態非定常流動問題,采用經驗模態分解(empirical mode decomposition,EMD)方法提取非穩定特征,定量分析其噪聲特性的變化。
艙門運動規律的選取由條件相似準則確定,而斯特勞哈爾數(St)是非定常流動中須考慮的相似準則,物理上代表非定常運動與定常運動的慣性力之比,其定義為:

式中,LD為艙門長度, ΔT為艙門每次開啟的總時間,V∞表示來流速度。戰斗機武器艙艙門開啟過程中,需1 ~2 s 打開100°[35]。為使得艙門運動規律對流動的影響與真實飛行情況一致,需要保證St數相同。因此,在來流條件相同的情況下,當模型縮比為1∶10 時,在縮比模型模擬過程中,艙門需要0.1~0.2s左右完成一次動作。本文選取0.1 s,艙門運動的角速度為900 rad/s;從艙門完全關閉到打開120°,模擬的時間步長取 Δt=2×10?5。
首先對艙門關閉狀態進行非定常計算,保證流動充分發展;然后艙門開始運動,進行完整的非穩態非定常流動模擬。圖18 給出了艙門從0°打開到100°過程中,空腔內動態測點的壓力系數變化。從圖中可以看出,整個流動大致分為三個階段:初期的小幅穩定階段(≤30°)、過渡階段(30°~60°)、充分發展階段(≥60°)。同時空腔前半區域各動態采集點(p3~p7)的壓力系數呈對稱擺動趨勢,而空腔后半區域動態采集點的壓力系數呈增長趨勢。以p11 測點為例,圖19展示了其壓力值的變化情況。

圖18 艙門打開過程中各動態采集點的壓力系數變化Fig. 18 Variations of pressure coefficients at different probes during the opening of cavity doors

圖19 艙門開啟過程中彈艙后壁測點p11 處壓力變化與EMD 模態曲線Fig. 19 Variation of pressure and its EMD modes at probe p11 during the opening of cavity doors
由于艙門運動使得空腔內流動呈現非穩態的變化,即壓力脈動(p′)需要使用局部平均值而不是簡單的代數平均。為了準確描述艙門運動過程的非穩態變化,本文提出采用希爾伯特-黃轉換(Hilbert-Huang Transform) 中 的 經 驗 模 態 分 解( Empirical Mode Decomposition,EMD)進行分析。該方法的核心在于尋找當前數據中的所有局部極大值以及局部極小值,然后分別采用三次樣條構建極大值包絡線和極小值包絡線,兩者的平均值即為EMD 模態;再對扣除EMD 模態的剩余數據進行分析,直到剩余數據為單調函數為止。
為證明該方法的有效性,采用艙門打開90°的準定常流動進行分析,并用代數平均的計算結果進行對比驗證。以測壓點p11 為例,圖20 給出了不同EMD 模態的分布曲線。分析結果表明,低階EMD 模態分析包含較多高頻信息,與原始數據分布具有一致性,無法直接用于捕捉非穩態的信息;高階的EMD 模態多數與低頻流場信息有關,可用于描述艙門運動引起的非穩態流場的變化。從圖20 中,EMD5 和EMD6模態的分布差異可以觀察到此特征。而更高階的EMD 模態(如EMD9)基本與代數平均值一致。從而說明EMD 方法具有與代數平均計算方式相同的功能。

圖20 艙門打開90°時p11 點壓力變化與不同EMD 模態曲線Fig. 20 Time histories of pressure and its EMD modes at probe p11 in the cavity with vertical doors
圖21 中給出了采用不同EMD 模態作為時均值時,p11 點總聲壓級的變化。結果表明,采用低階EMD 模態所得的OASPL 遠低于采用代數平均值方式計算的值。而提高采用的EMD 模態,OASPL 會逐漸收斂于代數平均計算的值(相差在0.5%以內)。這說明采用EMD 模態進行OASPL 計算是有效可行的。盡管EMD 模態的選取無統一標準,但可以采用OASPL 的收斂性判斷選取的EMD 模態。

圖21 艙門打開90°時考慮不同EMD 模態的p11 點OASPL 計算結果Fig. 21 OASPL computed using different EMD modes at probe p11 in the cavity with vertical doors
圖22 給出了EMD 模態對p11 測點功率譜的影響,主要特點為低頻區域有差別,高頻區域無明顯差異。由于EMD 本身就是對高頻濾波的結果,不同的EMD 曲線代表不同的低頻模態,所以會對低階頻率產生影響。計算結果中的典型Rossiter 頻率與采用平均值計算的方式一致,可見采用EMD 模態對壓力功率譜的計算方式是可行的。
采用經驗模態分解方法,對亞聲速空腔艙門開啟過程中的非穩態動態特性進行分析。圖23 給出了空腔后壁p11 測點的OASPL 隨著選取不同EMD 模態所產生的變化。當取到EMD5 模態時,計算得到的OASPL 開始趨于穩定,收斂值為艙門運動過程中的OASPL。值得注意的是,該值低于靜態艙門的計算值。這是因為,在艙門的開啟過程中,當艙門打開角度小時,流動的脈動較小,艙門打開后期的脈動變大(見圖19),空腔內流動與聲場逐漸產生相互作用并增強,進而產生了流動的振蕩現象。圖24 給出了采用艙門開啟到90°附近的數據進行分段動態分析的OASPL 計算結果。此時的結果與靜態艙門計算值相當,說明對靜態艙門打開90°的準定常計算跟真實運動過程中的狀態吻合,這對工程問題的分析具有指導意義。

圖23 艙門開啟過程中不同EMD 模態的p11 點OASPL 計算結果Fig. 23 OASPL computed using different EMD modes at probe p11 during the opening of cavity doors

圖24 艙門開啟到90°時不同EMD 模態的p11 點的OASPL 分析結果Fig. 24 OASPL computed using different EMD modes at probe p11 when cavity doors are opened to 90°
圖25 給出了采用不同EMD 模態時,對應的脈動壓力功率譜分析結果。艙門運動過程中,二階Rossiter 模態仍為主導模態,幅值明顯高于其他模態10 dB 以上,但低于靜態艙門打開90°的分析結果,此處與總聲壓級降低的原因一致。一階Rossiter 模態強度降低,而三、四階Rossiter 模態略有升高。

圖25 艙門開啟過程中不同EMD 模態的p11 點的功率譜比較Fig. 25 SPL computed using different EMD modes at probe p11 during the opening of cavity doors
采用EMD 方法對空腔內全部脈動測點進行OASPL 分析,如圖26 所示。計算結果表明,艙門開啟過程中,OASPL 分布狀態仍呈現“W”形狀,但其值比靜態大開度艙門整體降低5 dB 左右。這是由于艙門開啟過程中,存在由小幅脈動到強脈動的發展過程。小角度的艙門開啟過程中,流場的脈動較小,還沒有完全發展;強脈動的完全發展需要一定的時間。在艙門的整個運動過程中流動呈現過渡發展的特點,因此出現OASPL 的降低。

圖26 艙門開啟過程中空腔底部總聲壓級分布Fig. 26 OASPL at cavity floor during the opening of cavity doors
圖27 給出了艙門開啟過程中,艙門打開到不同角度時,物面壓力系數的分布情況。從圖中可以看出,壓力分布具有一定對稱性,兩側艙門承受的壓力相當。在艙門打開小角度時,艙門限制了空腔內部與外部流場的摻混與交換,壓力變化較小。隨著艙門的逐漸打開,空腔與外流場發生混合,空腔內壓力變化增大。艙門打開角度較大(60°以后)時,空腔內流動出現大幅度振蕩現象,流場與聲場出現強相互作用,從而造成流場脈動增強,噪聲特性變得顯著。與前文OASPL 分析和功率譜分析結果一致。

圖27 艙門開啟過程中物面壓力系數變化Fig. 27 Instantaneous wall pressure coefficients during the opening of cavity doors
針對考慮艙門的多尺度、強湍流、高雷諾數、強脈動的空腔流動問題,采用本文發展的嵌套重疊網格方法對帶運動艙門的空腔生成網格,用IDDES 湍流模擬方法對亞聲速(Ma= 0.6)空腔流動進行模擬,分析了靜、動態艙門對空腔復雜流動和噪聲特性的影響。具體結論如下:
1)發展的動態嵌套重疊網格方法,可自動高效生成大幅度運動物體的計算網格,用于運動艙門的空腔復雜流場模擬。
2)采用四階對稱重構和三階MUSCL 混合的格式,在不明顯增加計算量的前提下,明顯降低了數值耗散,提升了計算的準確度,可以捕捉空腔中的湍流結構。對標模M219 的計算結果與文獻給的結果吻合。對空腔標模C201 腔內噪聲的預測和對脈動壓力功率譜的分析結果與試驗數據一致。
3)亞聲速空腔標模C201 腔內底部壓強分布表現出三段分布:前段平緩負壓區;中段壓力增長過渡區;后段高度增加區,形成后壁的高壓區域。腔內主導峰值頻率為第二峰值。
4)與干凈空腔相比,當艙門為靜態時,艙門小開度(30°)時,時均流場中截面上腔內的流態變得更加復雜,并且腔內底部壓力系數的值均會明顯減小。艙門開度較大(60°以后)時,時均流場中截面的腔內渦結構與干凈空腔一致,底部壓力分布趨于干凈空腔的計算結果;但是艙門的存在,會使得底部后半部分的壓力較干凈空腔出現一定程度上的降低。艙門開度為90°時,腔內總聲壓級較干凈空腔有所升高,最大升高值近9 dB,呈現出“W”形態的分布;二階Rossiter模態占主導。艙門與空腔后半部分的不同流動結構之間相互作用,從而對空腔流動的脈動和噪聲特性產生影響。
5)針對艙門運動的非穩態非定常特點,本文采用經驗模態分解方法對噪聲特性進行分析,計算結果表明該方法有效可行。艙門開啟過程中,整體的OASPL有所下降,主導頻率依然為二階Rossiter 模態。對于小開度的艙門,靜態艙門模擬與運動艙門結果相差較大,因而此類情況需要考慮運動艙門的影響。在艙門開度較大(90°)時,兩者的分析結論一致,說明可采用準定常來分析此狀態的流場與噪聲特性,指導內埋彈艙的設計。