可壓縮渦流場中空泡運動規律及聲輻射特性研究*

葉 曦 姚熊亮 張阿漫 龐福振

(哈爾濱工程大學船舶工程學院，哈爾濱 150001)

(2012年12月18日收到;2013年2月11日收到修改稿)

1 引言

流體通過有限翼展的翼狀結構，會在其后釋放自由渦[1].由流場中存在的氣核，以及翼狀結構高速運動所形成的低壓區內生成的空泡，被渦流場捕捉后，會發生如延展、撕裂、射流等一系列復雜的運動[2-6].船舶高速航行時，螺旋槳旋轉在船后形成渦流場，而尾跡內存在大量的空泡[7]，同時空泡在坍塌、潰滅過程中會對螺旋槳表面產生剝蝕，并輻射大量噪聲，對于隱身性、安全性要求較高的軍用艦船，不僅會使其易被敵方偵測而暴露位置，且會對艦船的強度及操縱性造成負面影響.因此，研究渦流場中的空泡運動規律及聲輻射特性，有助于提高艦船航行的安全、隱蔽性，具有較高的工程實用價值.

空泡輻射噪聲的求解多是采用邊界積分方程.Kirchhoff邊界積分方程假定流場滿足線性波動方程，廣泛應用于內外聲場的求解[8,9].考慮到運動結構聲場求解的需要，Morgans[10]和Farassat[11]將廣義Kirchhoff方程進行擴展，由包圍結構的運動變形邊界，求得內外場聲輻射.然而，Farassat在求解過程中提出的坐標變換方法，對于某些形狀的邊界面，邊界積分方程中的被積函數是多值的，無法直接進行求解，因此對邊界面的形狀有所限制.

由于流體的可壓縮性，空泡的能量會在運動過程中不斷耗散[12]，以聲能的形式向外傳播.以往所采用的數值模型多是基于流場不可壓縮的假設[13,14]，流場中的氣泡運動及聲能不隨時間衰減，總能量不守恒.如Choi等人在文獻[2—4]中關于輻射聲壓的結果，由拉普拉斯方程對應的邊界積分方程直接求得，未計空泡能量的衰減和聲波傳播的延遲效應.魯傳敬、戚定滿[15,16]基于混合邊界元法，采用不可壓縮邊界元法求解空泡的運動規律，在空泡近場處設置一包圍空泡的固定虛擬面作為輻射聲源，采用定邊界積分方程求解流場中的聲壓分布.然而，為了便于將節點造成的擾動在測點處疊加，虛擬面上網格的劃分取決于測點與空泡的相對位置，同時單元長度與時間步相關;改變測點位置后需對虛擬面重新離散，多測點求解時較為不便，而減小時間步則會增加網格數量，提高計算耗費.Turangan等[17]采用自由拉格朗日法求解沖擊波誘導下的氣泡運動，通過定邊界FH-W方程[18,19]與Kirchhoff方程，同樣借助設置固定虛擬面求解流場中的輻射聲壓.基于固定虛擬面的輻射噪聲求解，實際上將空泡的運動與噪聲計算相分離.一方面虛擬面的布設位置不能距離空泡表面過遠，在使用定邊界的Kirchhoff方程時，虛擬面需包圍所有非線性因素[18]，故也不能距離空泡表面過近;另一方面，當空泡運動過程中的最大尺寸與最小尺寸、不同時刻的空泡位置相差過大，或是多空泡的輻射噪聲計算，使固定虛擬面與空泡表面相距較遠，無法獲得準確結果.可知，固定虛擬面在布設位置及求解工況上有較多限制.

本文基于可壓縮流體力學理論，借鑒DAA法[20-22]思想，采用計及流場可壓縮性的邊界積分方程，求解渦流場中的空泡運動規律.為了克服固定虛擬面的不足，本文將動邊界積分方程直接用于運動中的空泡表面，將空泡表面作為噪聲源，求解其聲輻射特性.采用單元離散映射的方式，可將動邊界積分方程用于任意形狀的邊界面.

2 可壓縮渦流場中的空泡運動基本理論

2.1 可壓縮渦流場中的邊界積分方程

