宣領(lǐng)寬, 龔京風(fēng), 張文平, 袁興軍

(1.中國(guó)艦船研究設(shè)計(jì)中心, 武漢 430064; 2.哈爾濱工程大學(xué) 動(dòng)力與能源工程學(xué)院, 哈爾濱 150001)

求解聲波動(dòng)方程的隱格式有限體積法

宣領(lǐng)寬1, 龔京風(fēng)1, 張文平2, 袁興軍1

(1.中國(guó)艦船研究設(shè)計(jì)中心, 武漢 430064; 2.哈爾濱工程大學(xué) 動(dòng)力與能源工程學(xué)院, 哈爾濱 150001)

摘要:為研究聲傳播問題,提出一種聲波動(dòng)方程的隱格式有限體積法,該方法將格點(diǎn)型有限體積法與Newmark格式相結(jié)合．模擬平面波的傳播過程,對(duì)比分析隱格式有限體積法和文獻(xiàn)中顯格式有限體積法的精度、穩(wěn)定性及計(jì)算消耗等方面的性能．?dāng)?shù)值結(jié)果表明:當(dāng)λ/Δx≥10時(shí),兩種算法均能得到滿足精度要求的解;采用無條件穩(wěn)定的隱格式算法,當(dāng)滿足ω0Δt≤0.3時(shí),預(yù)測(cè)聲壓的相對(duì)峰值誤差<1%;當(dāng)采用相同時(shí)間、空間步長(zhǎng)時(shí),隱格式算法精度高于顯格式算法;隱格式算法對(duì)吸收邊界的處理精度高于顯格式算法,但對(duì)全反射邊界的處理精度低于顯格式算法;兩種算法內(nèi)存消耗比較接近,顯格式算法的CPU耗時(shí)較少．

關(guān)鍵詞:有限體積法;差分格式;Newmark格式;聲波動(dòng)方程;平面波

近年來,有限體積法(FVM)被逐漸應(yīng)用于聲學(xué)問題,其目的是采用統(tǒng)一方法求解多物理場(chǎng)(如流固聲耦合)問題,避免混合方法帶來的一系列問題(如收斂困難)[1]．Fogarty等[2]嘗試將FVM應(yīng)用于一維、二維周期性材料及隨機(jī)材料等屬性變化劇烈介質(zhì)中的聲傳播問題．Del-Pino等[3]利用高階Lax-Wendroff格式提出高階FVM,并應(yīng)用于模擬地球大氣中的大尺度聲傳播問題．寧立方等[4]采用FVM求解流體Navier-Stokes方程,模擬一維諧振管內(nèi)的非線性駐波．宣領(lǐng)寬等[5]基于非結(jié)構(gòu)四面體網(wǎng)格將FVM應(yīng)用于求解三維非均勻介質(zhì)聲波動(dòng)方程,并基于時(shí)域脈沖法將其擴(kuò)展應(yīng)用于預(yù)測(cè)消聲器的傳遞損失．目前,基于FVM求解聲波動(dòng)方程的算法主要采用顯格式算法[2-5]．顯格式算法形成的矩陣通常非常稀疏,甚至不需要采用矩陣形式即可求解,其對(duì)存儲(chǔ)的消耗很低;同時(shí),不需要特定的線性方程組求解器,算法的程序?qū)崿F(xiàn)比較容易,且計(jì)算效率高. 劉恒等[6]指出提高計(jì)算效率的重要途徑是發(fā)展顯格式算法;但是,由于穩(wěn)定性條件的限制,顯格式算法需要的時(shí)間步長(zhǎng)會(huì)隨著網(wǎng)格的加密而減小,此時(shí)會(huì)造成顯格式算法的計(jì)算效率大大降低．隱格式算法的無條件穩(wěn)定卻有可能解決這一問題．

針對(duì)聲波動(dòng)方程,本文采用格點(diǎn)型FVM離散其空間項(xiàng),Newmark格式離散其時(shí)間項(xiàng),算法上采用隱格式算法．模擬平面波的傳播,對(duì)比分析本文隱格式Newmark算法與文獻(xiàn)中顯格式中心差分算法的精度、計(jì)算消耗、穩(wěn)定性及邊界條件處理正確性．

1基本方程

靜態(tài)非均勻流體中的聲波動(dòng)方程為

