基于Sobol’法的可變形狀空腔全局靈敏度分析

寧方立，張暢通，劉 哲，丁 輝，翟慶波，韋 娟

(1.西北工業(yè)大學(xué) 機電學(xué)院，西安 710072；2. 東莞市三航軍民融合創(chuàng)新研究院，廣東 東莞 523808；3. 西安電子科技大學(xué) 通信工程學(xué)院，西安 710071)

空腔結(jié)構(gòu)廣泛存在于航空飛行器中，比如戰(zhàn)斗機的內(nèi)埋彈艙、飛機的起落架艙、導(dǎo)彈的光學(xué)頭罩等，當飛行器處在高速飛行狀態(tài)時，由于氣流和空腔結(jié)構(gòu)的耦合作用，腔體內(nèi)會產(chǎn)生劇烈的壓力脈動，從而產(chǎn)生高強度噪聲，損壞空腔結(jié)構(gòu)和空腔內(nèi)的電子設(shè)備。在過去的數(shù)十年中，空腔噪聲問題受到了廣泛的關(guān)注[1-5]，眾多被動或主動控制方法用來研究抑制空腔噪聲。

周方奇等[6-8]試驗研究了前緣擾流板對空腔噪聲的影響，結(jié)果表明，擾流板能夠影響剪切層在空腔前緣處的發(fā)展，從而能夠有效地降低飛行器上空腔結(jié)構(gòu)的噪聲。Pey等[9]對不同后壁面傾斜角度的空腔流動進行了試驗研究，結(jié)果表明后壁面傾斜可以抑制空腔噪聲。吳繼飛等[10]試驗研究了后壁倒角對開式空腔氣動噪聲的抑制效果，試驗結(jié)果表明，后壁倒角能夠有效抑制艙底的總聲壓級，并且抑制效果隨著馬赫數(shù)的增大愈加明顯。楊黨國等[11-14]通過在空腔前緣安裝射流裝置對空腔噪聲進行主動控制，結(jié)果表明，射流會改善空腔內(nèi)的流場，對空腔噪聲有很好的抑制效果。

本文提出一種兩自由度的可變形狀空腔控制方法，并用Sobol’法研究這種控制方法的全局靈敏度。前期的工作[15-16]研究一自由度的可變形狀空腔的流場及聲場特性。二自由度的可變形狀空腔的示意圖，如圖1所示。二自由度的可變形狀空腔與單自由度的可變形狀空腔相比，在圖1的位置1處增加一個活動鉸鏈，可以實現(xiàn)腔體前、后壁面的同時傾斜，其中前壁面的傾斜角度為α，后壁面的傾斜角度為β。該控制方法具有以下優(yōu)點：① 同時利用空腔的前、后壁面及底板的傾斜來抑制空腔噪聲；②可以根據(jù)實際的工況來實時調(diào)整前、后壁面的傾斜角度。

圖1 可變形狀空腔 (二自由度)Fig.1 Deformable cavity (two degrees of freedom)

對二自由度的可變形狀空腔進行全局靈敏度分析。由于空腔底板結(jié)構(gòu)的固有頻率較低，文獻[17]中空腔底板的一階固有頻率為196.8 Hz，空腔噪聲的不同頻率成分對結(jié)構(gòu)振動的作用效果不同，在噪聲的主模態(tài)位置，結(jié)構(gòu)振動與表面噪聲載荷的主模態(tài)相關(guān)性最強，因此本文選擇空腔底板表面噪聲載荷的主模態(tài)幅值為預(yù)測對象，輸入變量為空腔的前、后壁面的傾斜角度。采用計算流體力學(xué)的方法計算空腔底板表面的壓力波動，再利用快速傅里葉變換得到空腔底板表面噪聲載荷的主模態(tài)幅值。利用Kriging模型動態(tài)建立可變形狀空腔的前、后壁面傾斜角度與空腔底板表面噪聲載荷的主模態(tài)幅值之間的映射關(guān)系。最后，根據(jù)已經(jīng)建立的代理模型，使用Sobol’法分析前、后壁面傾斜角度對空腔底板表面噪聲載荷的主模態(tài)幅值影響所占的比重。

