聲學(xué)數(shù)值計算的分區(qū)光滑徑向點插值無網(wǎng)格法

姚凌云， 于德介， 臧獻(xiàn)國

（湖南大學(xué) 汽車車身先進(jìn)設(shè)計制造國家重點實驗室，長沙 410082）

近年來，噪聲預(yù)測和控制技術(shù)在車輛、飛機(jī)、潛艇等工程領(lǐng)域中得到越來越廣泛的應(yīng)用。聲學(xué)數(shù)值計算是實現(xiàn)噪聲控制和預(yù)測的關(guān)鍵技術(shù)，是進(jìn)行結(jié)構(gòu)噪聲預(yù)測和優(yōu)化設(shè)計的基礎(chǔ)。針對Helmholtz聲學(xué)方程的數(shù)值計算方法研究是近年來聲學(xué)研究的一個熱點，其研究重點在于如何提高聲學(xué)數(shù)值計算的精度、計算效率、正確性以及算法的適用性。聲學(xué)數(shù)值方法主要有有限元法（Finite element method，F(xiàn)EM）［1-5］、邊界元法［6，7］和無網(wǎng)格（Meshfree）技術(shù)［8-10］等。

運(yùn)用FEM求解聲波方程時通常會遇到“數(shù)值色散”問題，即存在數(shù)值波的相位誤差，數(shù)值仿真結(jié)果在高波數(shù)時色散嚴(yán)重。為了獲得精度可靠的仿真結(jié)果，通常需要更精細(xì)單元和更高階的多項式近似函數(shù)，這樣勢必增加計算時間和存儲空間。為了降低色散效應(yīng)，Petersen等［11］將高階譜單元形函數(shù)應(yīng)用到內(nèi)部聲學(xué)分析中，提高其計算精度和效率；Harari等［12］在聲學(xué)波動方程的基本Galerkin形式下采用穩(wěn)定化有限元法（Stabilized finite element method）提高數(shù)值方法的穩(wěn)定性。

近年來，無網(wǎng)格技術(shù)在聲學(xué)波動方程計算中有很大的應(yīng)用和發(fā)展。Bouillard等［9］將無單元迦遼金法應(yīng)用到聲學(xué)波動方程的計算中，其結(jié)果的色散誤差明顯小于FEM結(jié)果。由于無網(wǎng)格法存在節(jié)點積分的不穩(wěn)定性，Chen 等［13］提出了光滑應(yīng)變技術(shù)。Liu［14］等將應(yīng)變光滑技術(shù)和徑向點插值法相結(jié)合，提出了節(jié)點光滑徑向點插值法（Node-based conforming radial point interpolation method，NS-RPIM）；將有限元法與無網(wǎng)格中分區(qū)應(yīng)變光滑技術(shù)相結(jié)合提出光滑有限元法（Smoothed finite element method，SFEM）［15］。Cui［16］用 SFEM 劃分光滑域的思想與徑向點插值法相結(jié)合，提出分區(qū)光滑徑向點插值法（Cell-based smoothed radial point interpolation method，CS-RPIM），該方法提供合適的模型剛度，能很好地解決力學(xué)計算問題。

在聲學(xué)數(shù)值計算中，由于“數(shù)值色散”效應(yīng)導(dǎo)致計算誤差隨著波數(shù)k的增加變大。針對此問題，本文將CS-RPIM推廣到聲學(xué)問題的數(shù)值計算中，推導(dǎo)了CSRPIM計算二維聲學(xué)方程的公式，利用分區(qū)光滑梯度技術(shù)對聲壓梯度進(jìn)行光滑處理，只需對光滑域邊界進(jìn)行積分，有效降低因數(shù)值色散而引起的計算誤差。為驗證CS-RPIM求解聲學(xué)波動方程的有效性，本文分析了二維管道聲腔和車內(nèi)聲腔模型兩個算例，結(jié)果顯示：CS-RPIM分析聲學(xué)問題時比FEM具有更高的精度和更好的收斂性，數(shù)值色散誤差更小，在較高波數(shù)下仍具有較高的精度。

1 二維聲學(xué)Helmholtz方程

結(jié)構(gòu)振動在理想流體介質(zhì)中引起的小振幅簡諧聲波滿足Helmholtz波動方程：

式中：▽為拉普拉斯算子；k=ω/c為波數(shù)，ω表示圓頻率；c表示聲速。