(1)

文中涉及兩大類邊界條件:(a)全反射邊界條件[7],包含“絕對(duì)硬”理想邊界Γv與“絕對(duì)軟”理想邊界Γp,可分別表示為p·n=0與p=0;(b)一階吸收邊界條件Γb做為對(duì)無窮遠(yuǎn)的近似,表示為[8]

2數(shù)值離散

本文采用隱格式法求解聲波動(dòng)方程,并與文獻(xiàn)[9]的顯格式法對(duì)比．二者共同點(diǎn)是聲波動(dòng)方程中的空間項(xiàng)都是基于格點(diǎn)型FVM離散,區(qū)別在于顯格式法中時(shí)間項(xiàng)采用中心差分格式離散,顯式推進(jìn)求解離散方程,而隱格式法中時(shí)間項(xiàng)采用Newmark格式離散,隱式迭代求解離散方程．式(1)的積分形式為

(2)

格點(diǎn)型FVM在離散過程中會(huì)用到單元的型函數(shù),常見的網(wǎng)格單元形式有三角形、四邊形、四面體、三棱柱、六面體[10]．本文以二維四邊形單元為例．

2.1空間離散

圖1為二維聲場(chǎng)控制體積示意圖,計(jì)算域用4節(jié)點(diǎn)四邊形單元?jiǎng)澐郑趦?nèi)部節(jié)點(diǎn)a周圍建立邊界Γint為3-4-5-6-7-8-9-10-3的控制體Ωa,在邊界節(jié)點(diǎn)b周圍建立邊界Γint為1-5-4-3-2-b-1的控制體Ωb．密度ρ、聲速c定義在單元中心,并假設(shè)其在單元內(nèi)均勻分布;待解變量聲壓p定義在節(jié)點(diǎn)上,并假設(shè)在控制體內(nèi)均勻分布．參考文獻(xiàn)[11],假設(shè)聲壓加速度在內(nèi)部控制體Ωa內(nèi)均勻分布,式(2)的左端表示為

(3)

式中:na為節(jié)點(diǎn)a周圍單元總數(shù),ρi、ci為節(jié)點(diǎn)a周圍第i個(gè)單元內(nèi)的密度、聲速,(Sa)i為節(jié)點(diǎn)a周圍第i個(gè)單元中參與控制體內(nèi)的面積,ga為節(jié)點(diǎn)a集中屬性．

式(2)的右端項(xiàng)采用高斯定理

∫Ω·(ρ-1p)dΩ=∫Γintρ-1p·ndΓ=

(4)

圖1　二維聲場(chǎng)控制體積示意

引入型函數(shù),單元內(nèi)的壓力梯度為

式中:N為型函數(shù), ?Nj/?x、?Nj/?y的表達(dá)式可參考文獻(xiàn)[9-12]．

由式(3)、(4),式(2)可變?yōu)?/p>

(5)

式中:下標(biāo)ij表示節(jié)點(diǎn)a周圍第i個(gè)單元內(nèi)第j個(gè)節(jié)點(diǎn),四邊形單元中kij的表達(dá)式見2.2節(jié)．

(6)

假設(shè)為“絕對(duì)硬”全反射邊界條件,則式(6)右端第二項(xiàng)為零,自然滿足．

假設(shè)為吸收邊界條件Γb,則式(6)右端第二項(xiàng)可表示為

(7)

式中:L2-b、Lb-1分別為邊2-b、b-1的長(zhǎng)度,(ρc)3、(ρc)5分別為單元abcd、abef內(nèi)的密度、聲速的乘積．

由式(6)、(7),式(5)可變?yōu)?/p>

式中,nbn為節(jié)點(diǎn)b周圍的節(jié)點(diǎn)總數(shù)(包含節(jié)點(diǎn)b), 則

(8)

2.2四邊形單元處理

式中,ρi為節(jié)點(diǎn)1周圍第i個(gè)四邊形單元C1234的密度,且有

式中(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)分別為四邊形單元節(jié)點(diǎn)1、2、3、4的坐標(biāo)．

(a)整體坐標(biāo)系

(b)局部坐標(biāo)系

2.3時(shí)間離散