1 數(shù)值方法

在研究空腔問題中，三維空腔的計算量要比二維空腔的計算量大很多，為降低計算量，本文研究的空腔是二維空腔，并選擇大渦模擬(Large Eddy Simulation, LES)湍流求解方法。雖然LES在二維計算中會過度預(yù)測渦的耗散，但是LES在預(yù)測二維空腔流動中仍然可以取得不錯的效果，很多學(xué)者在研究空腔問題中證明這一點[18-20]。

在前期的工作中，將本文的二維LES仿真結(jié)果和參考文獻[21]相對比，驗證本文選用的數(shù)值模型及計算結(jié)果是可靠的。

2 基于MSE準則的Kriging模型和Sobol’法

基于MSE準則的Kriging模型和Sobol’法，首先需要通過Kriging模型建立輸入?yún)?shù)和輸出參數(shù)之間的映射關(guān)系。然后，利用Sobol’法分析輸入?yún)?shù)對輸出參數(shù)影響所占的比重。

2.1 Kriging模型

Lophaven等[22]編寫Matlab工具箱DACE，該工具箱在工程領(lǐng)域得到廣泛的運用，采用模式搜索算法優(yōu)化相關(guān)參數(shù)。游龍海等[23]和陳鵬等[24]分別應(yīng)用遺傳算法和粒子群算法優(yōu)化Kriging模型的相關(guān)參數(shù)，使得相關(guān)參數(shù)的數(shù)值更加接近全局最優(yōu)解。因此本文選擇DACE工具箱，并用粒子群算法優(yōu)化相關(guān)參數(shù)，采用MSE加點準則[25- 27]動態(tài)建立Kriging代理模型。

Kriging代理模型由回歸項和隨機項兩部分組成

y(x)=F(x)+α(x)

(1)

式中：F(x)為回歸項，階次有0階、1階和2階，本文選擇1階；α(x)為隨機項，其期望為0；任意兩點w和x的隨機項的協(xié)方差的期望為

E[α(w)α(x)]=σ2R(θ,ω,x)

(2)

式中：σ2為樣本隨機項的均方；R(θ,ω,x)為w和x之間的高斯相關(guān)函數(shù)，其值為

(3)

式中：n為設(shè)計空間的維度；θj為未知的相關(guān)參數(shù)，其數(shù)值越大代表兩個樣本點之間的相關(guān)性隨著距離的增加衰減的越劇烈。

通過極大似然估計可以確定相關(guān)參數(shù)

(4)

式中：R為樣本點集的相關(guān)系數(shù)矩陣；m為樣本點的個數(shù)。

任意點處的MSE函數(shù)值為該點處的真實值與預(yù)測值差值平方的期望，基于MSE的加點準則需要在模型的MSE函數(shù)值的最大處增加樣本點以改善模型全局精度。

2.2 Sobol’法

Sobol’法是常見的一種基于方差的Monte Carlo法，是俄羅斯學(xué)者Sobol在20世紀90年代提出的一種有效的研究輸入?yún)?shù)的方差對輸出參數(shù)方差影響的方法，關(guān)于Sobol’法的詳細介紹可以參見文獻[28-29]。

因為可變形狀空腔有兩個自由度，所以只需要計算空腔前壁和后壁面的一階靈敏度，就可以得到前、后壁之間的二階靈敏度以及前壁和后壁面的全局靈敏度

S12=1-S1-S2

(5)

(6)

(7)

3 可變形狀空腔的計算模型

