基于移動粒子半隱式方法波傳播模型的聲傳播數值求解方法研究

周子棋,孫一頡,韓沛東,孫中國,席光

(西安交通大學能源與動力工程學院,710049,西安)

聲傳播可用波動方程[1]來表述,求解該方程的數值方法包括有限差分法、有限元法、譜元法和邊界元法等[2]。這些方法須基于網格求解,對網格數量和質量較為敏感。無網格法消除了傳統數值計算對于網格的依賴性,典型的無網格方法有無單元伽遼金法、雜交邊界點法、徑向基函數法、光滑粒子動力學方法、移動粒子半隱式法等[3],其中徑向基函數配點法求解橢圓型、拋物型和雙曲型微分方程獲得了廣泛應用[4]。該類方法具有不用網格剖分、無需積分運算的特點,但其多基于歐拉描述,不易于應用到計算域隨時間變化的非定常問題,因此采用無網格粒子方法建立聲傳播的計算模型較為可行。Hahn對光滑粒子動力學方法計算聲傳播問題時的計算參數和邊界條件設置進行了討論[5]。魏建國等對比了光滑粒子動力學與時域有限差分方法的計算結果,驗證了其正確性[6]。張詠鷗等在一系列論文中逐步完善了聲傳播模型在光滑粒子動力學中的表達形式和相關邊界條件[7-9]。

移動粒子半隱式法(MPS)[10]是一種較新的無網格粒子方法,適合用于求解帶有大變形、自由面、移動邊界等復雜非定常問題。在MPS方法中,聲傳播模型尚未提出,為了建立基于MPS的聲波傳播模擬方法,本文從流聲分離假設[1]出發,推導了聲傳播控制方程,實現了典型邊界條件,研究了CFL數及粒子間距等關鍵參數對計算精度的影響,并通過數值計算算例進行了驗證。

1 數值方法

1.1 控制方程

MPS方法中,不可壓縮流動下的質量守恒方程與動量守恒方程為

(1)

(2)

式中:ρ為流體密度;u為流體速度;p為壓力;μ為動力黏度;f為體積力。

1.2 核函數和離散算子

核函數描述了粒子之間的作用關系,需要滿足緊支性條件[11],本文采用的經典核函數及離散算子如下

(3)

(4)

(5)

式中:w(r)為核函數;r=|ri-rj|為i、j粒子之間距離;re一般取2.1l0;φ是任意標量函數;d為維度;n0是粒子數密度常數。在忽略連續域內核函數截斷誤差的情況下,拉普拉斯算子可簡化為

(6)

(7)

2 聲學計算模型

本文構建的聲學計算模型稱為MPS-WP模型,包括聲學控制方程的粒子離散表達方程組、經典四階Runge-Kutta時間積分法[12]和典型邊界條件模型。MPS-WP模型應用了MPS方法中的離散算子和粒子表達形式,可在MPS方法求解得到的流場信息的基礎上計算聲場信息,擴大了MPS方法的應用范圍。

2.1 聲學控制方程及方程離散

聲波傳播的介質多為空氣或水等黏性及導熱系數均較小的流體,其擾動的空間梯度相比擾動本身較小,因此黏性和熱傳導的影響可以忽略不計,擾動的傳播可以通過求解歐拉方程、連續方程和聲速方程來確定。其中,連續性方程為

(8)

歐拉(動量)方程為

(9)

聲速方程為

(10)

假設流體為單一組分且熱力學平衡,對于定常流,在沒有外力或質量增加的條件下,方程(8)～(10)可以改寫為方程組

(11)

方程組(11)描述了某一擾動在流體中的傳播過程。聲場的基本物理量可表示為流場物理量上疊加的微小擾動,其關系表示為

(12)

方程組(11)可以改寫為

(13)

(14)

(15)

方程(13)～(15)定義了流聲分離假設下聲學物理量在聲傳播過程中的隨體導數表達式。本文計算環境為不可壓縮靜止無黏流場,上式可簡化為

(16)

(17)

(18)

方程離散后可得到3個聲學物理量的導數表達式,對于聲波來說,通常為縱波,其擾動速度在x和y方向存在分量δu和δv。式(16)～(18)可寫為

