姚凌云， 于德介， 臧獻國

（湖南大學 汽車車身先進設計制造國家重點實驗室，長沙 410082）

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

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

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

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

1 二維聲學Helmholtz方程

結構振動在理想流體介質中引起的小振幅簡諧聲波滿足Helmholtz波動方程：

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

聲壓與簡諧聲波的振動速度關系：

內部聲場的邊界條件為：

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

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

2.1 徑向點插值法

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

式中：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ū)梯度光滑處理

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

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

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

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

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

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

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

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

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

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

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

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

M表示聲學質量矩陣，可以寫成集中質量的形式：

F表示聲壓載荷矢量：

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

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

3 數(shù)值算例

3.1 管道聲場分析

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

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

管道聲壓的精確解為：

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

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

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

圖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結果精度和收斂性的比較Fig.5 Comparison of accuracy and convergence property at different frequency values between CS-RPIM and FEM

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

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

3.2 轎車車內聲場分析

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

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

4 結論

本文推導了CS-RPIM求解聲學波動方程的計算公式，將該方法應用到管道和車內聲場的仿真分析中，研究結果表明：

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

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

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

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

［1］朱才朝，秦大同，李潤方.車身結構振動與車內噪聲聲場耦合分析與控制［J］.機械工程學報，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］程 昊，高 煌，張永斌，等.振動體聲學靈敏度分析的邊界元法［J］.機械工程學報，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.