二階時(shí)間項(xiàng)的離散可采用的隱格式算法有線性加速度法、Newmark法及Houbolt法等. 文獻(xiàn)[13]指出,線性加速度法可選擇的時(shí)間步長(zhǎng)較小,而Newmark法的計(jì)算精度比Houbolt法高．同時(shí),Newmark法是結(jié)構(gòu)動(dòng)力學(xué)中最流行的一種方法,文獻(xiàn)[11]指出,當(dāng)δ≥0.5,α≥0.25(0.5+δ)2時(shí),Newmark法是無條件穩(wěn)定的．文獻(xiàn)[14]指出當(dāng)ω0Δt<0.3時(shí)(ω0為固有頻率對(duì)應(yīng)的角頻率),數(shù)值阻尼<2%,系統(tǒng)固有頻率誤差<1%．因此本文嘗試采用Newmark法求解二階時(shí)間項(xiàng),并與格點(diǎn)型FVM相結(jié)合,將Newmark法用于求解式(8)(K-1/(αΔt2)G-δ/(αΔt)C)Pt+Δt=

(9)

依據(jù)式(9)可獲得節(jié)點(diǎn)聲壓p,節(jié)點(diǎn)聲壓速度、加速度可分別由式(10)、(11)更新:

(10)

(11)

由式(10)、(11),式(8)可變?yōu)橐怨?jié)點(diǎn)聲壓P為未知量的線性系統(tǒng),可采用迭代求解器或直接求解器進(jìn)行求解,本文采用迭代求解器代數(shù)多重網(wǎng)格法[15]求解．“絕對(duì)軟”全反射邊界條件,可采用置大數(shù)法或置“1”法等[16]進(jìn)行處理．

3兩種格式比較

為方便對(duì)比兩種格式,研究圖3所示的聲學(xué)計(jì)算域,上下為絕對(duì)硬邊界Γv,左側(cè)施加高斯壓力脈沖

其中:T′=2.5×10-3s,聲速c=500m/s,密度ρ=1kg/m3．

圖3　矩形聲場(chǎng)計(jì)算域

3.1空間步長(zhǎng)的影響

表1中4種均勻分布的四邊形單元用于網(wǎng)格無關(guān)性分析,時(shí)間步長(zhǎng)均為Δt=5×10-5s．監(jiān)測(cè)點(diǎn)M聲壓的時(shí)域響應(yīng)如圖4所示,Grid2、Grid3及Grid4的結(jié)果相互吻合良好,所以當(dāng)λ/Δx≥10,顯格式法與隱格式法均能夠獲得滿足精度要求的解,Grid3能夠用于進(jìn)一步分析．表1中給出兩種格式的內(nèi)存和CPU消耗,二者的內(nèi)存消耗區(qū)別不大,但時(shí)間步長(zhǎng)Δt相同時(shí),顯格式消費(fèi)的CPU時(shí)間更少．

表1　計(jì)算詳細(xì)信息

(a)顯格式

(b)隱格式

3.2時(shí)間步長(zhǎng)的影響

利用Grid 3分析時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響．表2中列出兩種格式點(diǎn)M聲壓的峰值誤差．針對(duì)顯格式算法,(c0Δt)/dx=0.98、1.00時(shí),計(jì)算結(jié)果穩(wěn)定且吻合良好;當(dāng)(c0Δt)/dx=1.02>1.00時(shí),計(jì)算結(jié)果發(fā)散,說明顯格式算法是條件穩(wěn)定的,且穩(wěn)定性條件與文獻(xiàn)[10]的推導(dǎo)結(jié)果基本一致．針對(duì)隱格式算法,ω0Δt=0.1、0.3時(shí),計(jì)算結(jié)果吻合良好,但ω0Δt=0.50、0.51及0.60時(shí)出現(xiàn)明顯衰減,其中ω0=2πf ′max,5種情況下的結(jié)果均沒有出現(xiàn)發(fā)散,表明隱格式算法是無條件穩(wěn)定的．當(dāng)ω0Δt=0.3時(shí),相對(duì)誤差<1%,對(duì)應(yīng)一個(gè)周期內(nèi)21個(gè)時(shí)間采樣點(diǎn),小于文獻(xiàn)[17]中33個(gè)時(shí)間采樣點(diǎn)．當(dāng)ω0Δt> 0.5時(shí),預(yù)測(cè)聲壓級(jí)差超過工程上±0.5dB的誤差需求,一個(gè)周期內(nèi)有13個(gè)時(shí)間采樣點(diǎn)．相對(duì)誤差、絕對(duì)誤差、聲壓級(jí)差分別由式(12)～(14)計(jì)算:

相對(duì)誤差=(p預(yù)測(cè)-p解析)/p解析×100%,

(12)

絕對(duì)誤差=|p預(yù)測(cè)-p解析|,

(13)

聲壓級(jí)差=20×lg10(p預(yù)測(cè)/p解析).

(14)

基于Grid3,采用3種時(shí)間步長(zhǎng),對(duì)比隱格式和顯格式方法得到點(diǎn)M聲壓的絕對(duì)誤差,見圖5．結(jié)果表明,基于相同的計(jì)算網(wǎng)格和時(shí)間步長(zhǎng),隱格式方法得到的結(jié)果與精確解吻合更好．

圖5　不同時(shí)間步長(zhǎng)時(shí)兩種格式絕對(duì)誤差的比較

3.3邊界條件的影響

分析兩種方法處理絕對(duì)軟邊界、絕對(duì)硬邊界、吸收邊界的正確性,基于Grid3,時(shí)間步長(zhǎng)Δt=1.25×10-4s,右側(cè)邊界分別為不同邊界條件時(shí)點(diǎn)M聲壓(見圖6)．當(dāng)右側(cè)邊界為絕對(duì)軟邊界時(shí),兩種格式均存在反相同幅值反射;當(dāng)右側(cè)邊界為絕對(duì)硬邊界時(shí),兩種格式均存在同相同幅值反射;當(dāng)右側(cè)邊界為吸收邊界時(shí),均將波成功透射出去．

圖7為不同邊界條件下點(diǎn)M聲壓的絕對(duì)誤差,隱格式算法得到的入射波精度高于顯格式算法;對(duì)全反射邊界條件的處理,顯格式算法的精度高于隱格式算法;對(duì)吸收邊界條件的處理,顯格式算法的精度低于隱格式算法．

(a)顯格式

(b)隱格式

圖7　不同邊界條件時(shí)兩種格式絕對(duì)誤差的比較

4結(jié)論

1) 討論了空間步長(zhǎng)的影響. 當(dāng)λ/Δx≥10后,顯格式算法與隱格式算法均能夠達(dá)到滿足精度要求的解,并且隨著網(wǎng)格的加密,計(jì)算結(jié)果保持不變;兩種格式在內(nèi)存消耗上比較接近,但在采用相同時(shí)間步長(zhǎng)的前提下,顯格式耗費(fèi)較少的CPU時(shí)間．

2) 討論時(shí)間步長(zhǎng)的影響. 顯格式算法是條件穩(wěn)定的,當(dāng)(c0Δt)/dx≤1.0時(shí),計(jì)算結(jié)果穩(wěn)定,而不滿足此條件時(shí)計(jì)算會(huì)快速發(fā)散,這與理論推導(dǎo)的結(jié)果一致. 隱格式Newmark算法是無條件穩(wěn)定的,當(dāng)ω0Δt≤0.3時(shí),所預(yù)測(cè)聲壓的相對(duì)峰值誤差<1%,當(dāng)ω0Δt>0.5時(shí),所預(yù)測(cè)聲壓級(jí)的相對(duì)峰值誤差會(huì)超出工程上±0.5 dB的誤差需求．

3) 對(duì)于文中提到的全反射邊界和吸收邊界條件,兩種算法均能夠提供峰值誤差<2%的數(shù)值解;顯格式算法所得的入射波精度高于隱格式算法;對(duì)全反射邊界條件的處理,顯格式算法處理的精度高于隱格式算法;對(duì)吸收邊界條件的處理,顯格式算法處理的精度低于隱格式算法．

參考文獻(xiàn)

[1] MANEERATANA K. Development of the finite volume method for non-linear structural applications[D]. London: the University of London, 2000.

[2] FOGARTY T R, LEVEQUE R J. High-resolution finite-volume methods for acoustic waves in periodic and random media[J]. Journal of the Acoustical Society of America, 1999, 106(1): 17-28.