ri)w(|rj-ri|)

(19)

xi)w(|rj-ri|)

(20)

yi)w(|rj-ri|)

(21)

ri)w(|rj-ri|)

(22)

即MPS-WP在不可壓縮靜止無黏流場下的表達形式。如需計算橫波,即振動方向與傳播方向垂直的波,在二維條件下則僅需要一個擾動速度分量,其具體表達式與式(19)～(22)類似,本文不再贅述。

2.2 時間積分方法

時間積分采用經典四階Runge-Kutta法,替代MPS常用的一階歐拉法yn+1=yn+hf(tn,yn),以保證其在聲波傳播計算中的高精度及收斂性。四階Runge-Kutta時間積分方程為

(23)

式中:f(tn,yn)為物理量yn在tn時刻的導數;h為虛擬時間步長。

2.3 邊界條件模擬

粒子方法常采用虛擬粒子來解決支持域截斷問題,在聲學計算中,虛擬粒子也可用于表征聲學邊界條件。對于聲學硬邊界,即在壁面處滿足條件

(24)

采用虛擬粒子在邊界處鏡像內部粒子得到

δp(i+1)=δp(i);δv(i+1)=0

(25)

對于吸收邊界,本文采用在FDTD中廣泛使用的Mur邊界[13]進行粒子離散表達,如圖1所示。

圖1 右側邊界粒子布置示意圖Fig.1 Boundary particles arrangement on right side

對于邊界條件粒子,滿足如下方程

(fn(i,j)-fn-1(i+1,j))

(26)

(fn(i,j)-fn-1(i+2,j))

(27)

2.4 計算參數選取

模型計算精度受粒子間距l0和CFL數的影響。定義CFL數為

(28)

式中:c0為聲波傳播的速度;dt為實際時間步長;l0為粒子間距。

針對100 m×100 m方形計算域中的高斯脈沖傳播過程模擬,當l0=0.5 m時,CFL數與聲壓的相對均方根誤差關系如圖2所示。當CFL數小于0.68時,不同時刻的計算誤差均較小。

圖2 不同CFL數下的相對均方根誤差Fig.2 Relative root mean square error at different CFL numbers

當CFL數為0.17時,粒子間距與聲壓的相對均方根誤差關系如圖3所示。當l0小于0.5 m時,不同時刻的計算誤差均較小。在本文的計算中,取CFL數為0.136,l0=0.25 m。

圖3 不同粒子間距下的相對均方根誤差Fig.3 Relative root mean square error at different particle spacings

3 結果與分析

3.1 二維高斯脈沖傳播

參考NASA標準算例[14]建立二維高斯脈沖計算。初始時刻聲壓分布為

(29)

任意時刻的聲壓解析解可以表示為

(30)

(a)t=0 s

(b)t=0.05 s

(c)t=0.10 s圖4 二維高斯脈沖傳播的聲壓分布Fig.4 Acoustic pressure distributions in two-dimensional Gauss pulse propagation

圖5為采用MPS-WP在中心位置沿x方向提取得到的聲壓分布以及與解析解的對比,結果表明MPS-WP的計算結果與解析解吻合較好。

圖5 中心位置處沿x方向的聲壓分布Fig.5 Acoustic pressure distributions at center along x-direction

3.2 硬邊界條件下的高斯脈沖傳播

(a)t=0.1 s

(b)t=0.15 s

(c)t=0.20 s圖6 聲學硬邊界條件下高斯脈沖傳播的聲壓分布Fig.6 Acoustic pressure distributions of Gaussian pulse propagation under acoustic hard boundary condition

聲學硬邊界表現為聲波在壁面處的全反射?；谑?25)施加了聲學硬邊界。高斯脈沖的傳播情況如圖6所示。由圖6可知,在0.15～0.20 s間,脈沖與壁面接觸產生反射和疊加。y=0處的聲壓分布如圖7所示,在聲學硬邊界條件下,計算結果與解析解吻合良好。

圖7 聲學硬邊界條件下聲壓的MPS-WP計算解與解析解的對比Fig.7 Acoustic pressure comparison under hard boundary condition between the MPS-WP solution and the analytical solution

3.3 吸收邊界條件下的高斯脈沖傳播

