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

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

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

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

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

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

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

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

1 數值方法

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

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

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

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

2.1 Kriging模型

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

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

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

(1)

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

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

(2)

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

(3)

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

通過極大似然估計可以確定相關參數

(4)

式中：R為樣本點集的相關系數矩陣；m為樣本點的個數。

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

2.2 Sobol’法

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

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

S12=1-S1-S2

(5)

(6)

(7)

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

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

圖2 可變形狀空腔的計算網格Fig.2 Mesh of deformable cavity

4 可變形狀空腔的流場激勵及聲學模態分析

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

圖3為基準算例流場一周期內的渦量云圖?？涨涣鲌鲋写嬖谥芷谛缘匿鰷u生成、移動、碰撞等過程。在空腔的前壁處，邊界層與空腔前緣發生分離，隨著流體流動，剪切層在空腔前部形成大尺度的渦并與剪切層發生分離，分離的大尺度渦繼續向下游運動，在空腔的中部位置與空腔底板發生分離，最終與空腔后壁發生碰撞，碰撞過后，漩渦破裂并繼續向下游運動，逐漸淡出空腔流場區域，在此期間又會有新的大尺度渦產生。正是源于空腔流場中存在著周期性的漩渦碰撞過程，使得空腔后壁處的靜態壓力快速振蕩變化，從而產生了空腔噪聲。

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

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

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

算例2流場一周期內的渦量云圖，如圖5所示?？涨涣鲌鲋腥匀淮嬖谥芷谛缘匿鰷u生成、移動、碰撞等過程，但是與算例1相比存在以下三點不同：① 在空腔的前部，空腔的平均深度更小，使得渦在“Y”方向的發展空間更小，從而使得從空腔前緣唇口處分離出的漩渦較矩形空腔更??；②空腔后壁面傾斜角度較小，使得渦與空腔后壁面的碰撞強度加大；③ 由于前壁面傾斜角度大于后壁面的傾斜角度，使得空腔底板傾斜方向朝著空腔內部，這使得渦在向下游運動時，更靠近空腔內部，使得渦與空腔的碰撞面積比算例1更大。

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

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

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

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

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

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

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

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

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

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

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

在設計空間內選擇4個樣本點用來檢驗模型的精度，對比這4個樣本點的CFD預測值和模型預測值，如表2所示。從表2可知，4個檢驗樣本點的模型預測值和CFD預測值相差不大，這說明模型的精度較好。

表2 檢驗樣本點的主模態峰值

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

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

6 結 論

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

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

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

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

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