聲壓與簡諧聲波的振動速度關(guān)系：

內(nèi)部聲場的邊界條件為：

式中：p為邊界處聲壓；n為腔邊界表面法線方向；An為聲導(dǎo)納系數(shù)。

2 聲學(xué)分區(qū)光滑徑向點插值法基本理論

2.1 徑向點插值法

對聲場域內(nèi)任意場點x聲壓p（x）采用耦合多項式基的徑向點插值進(jìn)行插值，有：

式中：Ri（x）是徑向基函數(shù)，n是徑向基函數(shù)的項數(shù)，qj（x）是多項式基函數(shù)，m是多項式基函數(shù)的項數(shù)，ai、bj為待定系數(shù).。令近似聲壓向量p（x）為節(jié)點聲壓P，可得：

為求解系數(shù)ai和bj， 令：

聯(lián)合方程式（7）和式（8）得：

求解方程式（9）并將求解的未知量代入式（6）：

式中：φi（x）為節(jié)點形函數(shù)，pi為節(jié)點向量。

2.2 分區(qū)梯度光滑處理

分析二維聲學(xué)問題時先將問題域Ω離散為容易生成的三角形背景單元，然后將三角形單元劃分為若干個光滑區(qū)域（Smoothing cells，SC），滿足 Ω1∪Ω2∪…ΩSC=Ω和Ω1∩Ω2∩…ΩSC=Ψ，其中Ψ表示空集；文獻(xiàn)［16］中對背景單元光滑域的個數(shù)和形狀進(jìn)行了研究，研究光滑域的個數(shù)越多，積分就越硬，而一般有限元法全積分相當(dāng)于光滑域為無窮多個，即SC=∞；當(dāng)光滑域為1時，相當(dāng)有限元法的縮減積分。在聲學(xué)數(shù)值計算中光滑域的個數(shù)一般取3為宜。三角形背景單元劃分光滑域個數(shù)1，3和4方式如圖1所示。

圖1 背景單元劃分光滑域示意圖Fig.1 Division of background cell into SC smoothing cells

由式（7）可知聲壓梯度光滑處理的實質(zhì)是速度光滑處理。在第C個光滑區(qū)域內(nèi)，光滑速度表示為：

本文采用常值的光滑函數(shù)，如式（12）所示：

將式（12）代入式（11）并利用分步積分公式，在光滑域內(nèi)的速度積分變?yōu)橛蜻吔缟下晧悍e分：

式中：為光滑域邊界；n為光滑域邊界法向單位矢量。將式（10）代入式（13），光滑速度寫成：

式中：N1為光滑域邊界段個數(shù)；xij為第i邊界段的第j個高斯積分點；nid為第i邊界段的外法向矢量；li為第i邊界段長度。

2.3 系統(tǒng)方程離散

CS-RPIM采用標(biāo)準(zhǔn)的Galerkin方法與聲壓梯度光滑處理相結(jié)合。在整個問題域里對Helmholtz方程及其邊界條件進(jìn)行離散處理，得到光滑的Galerkin弱形式：

將式（10）和式（14）代入式（16）中，得到聲域離散系統(tǒng)方程：

式中：Ncell為背景單元個數(shù)；Kij（k）為與單元Ωk相關(guān)聯(lián)的剛度矩陣，表示為：

C表示由Robin邊界條件得到的聲學(xué)阻尼矩陣：

M表示聲學(xué)質(zhì)量矩陣，可以寫成集中質(zhì)量的形式：

F表示聲壓載荷矢量：

P表示節(jié)點聲壓矢量：

由上面的推導(dǎo)可以看出，應(yīng)用梯度光滑處理可以使光滑聲學(xué)梯度矩陣的計算簡單，使區(qū)域的面積分變成邊界線積分。

3 數(shù)值算例

3.1 管道聲場分析

本節(jié)應(yīng)用CS-RPIM分析文獻(xiàn)［17，18］中算例模型，對一管道聲場進(jìn)行計算分析。調(diào)整管道聲場參數(shù)：管長1 m、寬0.1 m；聲場內(nèi)部由空氣流體填充，其密度ρ為1.225 kg/m3。聲波在空氣中的波速c=343 m/s，管子一端施加速度邊界條件：法向速度vn=0.1 sin（ωt）；另外一端為剛性壁，如圖2所示。

圖2 管道聲學(xué)模型Fig.2 The model of tube acoustic problem

管道聲壓的精確解為：

