水下平板流激聲輻射快速預報方法研究

胡昊灝, 周石頭, 趙 傲, 蔣哲倫

(江蘇科技大學 船舶與海洋工程學院,江蘇 鎮江 212000)

水下航行器中高航速行駛時,會在結構表面產生湍流脈動壓力,并引起流激振動噪聲[1]。準確預報航行體流激噪聲,對水下結構減振降噪和聲隱身有重要意義[2-3]。由于這類噪聲涉及到流場、結構、聲場等多重因素,且具有典型時空隨機性,所以研究起來存在難度。

Graham[4]通過模態展開法建立了流激平板聲輻射的解析理論,為掌握流激噪聲產生規律奠定了理論基礎。當水下結構形狀復雜時,解析公式不再適用,Chevalier等[5]總結了水下航行器流激自噪聲和輻射噪聲的工程數值預報方法,小尺度模型一般先采用計算流體動力學 (computational fluid dynamics,CFD)大渦模擬計算結構表面湍流脈動壓力,然后導入聲振有限元模塊計算流激聲輻射,而對于大尺度模型CFD計算效率無法滿足工程要求,一般用半經驗公式描述湍流脈動壓力,呂世金等[6]通過風洞測量獲取了回轉體表面的湍流脈動壓力,發現測試結果與Corcos模型經驗公式有較好的一致性。Maxit等[7]比較了采用空間-頻率譜和波數-頻率譜兩種形式的湍流脈動壓力經驗公式數值求解流激噪聲的計算效率,指出空間-頻率譜有限元法由于要構造大維數激勵力矩陣,對內存占用嚴重(即使采用矩陣分解法依然如此),從而導致計算效率低,而波數-頻率譜則能明顯降低內存占用,同時物理意義更清晰,然而關于波數-頻率譜形式的有限元法研究目前還很少。

Aucejo等[8]提出了一種用于替代風洞中測量流激振動的試驗方案,該方法用若干組空間隨機分布的非相關平面波疊加來替代湍流脈動壓力波數-頻率譜。本文在此基礎上,提出一種基于非相關平面波疊加的方法合成湍流脈動壓力,并結合有限元模塊建立波數-頻率域的流激噪聲數值預報方法,探索了相關參數對計算精度的影響。

1 平板流激聲輻射基本理論

如圖1所示,長為Lx,寬為Ly的矩形平板四邊簡支在無限大剛性障板上,在充分發展的湍流作用下,計算上半空間無限大聲場的流激輻射噪聲,這里假設彈性板的振動不會對流場產生影響。

圖1 流激平板聲輻射示意圖

此時彈性板的運動方程為

D?4w(x,y)-ρshω2w(x,y)=-pt(x,y)+pa(x,y)

(1)

式中:w為板法向位移;ρt為作用在板表面的湍流脈動壓力;pa為板表面聲壓負載。通過解析法可以求得彈性板振速譜Svv(x,ω)和輻射聲功率譜Πrad(ω)

(2)

Πrad(ω)=

(3)

式中:ΦPP(k,ω)為湍流脈動壓力波數-頻率譜;Hp(x,k,ω),Hp(x,k,ω)分別為振動傳遞函數以及聲壓傳遞函數,其物理意義為在波數為k單位平面波激勵下,x位置處的振速或聲壓響應。

2 非相關平面波流激聲輻射有限元法

2.1 非相關平面波疊加等效法

用非相關平面波疊加來合成湍流脈動壓力,是指通過一組空間隨機分布的廣義平面波疊加來等效湍流脈動壓力激勵力,該技術的關鍵是獲得平面波的等效幅值,具體的理論如下:

假設每一組平面波(這里指的廣義平面波)的表達式為

Pxy(x,y,t)=Axy(t)ejkxx+jkyy

(4)

式中:Axy(t)為平面波幅值;kx,ky分別為x,y方向的波數。

對空間任意兩點處的平面波作統計相關運算然后作傅里葉變換,得到空間-頻率譜函數表達式為

SPxyPxy(ξx,ξy,ω)=SAxyAxy(ω)ejkxξx+jkyξy

(5)

式中:ξx,ξy分別為空間兩點沿x軸和y軸的空間距離;SAxyAxy(ω)為平面波的自譜。

當多組形如(4)式的平面波疊加作用時,空間某點(x,y)處的壓力之和為

(6)

由于假設了各組平面波之間非相關,那么滿足當x≠x′,y≠y′時

SAxyAx′y′(ω)=0