[3] DEL-PINO S, DESPRéS B, HAVé P, et al. 3D Finite Volume simulation of acoustic waves in the earth atmosphere[J]. Computers & Fluids, 2009, 38(4): 765-777.[4] 寧方立, 張文治, 董梁. 圓柱形諧振管內(nèi)非線性駐波的有限體積計(jì)算方法[J]. 機(jī)械工程學(xué)報(bào), 2012, 48(16): 116-121.

[5] 宣領(lǐng)寬, 張文平, 明平劍, 等. 預(yù)測(cè)消聲器聲學(xué)性能的時(shí)域非結(jié)構(gòu)有限體積法[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2012, 33(2): 185-191.

[6] 劉恒, 廖振鵬. 波動(dòng)數(shù)值模擬的一種顯式方法——二維波動(dòng)[J]. 力學(xué)學(xué)報(bào), 2010, 42(6): 1104-1116.

[7] 杜功煥, 朱哲民, 龔秀芬. 聲學(xué)基礎(chǔ)[M]. 南京: 南京大學(xué)出版社, 2001.

[8] 李太寶. 聲場(chǎng)的方程和計(jì)算方法[M]. 北京: 科學(xué)出版社, 2005.

[9] 宣領(lǐng)寬, 張文平, 明平劍, 等. 基于不同時(shí)間步長(zhǎng)時(shí)域非結(jié)構(gòu)有限體積法模擬聲-彈性耦合問題[J]. 固體力學(xué)學(xué)報(bào), 2013, 34(2): 158-168.

[10]TAYLOR G A, BAILEY C, CROSS M. A vertex-based finite volume method applied to non-linear material problems in computational solid mechanics[J]. International Journal for Numerical Methods in Engineering, 2003, 56(4): 507-529.[11]王勖成. 有限單元法[M]. 北京: 清華大學(xué)出版社, 2003.[12]宣領(lǐng)寬, 明平劍, 龔京風(fēng), 等. 三維各向異性功能梯度材料的有限體積法[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2014, 46(7): 95-100.

[13]趙倩. 新型差分格式TDFEM及其腔體屏蔽效能分析研究[D]. 南京: 南京航空航天大學(xué), 2011.

[14]GILES M B. Stability and accuracy of numerical boundary conditions in aeroelastic analysis[J]. International Journal for Numerical Methods in Fluids, 1997, 24 (8): 739-757.[15]明平劍. 基于非結(jié)構(gòu)化網(wǎng)格氣液兩相流數(shù)值方法及并行計(jì)算研究與軟件開發(fā)[D]. 哈爾濱: 哈爾濱工程大學(xué), 2008.[16]李亞智, 趙美英, 萬小朋. 有限元法基礎(chǔ)與程序設(shè)計(jì)[M]. 北京: 科學(xué)出版社, 2004.

[17]倪大明. 非結(jié)構(gòu)化網(wǎng)格FVM在流體動(dòng)力聲學(xué)計(jì)算中的應(yīng)用研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2013.

(編輯楊波)

An implicit finite volume method for the acoustic wave equation

XUAN Lingkuan1, GONG Jingfeng1, ZHANG Wenping2, YUAN Xingjun1

(1. China Ship Development and Design Center, Wuhan 430064, China;2. College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China)

Abstract:An implicit finite volume method (FVM) with Newmark scheme is proposed for the sound propagation problem, and it is used to simulate the propagation of the plane wave. By comparing to the results from the FVM with that of explicit central difference algorithm, the numerical error, stability and computational consumption are analyzed. It is shown that both methods can acquire accurate numerical solutions when λ/Δx≥10. The relative peak error of the implicit scheme is less than 1% when ω0Δt≤0.3. With the same time step and spatial step, the numerical result of the implicit scheme agrees better with the exact result. The implicit method treats the absorbing boundary condition more accurate than the explicit method, but it is reverse for the total reflecting boundary condition. The two methods consume similar memory and the explicit method consumes less CPU time.

Keywords:finite volume method; difference scheme; Newmark scheme; acoustic wave equation; plane wave

doi：10.11918/j.issn.0367-6234.2016.07.024

收稿日期:2015-06-13

基金項(xiàng)目:國(guó)家自然科學(xué)基金(51509232)

作者簡(jiǎn)介:宣領(lǐng)寬(1987—),男,博士,工程師

通信作者:龔京風(fēng), gongjingfeng@126.com

中圖分類號(hào):O422

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):0367-6234(2016)07-0145-05