可變形狀空腔的底板長度為0.508 m，空腔的前壁面和后壁面的長度均為0.101 6 m，圖2(a)為前壁面傾斜5°、后壁面傾斜50°的可變形狀空腔(算例1)的網(wǎng)格，前壁面傾斜50°、后壁面傾斜5°的可變形狀空腔(算例2)的網(wǎng)格，如圖2(b)所示?？涨粌?nèi)的網(wǎng)格量均為380×200，網(wǎng)格的邊界條件和網(wǎng)格獨立性檢驗工作已經(jīng)在先前的研究中詳細給出。來流Ma=0.85，靜壓為62 940 Pa，靜溫為270.25 K，監(jiān)測點設(shè)置在空腔底板上，距離空腔后壁與空腔底板的交線0.05倍的空腔底板長度。

圖2 可變形狀空腔的計算網(wǎng)格Fig.2 Mesh of deformable cavity

4 可變形狀空腔的流場激勵及聲學(xué)模態(tài)分析

通過基準算例、算例1和算例2流場變化可以探究空腔噪聲的產(chǎn)生機理及可變形空腔的流場激勵對空腔噪聲的影響。

圖3為基準算例流場一周期內(nèi)的渦量云圖。空腔流場中存在著周期性的漩渦生成、移動、碰撞等過程。在空腔的前壁處，邊界層與空腔前緣發(fā)生分離，隨著流體流動，剪切層在空腔前部形成大尺度的渦并與剪切層發(fā)生分離，分離的大尺度渦繼續(xù)向下游運動，在空腔的中部位置與空腔底板發(fā)生分離，最終與空腔后壁發(fā)生碰撞，碰撞過后，漩渦破裂并繼續(xù)向下游運動，逐漸淡出空腔流場區(qū)域，在此期間又會有新的大尺度渦產(chǎn)生。正是源于空腔流場中存在著周期性的漩渦碰撞過程，使得空腔后壁處的靜態(tài)壓力快速振蕩變化，從而產(chǎn)生了空腔噪聲。

圖3 基準空腔的渦量圖Fig.3 Vorticity magnitude of standard cavity

算例1流場一周期內(nèi)的渦量云圖，如圖4所示?？涨涣鲌鲋腥匀淮嬖谥芷谛缘匿鰷u生成、移動、碰撞等過程，但是與矩形空腔流場相比存在以下三點不同：① 由于空腔的平均深度變小，使得渦在“Y”方向的發(fā)展空間變小，從空腔前緣唇口處分離出的漩渦較矩形空腔更小；②由于空腔后壁有一定的傾斜角度，使得后緣角與原矩形空腔相比“鈍”一些，同時由于前、后壁面傾斜，增大漩渦運動的距離，即增加能量耗散；③ 由于后壁面傾斜角度大于前壁面的傾斜角度，使得空腔底板傾斜方向朝著空腔外部，這導(dǎo)致從前緣分離的漩渦向下游運動時，有沿著腔體底板向上抬升的趨勢，由圖4(c)可知，漩渦幾乎完全處于腔體后緣平面之上，而圖3(c)中，漩渦與后緣角有較大面積的碰撞。

圖4 算例1的渦量圖Fig.4 Vorticity magnitude of case 1

算例2流場一周期內(nèi)的渦量云圖，如圖5所示。空腔流場中仍然存在著周期性的漩渦生成、移動、碰撞等過程，但是與算例1相比存在以下三點不同：① 在空腔的前部，空腔的平均深度更小，使得渦在“Y”方向的發(fā)展空間更小，從而使得從空腔前緣唇口處分離出的漩渦較矩形空腔更??；②空腔后壁面傾斜角度較小，使得渦與空腔后壁面的碰撞強度加大；③ 由于前壁面傾斜角度大于后壁面的傾斜角度，使得空腔底板傾斜方向朝著空腔內(nèi)部，這使得渦在向下游運動時，更靠近空腔內(nèi)部，使得渦與空腔的碰撞面積比算例1更大。

圖5 算例2的渦量圖Fig.5 Vorticity magnitude of case 2