運(yùn)用CS-RPIM分析背景單元尺寸為25 mm的不同頻率（波數(shù)）下的聲壓分布。為了便于比較，應(yīng)用FEM分析相同的網(wǎng)格模型。圖3和圖4分別表示管道中心線x方向的聲壓（虛部）分布。

為分析CS-RPIM結(jié)果的精度和收斂性，評價聲學(xué)計算的相對誤差：

為分析聲學(xué)網(wǎng)格的規(guī)則程度對CS-RPIM和FEM計算精度的影響，用不規(guī)則度參數(shù)來衡量網(wǎng)格的不規(guī)則程度。對于二維管道聲域問題，對節(jié)點坐標(biāo)進(jìn)行如下變換：

圖3 波數(shù)為10的聲壓虛部分布Fig.3 Spatial distribution of real part of the pressure at k=10

圖4 波數(shù)為30的聲壓虛部分布Fig.4 Spatial distribution of real part of the pressure at k=30

圖5 CS-RPIM和FEM結(jié)果精度和收斂性的比較Fig.5 Comparison of accuracy and convergence property at different frequency values between CS-RPIM and FEM

式中：x'和y'為節(jié)點不規(guī)則變動后的坐標(biāo)；x和y為原始坐標(biāo)；Δx和Δy為初始規(guī)則單元節(jié)點之間的x和y方向的距離；rC是由計算機(jī)在（0.0，1.0）之間產(chǎn)生的隨機(jī)數(shù)；βir為表示網(wǎng)格不規(guī)則度的參數(shù)，它在0.0和0.5之間。βir越大表明網(wǎng)格不規(guī)則程度越高。圖6（a）、圖6（b）分別表示βir為0和0.5時二維管道聲場有限元網(wǎng)格。

應(yīng)用CS-RPIM對兩種網(wǎng)格模型計算聲壓響應(yīng)。圖7為波數(shù)為20時管道聲場底部邊界線的聲壓分布。圖中CS-RPIM（re）表示規(guī)則背景網(wǎng)格模型下應(yīng)用CSRPIM的計算結(jié)果；CS-RPIM（irr）表示不規(guī)則背景網(wǎng)格模型下應(yīng)用CS-RPIM的計算結(jié)果；FEM-T3（re）表示規(guī)則三角形網(wǎng)格模型下應(yīng)用FEM計算結(jié)果；FEMT3（irr）表示不規(guī)則三角形網(wǎng)格模型的FEM結(jié)果。從圖中可以發(fā)現(xiàn)，CS-RPIM對聲學(xué)背景單元的網(wǎng)格不規(guī)則敏感度很低，在背景網(wǎng)格扭曲嚴(yán)重時仍然具有很高的精度。

3.2 轎車車內(nèi)聲場分析

本節(jié)運(yùn)用CS-RPIM分析轎車聲腔模型。將車內(nèi)聲腔簡化為二維模型，座椅簡化為內(nèi)部的一個空腔；邊界條件為在聲腔靠近發(fā)動機(jī)邊界上法向振動速度vn=0.1 m/s，如圖8所示。由HyperMesh自動生成三角形網(wǎng)格，單元數(shù)為347，節(jié)點數(shù)為219，如圖9所示。空氣密度為1.225 kg/m3；聲速為343 m/s。

應(yīng)用CS-RPIM和FEM分別對車內(nèi)聲腔網(wǎng)格進(jìn)行計算，并對不同頻率（波數(shù)）下的計算結(jié)果（聲壓虛部的分布）進(jìn)行對比。圖10和圖11分別表示波速k=5和10時路徑ab上的聲壓（虛部）分布。用FEM計算節(jié)點數(shù)為17 430的模型結(jié)果作為參考值。從圖中可以發(fā)現(xiàn)，與FEM相比，運(yùn)用CS-RPIM的計算結(jié)果精度隨波數(shù)的增加變化不大，即CS-RPIM計算精度結(jié)果對波數(shù)k的靈敏度比FEM要小。

4 結(jié)論

本文推導(dǎo)了CS-RPIM求解聲學(xué)波動方程的計算公式，將該方法應(yīng)用到管道和車內(nèi)聲場的仿真分析中，研究結(jié)果表明：

（1）CS-RPIM在計算聲學(xué)問題時能提供合適的模型硬度，有效降低因數(shù)值色散而引起的計算誤差，因此其收斂性和結(jié)果精度比FEM更高。