聲場計算中,通常采用吸收邊界來減小反射波,以消除人工截斷邊界對于計算域的影響。圖8為在Mur吸收邊界條件下的高斯脈沖傳播。在0.15～0.20 s時,聲波在邊界處完全消失,無明顯反射波存在。

(a)t=0.1 s

(b)t=0.15 s

(c)t=0.20 s圖8 Mur吸收邊界條件下高斯脈沖傳播的聲壓分布Fig.8 Acoustic pressure distributions in propagation of Gauss pulse under Mur absorbing boundary condition

同樣地,提取在y=0處的聲壓分布(見圖9),可以看到0.15 s時MPS-WP在Mur吸收邊界條件下的計算結果和解析解吻合良好。

圖9 Mur吸收邊界條件下聲壓的MPS-WP計算解與解析解的對比Fig.9 Acoustic pressure comparison under Mur absorbing boundary condition between the MPS-WP solution and the analytical solution

3.4 兩高斯脈沖疊加

(a)t=0 s

(b)t=0.05 s

(c)t=0.10 s圖10 兩高斯脈沖傳播的聲壓分布Fig.10 Acoustic pressure distributions in propagation of two Gauss pulses

圖10為不同時間點聲壓分布的情況,可知脈沖在計算域內的傳播情況較好。在t=0.05 s時兩脈沖邊緣接觸,出現了聲壓幅值的疊加;在t=0.10 s時,交點位置因為聲壓疊加達到最高值,脈沖向前傳播后抵達吸收邊界。

同樣地,采用MPS-WP提取y=0位置的聲壓分布與解析解對比,如圖11所示,采用MPS-WP計算的聲壓分布與解析解吻合較好。

圖11 兩高斯脈沖傳播聲壓的MPS-WP計算解與解析解的對比Fig.11 Acoustic pressure comparison in propagation of two Gauss pulses between the MPS-WP solution and the analytical solution

3.5 流動條件下的聲波傳播

脈沖聲源在運動流體中傳播存在多普勒效應。對于均勻流動下的聲傳播問題,有

p0=0;·v0=0

(31)

圖12 均勻流動下聲場傳播計算域布置Fig.12 Computational domain arrangement of acoustic wave propagation in uniform flow

(a)t=0 s

(b)t=0.01 s

(c)t=0.02 s

(d)t=0.03 s圖13 均勻流動下多脈沖傳播的聲壓分布Fig.13 Acoustic pressure distributions in propagation of multiple Gauss pulses in uniform flow

則依然滿足方程(13)～(15)、(16)～(18)的簡化過程,故可以將該模型拓展到均勻流動的聲傳播計算中。算例中高斯脈沖位于方腔中心,方腔上下邊界為Mur吸收邊界,左右分別為入口和出口邊界;對于流動,進出口分別為均勻速度進口以及出口邊界,兩側壁面為完全滑移邊界,壁面對流動不產生影響,則整個流場的均勻性保持不變。在均勻流動中速度Vr=100 m/s,高斯脈沖的產生頻率為200 Hz,其余參數與上文設置相同,如圖12所示。計算結果如圖13所示,隨著時間推進,聲壓在脈沖源的均勻來流方向疊加,脈沖在下游方向的傳播距離較上游方向遠。t=0.02s時多次脈沖疊加的結果如圖14所示,計算結果與解析解一致。

圖14 均勻來流條件下聲壓的MPS-WP計算解與解析解的對比Fig.14 Comparison of acoustic pressure distributions in uniform flow between the MPS-WP solution and the analytical solution

4 結 論

本文基于無網格MPS粒子方法,在拉格朗日框架下建立了用于計算聲傳播問題的計算模型MPS-WP,實現了硬邊界和Mur吸收邊界兩類邊界條件的離散表達。在本文的算例中,研究了CFL數和粒子間距對計算結果的影響并確定了適用范圍,模擬了二維高斯脈沖在硬邊界和Mur吸收邊界下的傳播過程,驗證了模型和邊界條件的正確性。本文還對靜場中脈沖疊加及均勻來流條件下的連續脈沖傳播過程進行了數值模擬和研究,計算結果與解析解吻合良好,為在復雜流場計算聲傳播問題提供了一種有效途徑。