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

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

（湖南大學(xué) 汽車(chē)車(chē)身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410082）

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

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

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

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

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

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

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

聲壓與簡(jiǎn)諧聲波的振動(dòng)速度關(guān)系：

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

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

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

2.1 徑向點(diǎn)插值法

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

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

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

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

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

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

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

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

圖1 背景單元?jiǎng)澐止饣蚴疽鈭DFig.1 Division of background cell into SC smoothing cells

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

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

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

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

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

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

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

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

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

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

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

F表示聲壓載荷矢量：

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

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

3 數(shù)值算例

3.1 管道聲場(chǎng)分析

本節(jié)應(yīng)用CS-RPIM分析文獻(xiàn)［17，18］中算例模型，對(duì)一管道聲場(chǎng)進(jìn)行計(jì)算分析。調(diào)整管道聲場(chǎng)參數(shù)：管長(zhǎng)1 m、寬0.1 m；聲場(chǎng)內(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分別表示管道中心線(xiàn)x方向的聲壓（虛部）分布。

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

為分析聲學(xué)網(wǎng)格的規(guī)則程度對(duì)CS-RPIM和FEM計(jì)算精度的影響，用不規(guī)則度參數(shù)來(lái)衡量網(wǎng)格的不規(guī)則程度。對(duì)于二維管道聲域問(wèn)題，對(duì)節(jié)點(diǎn)坐標(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é)點(diǎn)不規(guī)則變動(dòng)后的坐標(biāo)；x和y為原始坐標(biāo)；Δx和Δy為初始規(guī)則單元節(jié)點(diǎn)之間的x和y方向的距離；rC是由計(jì)算機(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時(shí)二維管道聲場(chǎng)有限元網(wǎng)格。

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

3.2 轎車(chē)車(chē)內(nèi)聲場(chǎng)分析

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

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

4 結(jié)論

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

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

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

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

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

［1］朱才朝，秦大同，李潤(rùn)方.車(chē)身結(jié)構(gòu)振動(dòng)與車(chē)內(nèi)噪聲聲場(chǎng)耦合分析與控制［J］.機(jī)械工程學(xué)報(bào)，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］程 昊，高 煌，張永斌，等.振動(dòng)體聲學(xué)靈敏度分析的邊界元法［J］.機(jī)械工程學(xué)報(bào)，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.