求解聲波動方程的隱格式有限體積法

宣領寬, 龔京風, 張文平, 袁興軍

(1.中國艦船研究設計中心, 武漢 430064; 2.哈爾濱工程大學 動力與能源工程學院, 哈爾濱 150001)

求解聲波動方程的隱格式有限體積法

宣領寬1, 龔京風1, 張文平2, 袁興軍1

(1.中國艦船研究設計中心, 武漢 430064; 2.哈爾濱工程大學 動力與能源工程學院, 哈爾濱 150001)

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

關鍵詞:有限體積法;差分格式;Newmark格式;聲波動方程;平面波

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

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

1基本方程

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

(1)

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

2數值離散

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

(2)

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

2.1空間離散

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

(3)

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

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

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

(4)

圖1　二維聲場控制體積示意

引入型函數,單元內的壓力梯度為

式中:N為型函數, ?Nj/?x、?Nj/?y的表達式可參考文獻[9-12]．

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

(5)

式中:下標ij表示節(jié)點a周圍第i個單元內第j個節(jié)點,四邊形單元中kij的表達式見2.2節(jié)．

(6)

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

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

(7)

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

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

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

(8)

2.2四邊形單元處理

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

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

(a)整體坐標系

(b)局部坐標系

2.3時間離散

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

(9)

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

(10)

(11)

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

3兩種格式比較

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

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

圖3　矩形聲場計算域

3.1空間步長的影響

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

表1　計算詳細信息

(a)顯格式

(b)隱格式

3.2時間步長的影響

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

相對誤差=(p預測-p解析)/p解析×100%,

(12)

絕對誤差=|p預測-p解析|,

(13)

聲壓級差=20×lg10(p預測/p解析).

(14)

基于Grid3,采用3種時間步長,對比隱格式和顯格式方法得到點M聲壓的絕對誤差,見圖5．結果表明,基于相同的計算網格和時間步長,隱格式方法得到的結果與精確解吻合更好．

圖5　不同時間步長時兩種格式絕對誤差的比較

3.3邊界條件的影響

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

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

(a)顯格式

(b)隱格式

圖7　不同邊界條件時兩種格式絕對誤差的比較

4結論

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

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

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

參考文獻

[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] 寧方立, 張文治, 董梁. 圓柱形諧振管內非線性駐波的有限體積計算方法[J]. 機械工程學報, 2012, 48(16): 116-121.

[5] 宣領寬, 張文平, 明平劍, 等. 預測消聲器聲學性能的時域非結構有限體積法[J]. 哈爾濱工程大學學報, 2012, 33(2): 185-191.

[6] 劉恒, 廖振鵬. 波動數值模擬的一種顯式方法——二維波動[J]. 力學學報, 2010, 42(6): 1104-1116.

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

[8] 李太寶. 聲場的方程和計算方法[M]. 北京: 科學出版社, 2005.

[9] 宣領寬, 張文平, 明平劍, 等. 基于不同時間步長時域非結構有限體積法模擬聲-彈性耦合問題[J]. 固體力學學報, 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]. 北京: 清華大學出版社, 2003.[12]宣領寬, 明平劍, 龔京風, 等. 三維各向異性功能梯度材料的有限體積法[J]. 哈爾濱工業(yè)大學學報, 2014, 46(7): 95-100.

[13]趙倩. 新型差分格式TDFEM及其腔體屏蔽效能分析研究[D]. 南京: 南京航空航天大學, 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]明平劍. 基于非結構化網格氣液兩相流數值方法及并行計算研究與軟件開發(fā)[D]. 哈爾濱: 哈爾濱工程大學, 2008.[16]李亞智, 趙美英, 萬小朋. 有限元法基礎與程序設計[M]. 北京: 科學出版社, 2004.

[17]倪大明. 非結構化網格FVM在流體動力聲學計算中的應用研究[D]. 哈爾濱: 哈爾濱工程大學, 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

基金項目:國家自然科學基金(51509232)

作者簡介:宣領寬(1987—),男,博士,工程師

通信作者:龔京風, gongjingfeng@126.com

中圖分類號:O422

文獻標志碼:A

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