周子棋,孫一頡,韓沛東,孫中國(guó),席光

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,710049,西安)

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

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

1 數(shù)值方法

1.1 控制方程

MPS方法中,不可壓縮流動(dòng)下的質(zhì)量守恒方程與動(dòng)量守恒方程為

(1)

(2)

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

1.2 核函數(shù)和離散算子

核函數(shù)描述了粒子之間的作用關(guān)系,需要滿足緊支性條件[11],本文采用的經(jīng)典核函數(shù)及離散算子如下

(3)

(4)

(5)

式中:w(r)為核函數(shù);r=|ri-rj|為i、j粒子之間距離;re一般取2.1l0;φ是任意標(biāo)量函數(shù);d為維度;n0是粒子數(shù)密度常數(shù)。在忽略連續(xù)域內(nèi)核函數(shù)截?cái)嗾`差的情況下,拉普拉斯算子可簡(jiǎn)化為

(6)

(7)

2 聲學(xué)計(jì)算模型

本文構(gòu)建的聲學(xué)計(jì)算模型稱為MPS-WP模型,包括聲學(xué)控制方程的粒子離散表達(dá)方程組、經(jīng)典四階Runge-Kutta時(shí)間積分法[12]和典型邊界條件模型。MPS-WP模型應(yīng)用了MPS方法中的離散算子和粒子表達(dá)形式,可在MPS方法求解得到的流場(chǎng)信息的基礎(chǔ)上計(jì)算聲場(chǎng)信息,擴(kuò)大了MPS方法的應(yīng)用范圍。

2.1 聲學(xué)控制方程及方程離散

聲波傳播的介質(zhì)多為空氣或水等黏性及導(dǎo)熱系數(shù)均較小的流體,其擾動(dòng)的空間梯度相比擾動(dòng)本身較小,因此黏性和熱傳導(dǎo)的影響可以忽略不計(jì),擾動(dòng)的傳播可以通過(guò)求解歐拉方程、連續(xù)方程和聲速方程來(lái)確定。其中,連續(xù)性方程為

(8)

歐拉(動(dòng)量)方程為

(9)

聲速方程為

(10)

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

(11)

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

(12)

方程組(11)可以改寫為

(13)

(14)

(15)

方程(13)～(15)定義了流聲分離假設(shè)下聲學(xué)物理量在聲傳播過(guò)程中的隨體導(dǎo)數(shù)表達(dá)式。本文計(jì)算環(huán)境為不可壓縮靜止無(wú)黏流場(chǎng),上式可簡(jiǎn)化為

(16)

(17)

(18)

方程離散后可得到3個(gè)聲學(xué)物理量的導(dǎo)數(shù)表達(dá)式,對(duì)于聲波來(lái)說(shuō),通常為縱波,其擾動(dòng)速度在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在不可壓縮靜止無(wú)黏流場(chǎng)下的表達(dá)形式。如需計(jì)算橫波,即振動(dòng)方向與傳播方向垂直的波,在二維條件下則僅需要一個(gè)擾動(dòng)速度分量,其具體表達(dá)式與式(19)～(22)類似,本文不再贅述。

2.2 時(shí)間積分方法

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

(23)

式中:f(tn,yn)為物理量yn在tn時(shí)刻的導(dǎo)數(shù);h為虛擬時(shí)間步長(zhǎng)。

2.3 邊界條件模擬

粒子方法常采用虛擬粒子來(lái)解決支持域截?cái)鄦?wèn)題,在聲學(xué)計(jì)算中,虛擬粒子也可用于表征聲學(xué)邊界條件。對(duì)于聲學(xué)硬邊界,即在壁面處滿足條件

(24)

采用虛擬粒子在邊界處鏡像內(nèi)部粒子得到

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

(25)

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

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

對(duì)于邊界條件粒子,滿足如下方程

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

(26)

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

(27)

2.4 計(jì)算參數(shù)選取

模型計(jì)算精度受粒子間距l(xiāng)0和CFL數(shù)的影響。定義CFL數(shù)為

(28)

式中:c0為聲波傳播的速度;dt為實(shí)際時(shí)間步長(zhǎng);l0為粒子間距。

針對(duì)100 m×100 m方形計(jì)算域中的高斯脈沖傳播過(guò)程模擬,當(dāng)l0=0.5 m時(shí),CFL數(shù)與聲壓的相對(duì)均方根誤差關(guān)系如圖2所示。當(dāng)CFL數(shù)小于0.68時(shí),不同時(shí)刻的計(jì)算誤差均較小。

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

當(dāng)CFL數(shù)為0.17時(shí),粒子間距與聲壓的相對(duì)均方根誤差關(guān)系如圖3所示。當(dāng)l0小于0.5 m時(shí),不同時(shí)刻的計(jì)算誤差均較小。在本文的計(jì)算中,取CFL數(shù)為0.136,l0=0.25 m。

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

3 結(jié)果與分析

3.1 二維高斯脈沖傳播

參考NASA標(biāo)準(zhǔn)算例[14]建立二維高斯脈沖計(jì)算。初始時(shí)刻聲壓分布為

(29)

任意時(shí)刻的聲壓解析解可以表示為

(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方向提取得到的聲壓分布以及與解析解的對(duì)比,結(jié)果表明MPS-WP的計(jì)算結(jié)果與解析解吻合較好。

圖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 聲學(xué)硬邊界條件下高斯脈沖傳播的聲壓分布Fig.6 Acoustic pressure distributions of Gaussian pulse propagation under acoustic hard boundary condition

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

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

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

聲場(chǎng)計(jì)算中,通常采用吸收邊界來(lái)減小反射波,以消除人工截?cái)噙吔鐚?duì)于計(jì)算域的影響。圖8為在Mur吸收邊界條件下的高斯脈沖傳播。在0.15～0.20 s時(shí),聲波在邊界處完全消失,無(wú)明顯反射波存在。

(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時(shí)MPS-WP在Mur吸收邊界條件下的計(jì)算結(jié)果和解析解吻合良好。

圖9 Mur吸收邊界條件下聲壓的MPS-WP計(jì)算解與解析解的對(duì)比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為不同時(shí)間點(diǎn)聲壓分布的情況,可知脈沖在計(jì)算域內(nèi)的傳播情況較好。在t=0.05 s時(shí)兩脈沖邊緣接觸,出現(xiàn)了聲壓幅值的疊加;在t=0.10 s時(shí),交點(diǎn)位置因?yàn)槁晧函B加達(dá)到最高值,脈沖向前傳播后抵達(dá)吸收邊界。

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

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

3.5 流動(dòng)條件下的聲波傳播

脈沖聲源在運(yùn)動(dòng)流體中傳播存在多普勒效應(yīng)。對(duì)于均勻流動(dòng)下的聲傳播問(wèn)題,有

p0=0;·v0=0

(31)

圖12 均勻流動(dòng)下聲場(chǎng)傳播計(jì)算域布置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 均勻流動(dòng)下多脈沖傳播的聲壓分布Fig.13 Acoustic pressure distributions in propagation of multiple Gauss pulses in uniform flow

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

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

4 結(jié) 論

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