本文計算模型如圖1所示，空泡周圍為可壓縮流體，空泡運動過程絕熱，且不計及粘性的影響，流場有一沿z軸正向的穩定流速uw，考慮球形氣核被渦流捕捉后的空泡運動特性.由Helmholtz分解定理可知，流場中的全速度矢量uf可表示為uf=ub+wf=?φ+wf.其中，ub=?φ由空泡誘導產生，φ為流場速度勢，wf為有旋流場速度.假定流場的擾動速度僅由空泡誘導產生，同時空泡的存在及其運動不影響背景渦流場，而渦流場的存在會對空泡運動產生影響，即空泡與渦流場之間為單向耦合[2].考慮小幅擾動，將渦流場表示為某矢量勢β的旋度wf=?×β，由連續性方程可知

其中，ρ0為靜止流體的密度，ρ′為密度擾動.由于?·(?×β)=0，可知

其中，對于小振幅擾動而言c≈1500m·s-1.利用波動方程的格林函數[23]并結合格林第二定理可得Kirchhoff延遲勢方程

其中，α為立體角;rp為場點位置矢量，rq為源點位置矢量;nq為源點處的單位法向量;τ=t-rpq/c為源點時間，t為場點時間.為了求解空泡的運動規律，將場點p及源點q都布置在空泡表面上.當空泡運動時間較小時，表面某點的物理量僅與周圍極小的球面微元相關[20]，由此可得空泡的前期局部近似

其中，κ為氣泡表面的曲率.隨著空泡的運動，表面上某點的速度勢開始與整個表面相關[22]，此時將(4)式進行泰勒展開，并保留至前2階，可得空泡運動的后期全局近似

其中，?，Θ，?為對應的系數矩陣.將(5)式與(7)式進行匹配，可得可壓縮流場中的邊界積分方程[24,25]

其中，Φ為空泡表面的速度勢矩陣;κ為空泡表面的局部曲率矩陣;Π =(Θ-1?-κ)(δ-Θ-1?Θ-1?)-1，δ為單位矩陣.令聲速無窮大，(8)式即簡化為不可壓縮流場中的邊界積分方程.

圖1 渦流場及空泡模型

2.2 渦流場中空泡運動邊界條件

在Lagrange觀點下，本文中空泡表面滿足的運動學邊界條件為

其中，xb為空泡表面節點的位置矢量.將uf代入Navier-Stokes方程，由于單向耦合的假設，背景渦流場仍滿足Navier-stokes方程，同時忽略密度的小幅度擾動，可得到修正的伯努利方程

其中，P，Pvor分別為流場中的壓力與渦流場誘導壓力.采用高斯渦模型[26]描述渦流場，流場中的渦流誘導壓力為

空泡表面處的壓力分布Pvorb采用SAP(Surface Average Pressure)模型[2]與NP(Node Pressure)模型描述.其中，SAP模型可用于渦流場中球狀空泡運動規律的求解，其對應的Pvorb為空泡表面壓力的面積加權平均，即表面流場一側各處的壓力相同，當氣泡半徑較小，空泡表面壓力的梯度相比空泡表面張力可以忽略不計時，采用SAP模型計算所得運動規律與實際情況極為相近;而NP模型可用于非球狀空泡的求解，對應的Pvorb表達式如(11)式所示.

假定渦流場僅存在于軸線方向，且在無窮遠處空泡對流場的擾動衰減為0，由此將(10)式簡化，并采用Lagrangian形式表示，可得空泡表面動力學邊界條件

其中，P0，V0為空泡內部的初始壓力與初始體積，Pc為飽和蒸汽壓，γ為氣體比熱，ι為表面張力系數.注意到流場中靜壓P0=0.5ρ0u2wσ∞+Pc，其中σ∞為空泡數.聯合求解(8)，(9)，(12)式以及(13)式即可得獲得可壓縮渦流場中空泡的運動規律.

3 空泡運動誘導輻射聲壓的求解

3.1 求解運動邊界聲輻射的邊界積分方程

在空泡運動及其表面上物理量可知的基礎上，求解其在流場中誘導的聲輻射.用于運動可變形邊界聲輻射求解的邊界積分方程為[11]

3.2 動邊界積分方程與渦流場中空泡的耦合