（2）與FEM相比，CS-RPIM對計算波數(shù)（即頻率）的敏感性較低，在較高波數(shù)時計算結(jié)果精度仍然很好。

（3）CS-RPIM對背景單元的質(zhì)量要求很低，可以對扭曲嚴(yán)重的背景網(wǎng)格進(jìn)行計算分析，因此可以在前處理分析中節(jié)省大量的人力和時間，具有良好的工程應(yīng)用前景。

（4）由于CS-RPIM只是改進(jìn)了聲學(xué)剛度矩陣的計算，降低數(shù)值色散效應(yīng)，因此CS-RPIM的計算結(jié)果受到計算頻率、自由度和階次的影響。

［1］朱才朝，秦大同，李潤方.車身結(jié)構(gòu)振動與車內(nèi)噪聲聲場耦合分析與控制［J］.機(jī)械工程學(xué)報，2002，38（8）：54-58.

［2］ Bales L，Lasiecka I.Continuous finite elements in space and time for the nonhomogeneous wave equation［J］.Computers＆ Mathematics with Applications，1994，27 91-102.

［3］ Ihlenburg F，Babu?ka I.Finite element solution of the helmholtz equation with high wave number.Part 1：the hversion of the FEM［J］.Computers＆ Mathematics with Applications，1995，30：9-37.

［4］Richter G R.An explicit finite element method for the wave equation ［J］. AppliedNumericalMathematics，1994，16（1-2）：65-80.

［5］ Jenkins E W，Riviěre B，Wheeler M F.A priori error estimates for mixed finite element approximations of the acoustic wave equation［J］.SIAM Journal on Numerical Analysis，2002，40（5）：1698-1715.

［6］Chen J T，Chen K H，Chyuan S W.Numerical experiments for acoustic modes of a square cavity using the dual boundary element method［J］.Applied Acoustic，1999，57：293-325.

［7］程 昊，高 煌，張永斌，等.振動體聲學(xué)靈敏度分析的邊界元法［J］.機(jī)械工程學(xué)報，2008，44（7）：45-51.

［8］Wenterodt C，Estorff O V.Dispersion analysis of the meshfree radial point interpolation method for the Helmholtz equation［J］. InternationalJournalforNumericalMethods in Engineering，2009，77：1670-1689.

［9］Bouillard P，Suleau S.Element-free Garlekin solutions for Helmholtz problems：formulation and numerical assessment of the pollution effect［J］.Computer Methods in Applied Mechanicals and Engineering，1998，162：317-335.

［10］ Suleau S，Deraemaeker A，Bouillard P.Dispersion and pollution of meshless solution for the Helmholtz equation［J］.Computer Methods in Applied Mechanicals and Engineering，2000，190：639-657.

［11］ Petersen S，Dreyer D，Estorff O V.Assessment of finite and spectralelementshape functions or efficientiterative simulations of interior acoustics［J］.Computer Methods in Applied Mechanicals and Engineering， 2006， 195：6463-6478.

［12］ Harari I，Magoules F.Numerical investigations of stabilized finite element computations for acoustics［J］.Wave Motion，2004，39：339-349.

［13］ Chen J S，Wu CT，Yoon S，et al.A stabilized conforming nodal integration forGalerkin meshfree methods［J］.International Journal for Numerical Methods in Engineering，2001，50：435-466.

［14］ Liu G R，Li Y，Dai K Y，et al.A linearly conforming radial point interpolation method for solid mechanics problems［J］.International Journal of Computational Methods，2006，3（4）：401-428.

［15］ Liu G R，Nguyen T T，Dai K Y，et al.Theoretical aspects of thesmoothed finite elementmethod （SFEM）［J］.International Journal for Numerical Methods in Engineering，2007，71：902-930.

［16］ Cui X Y，Liu G R，Li G Y，et al.A cell-based smoothed radial point interpolation method（CS-RPIM）for static and free vibration of solids［J］.Engineering analysis with boundary elements，2010，34（2）：144-1570.

［17］ Bouillard P. Admissible fields and error estimation for acousticFEA with low wavenumbers.Computersand Structures，1999，73：227-237.

［18］ He Z C，Liu G R，Zhong Z H，et al.An edge-based smoothed finite elementmethod （ES-FEM）foranalyzingthreedimensional acoustic problems［J］.Computer Methods in Applied Mechanicals and Engineering，2009，199：20-33.