(7)

相應的總壓力空間-頻率譜密度函數為

(8)

式(4)～式(8)給出了多組平面波疊加后的總壓力空間-頻率譜,下面要將這種譜與湍流脈動壓力波數-頻率譜建立聯系。

由空間傅里葉變換關系可知,湍流脈動壓力的空間-頻率譜與波數-頻率譜滿足如下關系

SPP(ξx,ξy,ω)=

(9)

式中,ΦPP(kx,ky,ω)可由脈動壓力的波數-頻率譜經驗公式給出。對式(9)的積分進行波數域離散得到

SPP(ξx,ξy,ω)=

(10)

對比式(8)和式(10)可知,當自譜SArsArs(ω)滿足式(11)時,湍流脈動壓力波數-頻率譜可由多組非相關平面波疊加來等效

(11)

2.2 流激聲輻射有限元法流程

由2.1節分析可知合成的湍流脈動壓力表達式為

(12)

式中,φij為均勻分布在[0,2π]上的空間隨機相位,表示平面波是空間非相關的。

(13)

(14)

式中,α1=0.116,α3=0.7分別為沿流向和展向的衰減系數。

表達式(12)將湍流脈動壓力隨機激勵轉化為一種確定性激勵,此時可以通過常規有限元聲振耦合技術來計算流激振動聲輻射,對應的有限元方程為

(15)

式中:ul為結構位移;pl為輻射聲壓;K,H,M分別為剛度矩陣、耦合矩陣和質量矩陣;下標s和f分別為結構和流體。

本文對式(15)的求解過程是在商業有限元軟件Comsol平臺實現的,得到ul和pl后,結構的流激振動自譜以及輻射聲功率譜分別為

Svv=-ω2E(uu*)

(16)

(17)

式中,E為數學期望。

計算流程如圖2所示,首先利用式(12)得到等效平面波的幅值,將多組空間隨機分布的平面波通過接口加載到Comsol平臺建立聲振耦合模型,將各組平面波激勵下的響應進行平均即可得到流激振動與聲輻射。

圖2 計算流程圖

3 算法驗證

由前邊理論可知,湍流脈動壓力是由多組隨機平面波(見式(12))所合成,因此最終的流激振動聲輻射響應由各組平面波激勵下的響應取平均值得到,圖3展示了各組平面波激勵下的響應曲線,圖中黑色實線為各組響應的均值,后續計算分析中均采用此法。

例如，在學習《多位數乘一位數》這部分內容的時候，我為了讓學生能夠熟練地口算各種算式，我把學生分成了男女兩隊，讓學生進行搶答對抗競賽，男女兩隊在活動中，表現得特別踴躍，比賽狀況十分膠著，兩隊都鉚足了勁兒，要比出個上下來，通過組織這樣的對抗競賽，活躍了學生的思維的同時，學生為了取得好成績，都積極地開動腦筋搶答問題，而且在競賽中學生之間互相合作、團結互助，促進了學生之間的友誼，我對學生在競賽中的良好表現做出了表揚，學生學習的勁頭更足了，這樣一來，這次的趣味性教學獲得了圓滿成功。

圖3 多組激勵下的平均響應

為了驗證第2節中算法的正確性,在有限元軟件Comsol中建立如圖4所示的四邊簡支平板的有限元聲振耦合模塊,將多組等效平面波激勵分別加載到彈性板表面,利用圖3中的處理方法得到最終的平均響應。將有限元數值計算結果與流激平板解析值進行對比,來驗證本文算法的正確性。

圖4 有限元模型示意圖

彈性板尺寸為Lx=Ly=0.5 m,厚度h=5 mm,材料為鋼板,來流速度為U∞=5 m/s,后續計算分析均用此組參數。從圖5可看出,本文算法與解析算法有較好的一致性,說明湍流脈動壓力平面波等效合成方法結合有限元計算能夠準確的預報流激聲輻射。圖中低頻時兩者有一定的誤差,可能是因為解析解假設了無限大剛性障,而有限元法只能取有限尺寸的障板。

圖5 算法驗證

4 參數分析

由合成的湍流脈動壓力表達式(12)可知,影響流激振動聲輻射精度的因素包括截斷波數范圍、波數分辨率、平面波組數,另外有限元網格的劃分準則也會影響計算精度,本節就這些問題展開討論。

4.1 截斷波數選取原則

圖6 湍流脈動壓力波數頻率譜

圖7 平板振動傳遞函數波數頻率譜