對于渦流場中的空泡而言，其表面本身由流體質點構成，因此表面運動變形的位移、速度與表面處流體的運動位移、速度一致.本文中，流場有一向前的流速uw，而空泡隨著流體一同向前運動，即空泡與流體仍保持相對靜止，故此處馬赫數M 中的邊界面運動速度取為v=?φ，φ取為空泡表面速度勢，同時邊界面上各點的位置矢量取為x=xb.

Farassat[11]的求解方式對于某些邊界面形狀可能會使Kirchhoff方程的被積函數為多值.為了將(14)式所示的動邊界積分方程直接應用于空泡運動誘導的聲輻射求解，本文采用表面離散映射的方式求解(14)式.假定流場中觀測點在某觀測時刻的一部分擾動是由空泡表面的節點n在空泡的運動時刻τ?造成，稱節點n為測點在此觀測時刻的擾動節點，其周邊單元離散及映射方式如圖2所示.將各三角形單元映射到ξ-η坐標系中，A，B，C三個頂點坐標分別為 (1，0)，(0，1)，(0，0)，單元上任意位置處的物理量λ可用線性插值表示：

圖2 擾動節點周邊的離散及單元映射

4 數值驗證

單個空泡脈動聲場的解析表達式如下所示[28]，將空泡視作小脈動單極子聲源，獲得其遠場聲輻射：

其中，rm為測點位置與空泡中心間的距離，由于空泡隨流場移動，rm將隨時間改變;rm/c表示延遲時間;R=R(τ?)為空泡半徑，由球狀空泡的Gilmore方程[29]求得.

為驗證本文方法的正確性，由可壓縮流場中的邊界積分方程(8)求解空泡運動，由動邊界積分方程計算空泡運動過程中的輻射噪聲，并將數值結果與解析結果進行比較.計算參數為：R0=1 mm，ac=5.6 mm，uw=10 m/s，Γ =0.28 m2/s，空泡數σ∞=3，空泡初始位置為 (0，0，0)，流場中測點位置(15 mm，0，0).為便于和解析解進行比較，數值計算中表面壓力采用SAP模型.對比結果如圖3和圖4所示，可知本文方法所得結果與解析解符合較好，可用于可壓縮渦流場中空泡運動規律及其聲輻射的求解.由圖3可知，由于計及流場的可壓縮性，空泡的半徑隨運動周期發生衰減.由圖4可知，空泡隨著流體向前運動，與測點間的距離增加，因此聲壓幅值不斷減小.考慮可壓縮性的影響時，由于空泡運動幅度減弱，其誘導聲壓隨時間的衰減與不可壓縮時相比更為明顯.空泡脈動的周期約為0.36 ms，可知其脈動所形成的聲波波長遠大于空泡脈動過程中的最大尺寸.

為進一步驗證本文方法可用于渦流場中非球狀空泡運動的求解.采用文獻[5]中的實驗參數，由可壓縮邊界積分方程計算梢渦場中的空泡運動歷程.其中，流場相關參數Γ=0.2123 m2/s，ac=4.51 mm，uw=10 m/s，σ∞=1.72，空泡初始半徑R0=750μm，保持空泡內部壓力不變，而空泡表面流場一側的壓力分布采用NP模型.如圖6(a)所示，為可壓縮邊界積分方程所得結果，6(b)為對應的試驗結果.圖中顏色表示速度勢的分布，紅色表示高速度勢，藍色表示低速度勢.試驗中空泡變形后的最大長度為20 mm，數值計算中空泡的最大長度為20.19 mm.由此可知，對于梢渦流場中的非球狀空泡，本文數值方法仍能獲得正確的計算結果.

圖3 空泡半徑對比曲線

圖4 空泡脈動輻射聲壓對比曲線

圖5 渦流場中的空泡形態(單位：ms)

圖6 數值結果與試驗結果對比 (a)數值結果;(b)實驗結果(文獻[5])

如圖5所示，初始為球形的靜止空泡，由于渦流場中誘導壓力的作用，球形空泡兩端處流場的壓力低于中部的壓力，空泡開始伸長.t=3.74 ms后在空泡兩端發生頸縮，形成球形腔體.t=5.20 ms時，空泡頸縮處與渦軸線間的距離極小，將發生分裂，形成三個子空泡.

5 數值計算及結果討論

5.1 可壓縮渦流場中的空泡運動規律及聲輻射特性