上述分析表明算例2改善空腔流場的作用沒有算例1明顯，但是仍然對空腔內(nèi)的流場有改善作用，這是由于前壁面傾斜可以抑制渦在空腔前部的發(fā)展，并且前、后壁面傾斜會使空腔的長度伸長，這增加渦的運動距離，即增加能量耗散，使得運動至后緣的漩渦與后緣的碰撞力度減小。為了進一步探究可變形空腔聲學(xué)模態(tài)的變化，對監(jiān)測點處基準算例、算例1和算例2的噪聲頻譜進行詳細分析。

算例1、算例2和基準監(jiān)測點的噪聲頻譜，如圖6所示。因為空腔噪聲的能量主要集中在低頻成分中，所以只給出了0～600 Hz這個頻率段的噪聲頻譜。圖中n=1和n=2這兩條豎線分別代表Rossiter公式預(yù)測的第一階模態(tài)和第二階模態(tài)的頻率位置，基準仿真的一階模態(tài)頻率位置與Rossiter公式預(yù)測位置相差1.397 4 Hz，二階模態(tài)相差38.000 2 Hz，說明了本文選用的數(shù)值模型及計算結(jié)果是可靠的。一階模態(tài)是監(jiān)測點噪聲頻譜的主模態(tài)，算例1和算例2的主模態(tài)頻率位置相較于基準算例發(fā)生了右移，在先前的工作中已經(jīng)對這種現(xiàn)象給予了充分的分析。算例1和算例2的主模態(tài)峰值相較于基準分別降低了6.67 dB、3.75 dB，說明了可變形狀空腔通過改變空腔的形狀使得空腔內(nèi)的流場發(fā)生改變，進而作用于空腔聲場，使得空腔噪聲得到了有效的抑制。一階模態(tài)作為主模態(tài)，能量最高，并且最接近空腔的固有頻率，因此可以將其與可變形狀空腔前、后壁面的傾斜角度的關(guān)系函數(shù)作為靈敏度分析的目標函數(shù)。

圖6 算1和算例2的噪聲頻譜Fig.6 Noise spectrum of case 1 and case 2

5 可變形狀空腔的靈敏度分析

可變形狀空腔靈敏度分析的過程，如圖7所示。空腔前、后壁面的傾斜角度在5°～60°變化，以監(jiān)測點處噪聲頻譜的主模態(tài)幅值作為響應(yīng)值，設(shè)定模型的迭代終止條件為最大MSE函數(shù)值在0.5 dB2以下。在設(shè)計空間內(nèi)隨機選擇30個樣本點，構(gòu)建擬合曲面和MSE函數(shù)曲面，如圖8所示。從圖8(a)可知，改變可變形狀空腔的前、后壁面的傾斜角度均可以抑制空腔噪聲的主模態(tài)峰值，并且空腔噪聲的主模態(tài)峰值整體上隨著前、后壁面傾斜角度的增加而降低。從圖8(b)中可以看出，距離樣本點較遠的位置MSE函數(shù)值較大。選擇粒子群算法搜索初始模型的MSE函數(shù)，得到最大MSE函數(shù)值為1.602 4 dB2，對應(yīng)的傾斜角度為前壁、后壁都傾斜60°，超過迭代終止條件，因此需要添加樣本點以提高模型的精度。

圖7 可變形狀空腔靈敏度分析過程Fig.7 Sensitivity analysis of deformable cavity

圖8 可變形狀空腔的擬合曲面和MSE函數(shù)曲面(30個樣本點)Fig.8 The surface of deformable cavity and MSE (30 samples)