從圖6可看出,湍流脈動壓力最大值分布在遷移波數附近(白色虛線)。從圖7可看出,平板振動傳遞函數響應最大值主要分布在小于等于彎曲波數的范圍內。從圖8可看出,流激振動響應的最大值也分布在彎曲波數附近。事實上只有當流速非常高時,流體遷移波數才能與結構彎曲波數耦合。因此這里初步將截斷波數上限取定為結構彎曲波數。

圖與乘積波數頻率譜

為了精細確認截斷波數對流激聲輻射的影響,圖9給出了0.5kf和1kf,2kf(kf為彎曲波數)四種截斷波數下的輻射聲功率曲線,可以看出截斷波數取兩倍彎曲波數Kcutoff時,計算結果收斂。圖10給出了1 000 Hz時兩種截斷波數下(0.5kf與4kf)板表面湍流脈動壓力分布,4kf時的脈動壓力場更精細。

圖9 截斷波數的影響

圖10 不同截斷波數對湍流脈動壓力場的影響

4.2 波數分辨率的影響

由式(12)可知,波數分辨率取值會影響湍流脈動壓力譜的精度,原則上分辨率越高(對應波數分辨率數值越小)計算結果越準確。圖11給出了幾種不同波數分辨率取值對流激聲輻射的影響,可以看出,波數分辨率對低頻聲輻射影響不大,對高頻聲輻射有一定的影響。綜合考慮計算量和計算精度,這里取波數分辨率為。

圖11 波數分辨率的影響

4.3 平面波組數的影響

在合成湍流脈動壓力場時,理論上需要無窮多組形如式(12)的平面波疊加,然而在實際仿真時只能取有限組數,圖12給出了不同組數非相關平面波疊加后的平板流激聲輻射對比,可以看出隨著平面波組數的增加,計算結果趨于收斂,平面波取10組與取20組的計算結果基本一致,綜合考慮計算量的因素,這里取10組非相關平面波來合成湍流脈動壓力場。

圖12 平面波組數的影響

4.4 網格密度的影響

在采用有限元法計算聲振耦合問題時,網格密度是影響計算精度的關鍵因素,尤其對于流激聲振耦合問題的網格密度,不僅要考慮彈性結構、聲場因素,還要考慮湍流脈動壓力的空間分布,文獻[9]指出,要獲得收斂的流激振動聲輻射結果,網格劃分密度要遠超常規的一個波長6個節點的準則,因此計算時間成本很大,無法用于分析一般工程問題。

另一方面,由于流激振動聲輻射是隨機問題,所以湍流脈動壓力激勵要考慮空間各點的互譜,從而導致激勵力矩陣維度很大(需要進行Choslesky矩陣分解),進一步影響了計算效率。而本文提出的方法,則是將隨機聲振耦合問題轉化為一般的確定性問題,降低了激勵力矩陣維度,能有效提高計算效率。

圖13比較了傳統Choslesky分解法與本文算法的計算結果,可以看出采用同樣的網格密度,本文算法與解析解有較好的一致性,而采用Choslesky分解法在中高頻有較大的誤差。要指出的是本文算法網格取0.5 mm時對應于一個波長4個節點的劃分原則,即網格取π/2kf時,計算結果足夠收斂。表1給出了不同算法單頻的計算時間,進一步表明了本文算法的計算優勢(所使用計算機內存為32 GB,CPU為酷睿i7 9700)。

表1 不同工況求解時間對比

圖13 網格密度的影響

5 結 論

本文以非相關平面波疊加技術為基礎獲得等效湍流脈動壓力激勵,并結合有限元聲振耦合模塊計算結構流激噪聲,將隨機聲學問題轉化為確定性問題。在深入分析流激噪聲產生機理的基礎上給出了影響該方法精度的準則,并得到如下結論:

(1)本文提出的數值預報方法能夠準確預報結構流激振動噪聲,并能顯著提高計算效率。

(2)截斷波數范圍、波數分辨率、平面波組數等參數是影響湍流脈動壓力激勵力準確性的關鍵因素,需要根據文中的準則進行合理設。

(3)采用非相關平面波疊加法計算流激噪聲時,聲振耦合模型有限元網格密度要比傳統方法稀疏,從而減少了計算負擔。

為了使物理意義清晰明了,本文在進行流程描述和參數分析時均以矩形平板為例,后續研究將探索采用非相關平面波疊加法計算任意復雜結構的流激噪聲。