如圖7所示為渦流場中，初始靜止球形空泡延展、撕裂、射流的完整過程.空泡初始半徑為R0=2 mm，流場參數為Γ =0.3 m2/s，ac=5 mm，uw=10 m/s，σ∞=2.8.空泡運動規律以及輻射聲壓的求解截止于射流沖擊之前.

運動初期，空泡在渦流場誘導壓力的作用下發生延展，如圖7(a)—(c)所示.隨后，空泡在延展的同時，由于中部流場壓力較大，故此處發生頸縮，逐漸形成兩個連通的球形腔，如圖7(d)—(f).在t=0.89 ms時，空泡發生撕裂，形成兩個子空泡.撕裂后的空泡運動如圖7(g)—(j)所示，子空泡間的距離不斷增加，同時分別自撕裂處朝向子空泡內部形成射流.

如圖8和圖9所示為采用邊界積分方程由SAP模型所得球狀空泡與NP模型所得非球狀空泡體積歷程曲線，以及在(0，8m m，0)處的輻射聲壓曲線.運動初期兩者體積曲線基本重合;自t=0.3 ms開始，隨著非球狀空泡的不斷延展，兩者體積的差異開始增大，非球狀空泡的體積大于球狀空泡.在空泡達到最小體積附近，非球狀空泡發生撕裂，子空泡無明顯延展，總體積與球狀空泡相近.

圖7 渦流場中空泡延展、頸縮、撕裂的變形歷程(單位：ms)

圖8 SAP模型與NP模型體積對比曲線

如圖9所示，當空泡運動達到最大體積時，其誘導的輻射聲壓達到波谷;而隨著空泡的潰滅，輻射聲壓在空泡達到最小體積時取到波峰.撕裂前，非球狀空泡與球狀空泡輻射聲壓特性相近，但球狀空泡的聲壓波動周期短于非球狀空泡.非球狀空泡撕裂時，對應的輻射聲壓發生突變.除空泡撕裂對流場造成的擾動外，空泡撕裂瞬時的數值處理對聲壓突變也有所貢獻[3].撕裂后，空泡脈動所產生的聲波波長遠大于空泡尺寸，因此子空泡間的聲散射可忽略不計[30].由于空泡隨流場向前運動，空泡與測點間的距離不斷增大，且考慮到流場可壓縮性對空泡能量的耗散，除撕裂造成的峰值外，輻射聲壓的峰值按周期遞減.

圖9 SAP模型與NP模型輻射聲壓對比曲線

圖10為空泡運動過程中周向聲壓分布情況，測點布置在半徑為8 mm且平行于來流方向圓上，隨著空泡一同向前運動.初始時刻，由于空泡以延展為主，兩端的輻射聲壓幅值較大，如圖10(a)所示.隨著空泡中部形成頸縮，對流場的擾動增大，導致中部的聲壓幅值增大，如圖10(b)所示.空泡的延展、頸縮使得空泡兩端與測點的距離減小，而中部與測點的距離增大，因此撕裂前，空泡兩端的聲壓幅值再次增大，如圖10(c)所示.空泡頸縮至極限后，撕裂成兩個子空泡，分別朝向位于最前和最后處的測點運動，并形成射流，此時空泡與兩端測點距離進一步減小，故空泡輻射聲壓的極大值仍位于兩端，如圖10(d)所示.由圖10可知，空泡運動過程中的周向輻射聲壓幅值的差距并不大，這是由于空泡延展頸縮直至撕裂形成子空泡的過程中，都類似于單個小球源或是同相小球源的脈動，小球源之間的距離遠小于聲波波長，故聲壓指向性較弱[30].

圖10 空泡輻射聲壓指向性(單位：ms)

5.2 渦流場參數對空泡運動規律及聲輻射特性的影響

渦流場的渦通量、空泡數、以及流速等會對空泡運動及其聲輻射特性產生影響.以下分析中空泡初始半徑R0=1 mm，流場中uw=10 m/s，ac=4.5 mm.由(11)式可知，渦通量的增大，空泡數以及流速的減小，會使渦流場誘導壓力降低.如圖11所示，令渦通量Γ=0.28 m2/s保持不變，隨著空泡數的增大，流場壓力增加，空泡運動過程中的最大等效半徑(與空泡體積相等的球體半徑，已用初始半徑無量綱化)逐漸減小.空泡數較小時，非球狀空泡的最大半徑高于球狀空泡.高空泡數時，由于空泡半徑較小，SAP模型采用的平均壓力與非球狀空泡各位置處的壓力差距較小，因此兩者的最大等效半徑基本相同.