基于MSE的加點準則需要在模型的MSE函數(shù)值的最大處增加樣本點以改善模型全局精度，因此添加的第一個樣本點的位置為前壁、后壁都傾斜60°。通過仿真計算得到前壁、后壁都傾斜60°算例的主模態(tài)峰值，因此得到了第一個動態(tài)添加的樣本點的值，再結(jié)合初始模型的30個樣本點構(gòu)建這31個樣本點的Kriging模型和MSE函數(shù)。繼續(xù)通過粒子群算法搜索新MSE函數(shù)，最大MSE函數(shù)值仍然超過0.5 dB2，因此需要繼續(xù)添加樣本點，重復(fù)上述操作，直到滿足最新的MSE函數(shù)的最大值不超過0.5 dB2這一迭代終止條件，最終動態(tài)添加了6個樣本點的信息，如表1所示。從表1可知,每次添加新的樣本點后最大MSE函數(shù)值在逐漸降低，這意味著Kriging模型的精度在逐漸提高，添加完6個樣本點后，最大MSE函數(shù)值為0.465 4 dB2，所以動態(tài)添加6個樣本點后，模型的最大MSE函數(shù)值在0.5 dB2以下，可以終止添加樣本點。最終建立的可變形狀空腔的擬合曲面和MSE函數(shù)曲面，如圖9所示。從圖9(a)可知，前、后壁面傾斜對空腔底板表面噪聲載荷的主模態(tài)幅值的影響規(guī)律沒有發(fā)生變化，空腔噪聲的主模態(tài)峰值整體上隨著前、后壁面傾斜角度的增加而降低，從圖9(b)可知，MSE函數(shù)值曲面在0.5 dB2以下。

表1 可變形狀空腔增加的6個樣本點信息

圖9 可變形狀空腔的擬合曲面和MSE函數(shù)曲面(36個樣本點)Fig.9 The surface of deformable cavity and MSE (36 samples)

在設(shè)計空間內(nèi)選擇4個樣本點用來檢驗?zāi)Ｐ偷木?，對比這4個樣本點的CFD預(yù)測值和模型預(yù)測值，如表2所示。從表2可知，4個檢驗樣本點的模型預(yù)測值和CFD預(yù)測值相差不大，這說明模型的精度較好。

表2 檢驗樣本點的主模態(tài)峰值

對可變形狀空腔進行Sobol’全局靈敏度分析，程序采用LHS法采集10萬個樣本點，得到的靈敏度分析結(jié)果如表3所示。從表3可知，后壁面傾斜的全局靈敏度是前壁面傾斜的1.72倍，說明前、后壁面傾斜對抑制空腔噪聲的主模態(tài)峰值均有效果，后壁面傾斜的抑制效果要比前壁面傾斜大很多。

表3 可變形狀空腔的靈敏度分析結(jié)果

6 結(jié) 論

本文用LES仿真法對可變形狀空腔的傾斜角度為5°～60°進行研究，證明了二自由度可變形狀空腔可以有效抑制空腔噪聲。以空腔前、后壁面傾斜角度作為設(shè)計變量，以空腔底板表面噪聲載荷的主模態(tài)幅值作為響應(yīng)，使用基于MSE準則的Kriging模型和Sobol’法研究了可變形狀空腔的靈敏度，得到以下研究結(jié)論：

(1) 隨著可變形狀空腔的前、后壁面傾斜角度的逐漸增大，空腔底板表面噪聲載荷的主模態(tài)幅值整體上是下降趨勢，后壁面傾斜對空腔底板表面噪聲載荷的主模態(tài)幅值的抑制效果更加明顯。

(2) 可變形狀空腔可以增加空腔的長度并減小空腔的平均深度，并且使得空腔后緣角變“鈍”，這有利于改善空腔流動。

(3) 后壁面傾斜的Sobol’全局靈敏度是前壁面傾斜的1.72倍，前、后壁面傾斜的二階靈敏度只有0.051 5，說明前、后壁面傾斜對空腔底板表面噪聲載荷的主模態(tài)幅值的交互影響小。

本文的研究對工程上利用空腔壁面傾斜抑制空腔噪聲提供了理論分析方法，具有工程應(yīng)用前景。