如圖12所示為空泡數σ∞=2.5不變，渦通量增加時，空泡撕裂時刻的總長度變化曲線.由(11)式可知，隨著渦通量的增大，流場同一位置處沿空泡延展方向的壓力梯度增大，空泡撕裂時刻的總長度不斷增加.當Γ＞0.34 m2/s時，總長度隨渦通量近似線性增長.以Γ≈0.34 m2/s為界，當Γ＜0.34 m2/s時，空泡撕裂后形成兩個子空泡，而Γ＞0.34 m2/s時，空泡撕裂后將形成三個子空泡，即隨著渦通量的增大，子空泡的數量增多.

圖11 空泡數對空泡最大半徑的影響

圖12 渦通量對最大空泡長度的影響

圖13為Γ=0.35 m2/s時的空泡運動歷程.運動初期空泡在體積膨脹的同時不斷延展，如圖13(a)，(b)所示.隨后空泡中部發生頸縮，同時仍不斷的進行延展，如圖13(c)—(e)所示.隨著延展的進行，空泡發生兩次頸縮，在空泡兩端形成兩個球形腔體，中間則為一梭形腔體，如圖13(f)，(g)所示.t=1.82 ms時，空泡發生撕裂，形成三個子空泡，如圖13(h)所示.撕裂后，子空泡間的距離不斷增大，同時在撕裂處，兩側的球形子空泡向各自內部形成射流，而中部的梭形子空泡形成對射流，如圖13(i)，(j)所示.

如圖14所示為不同渦通量、空泡數下，空泡撕裂瞬時的形態及射流形態.渦通量較小或空泡數較大時，空泡撕裂瞬時兩端球形腔體較為飽滿，腔體的形心偏向于撕裂點.隨著渦通量的增大或是空泡數的減小，撕裂瞬時的空泡被逐漸“拉長”，球形腔體的形心逐漸遠離撕裂點.且渦通量較小時(或空泡數較大時)，子空泡射流前端呈明顯的水滴狀.圖15 為測點 (0，8 mm，0)處，Γ =0.25,0.28,0.32 m2/s時的空泡輻射聲壓比較曲線.隨著渦通量的增大，空泡脈動周期與輻射聲壓的波動周期延長，撕裂時刻滯后，且空泡撕裂后直至射流沖擊的時間也逐漸減少.同時，空泡運動初期輻射聲壓值增大，而撕裂前空泡輻射聲壓的峰值逐漸滯后并減小.

圖13 渦流場中空泡撕裂為多個子空泡(單位：ms)

圖 14 渦流場參數對空泡形態的影響 (單位：ms)(a)Γ =0.25 m2/s，σ∞=2.5;(b)Γ =0.28 m2/s，σ∞=2.5;(c)Γ =0.32 m2/s，σ∞ =2.5;(d)Γ =0.28 m2/s，σ∞ =2.1;(e)Γ =0.28 m2/s，σ∞ =3

6 結論

本文基于可壓縮流場中的邊界元法，采用SAP模型與NP模型，求解了可壓縮渦流場中的空泡運動規律;通過單元離散及映射的方式，將動邊界積分方程直接用于空泡表面，求得流場中的輻射聲壓分布;分析比較了渦流場參數對空泡運動形態及聲輻射的影響.得出以下結論：

圖15 渦通量對空泡輻射聲壓的影響

1.計及流場可壓縮性，空泡運動會造成流場中的聲輻射;空泡脈動幅度隨時間不斷減小，流場中某點的輻射聲壓幅值也隨之減小.

2.空泡在渦流場中發生延展，同時形成頸縮，最終分裂成若干個子空泡，并分別形成朝向子空泡內部的射流.隨著渦通量、空泡數的變化，空泡周圍流場中的壓力會隨之改變.當流場中的壓力減小(空泡數減小或渦通量增大)時，空泡運動過程中的最大半徑及撕裂前的最大長度逐漸增加，且當流場中壓力較小時，空泡撕裂時形成的子空泡增多.

3.空泡誘導的輻射噪聲分別在最大體積時取到峰值，而在最小體積時取到谷值，同時非球狀空泡撕裂會使輻射聲壓產生突變，形成極大峰值;球形空泡模型的輻射聲壓脈動周期短于非球狀空泡模型.渦通量越大，運動初期的空泡輻射聲壓越高，且隨著渦通量的增大，空泡脈動周期也隨之增加，輻射聲壓峰值逐漸滯后并減小.撕裂之前，空泡在流場中的脈動類似于脈動小球源，而撕裂后則類似于多個同相脈動小球源，因此輻射聲壓的指向性較弱.

[1]HsiaoC T，Pauley L L 1999 J.Fluids Eng.121 198

[2]ChoiJ K，Chahine G L 2002 InternationalAssociation for Boundary Element Method Austin，TX，USA，2002，May 28—30，2002 p1

[3]ChoiJ K，Chahine G L 2003 The 8thInternationalConference on NumericalShip Hydrodynamics Busan，Korea，September 22—25 2003

[4]HsialC T，Pauley L L 2003 J.Fluids Eng.125 53

[5]Rebow M，ChoiJ，ChoiJ K，Chahine G L，CeccioS L 2004 11thInternationalSynposium on Flow Visualization Indiana，USA，August 9—12 2004 p1

[6]NiB Y，Zhang A M 2012 Appl.Math Mech.33 701

[7]Carrica M 1999 Int.J Multiphas.Flow 25 257

[8]Pierce A D 1981 Acoustic:An Introduction toIts PhysicalPrinciples and Applications(New York：McGraw-Hill)p180

[9]Hawkings D L 1979 Mechanics of sound generation in flows Goettingen，West Germany，August 28—31 1979 p294

[10]Morgans W R 1930 Philos.Mag.9 141

[11]Farassat F 1988 J.Sound Vib.123 451

[12]Wang S P 2011 Ph.D.Dissertation(Harbin：Harbin Engineering University)(in Chinese)[王詩平2011博士學位論文(哈爾濱：哈爾濱工程大學)]

[13]Zhang A M，YaoX L 2008 Acta Phys.Sin.57 339(in Chinese)[張阿漫，姚熊亮2008物理學報57 339]

[14]Zhang A M，YaoX L，LiJ 2008 Acta Phys.Sin.57 1672(in Chinese)[張阿漫，姚熊亮，李佳2008物理學報57 1672]

[15]QiD M，Lu C J 2001 J.Hydrodyn 16 9(in Chinese)[戚定滿，魯傳敬2001水動力學研究與進展16 9]

[16]QiD M 1999 Ph.D.Dissertation(Shanghai：ShanghaiJiaotong University)(in Chinese)[戚定滿1999博士學位論文(上海：上海交通大學)]

[17]Jamaluddin A R，Turangan C K 2011 J.Fluid Mech.677 305

[18]FrancescantonioDi1997 J.Sound Vib.202 491

[19]Farassat F 1983 Vertica 7 309

[20]Geers T L 1978 J.Acoust.Soc.Am.64 1500

[21]Geers T L 1971 J.Acoust.Soc.Am.49 1505

[22]Geers T L 1980 J.Acoust.Soc.Am.173 1152

[23]Liang K M 2010 Methods of MathematicalPhysics(Beijing：Higher Education Press)(in Chinese)[梁昆淼2010數學物理方法(北京：高等教育出版社)]

[24]Wang S P，Sun S L，Zhang A M 2012 Chinese Journalof Theoreticaland Applied Mechanics 44 513(in Chinese)[王詩平，孫士麗，張阿漫2012力學學報44 513]

[25]Zhang A M，Wang S P，Wu G X 2013 Eng.Anal.Bound.Elem.DOI：10.1016/j.enganabound.2013.04.013

[26]Saffman P G 1992 Vortex Dynamics(Cambridge：Cambridge University Press)

[27]Best J P 1993 J.Fluid Mech.251 79

[28]Rose D 1976 Mechanices of Underwater Noise(Pergamon)p62

[29]Gilmore F G 1952 HydroLab California Institute TechnicalReport 26 117

[30]Du G H，Zhu Z M 2001 Acoustics Foundation(Nanjing：Nanjing University)(in Chinese)[杜功煥，朱哲民，2001聲學基礎(南京：南京大學出版社)]