陳 寧, 于德介, 呂 輝, 夏百戰(zhàn)

(湖南大學(xué) 汽車車身先進(jìn)設(shè)計制造國家重點(diǎn)實(shí)驗(yàn)室,長沙 410082)

幾乎所有的聲學(xué)問題都與結(jié)構(gòu)-聲場耦合系統(tǒng)有關(guān)，因此對結(jié)構(gòu)-聲場耦合系統(tǒng)的研究具有重要的工程意義。結(jié)構(gòu)-聲場耦合系統(tǒng)的分析能為結(jié)構(gòu)件的優(yōu)化提供重要信息，特別是容易受到聲壓激勵而產(chǎn)生振動的彈性薄壁結(jié)構(gòu)件。板結(jié)構(gòu)-聲場耦合系統(tǒng)包含板結(jié)構(gòu)域、聲場域和兩個域之間的耦合作用。通常，在結(jié)構(gòu)域用位移描述結(jié)構(gòu)的狀態(tài)，在聲場域用聲壓描述聲場的狀態(tài)[1-2]。

目前，有限元/有限元法(Finite Element Method, FEM)和有限元(FEM)/邊界元法(Boundary Element Method, BEM)是分析結(jié)構(gòu)-聲場耦合系統(tǒng)最常用的數(shù)值計算方法。FEM/ BEM與FEM/ FEM相比，能減少模型的單元數(shù)量，但其矩陣為非稀疏矩陣，并不一定能提高計算效率；而對于內(nèi)聲場問題，F(xiàn)EM/ FEM 計算效率往往更高一些；同時兩者的計算結(jié)果受模型網(wǎng)格尺寸影響較大，因此FEM/ FEM 和FEM/ BEM 主要用于求解中低頻的結(jié)構(gòu)-聲場耦合 問題[3-4]。Yao等[4]提出的SFEM(Smoothed Finite Element Method)/FEM方法在結(jié)構(gòu)-聲場耦合問題的分析中能取得比FEM/FEM更高的精度。

有限元法雖然得到了廣泛的應(yīng)用，但仍存在一些不足，如等參單元對于網(wǎng)格幾何變形非常敏感，為了獲得可靠的計算結(jié)果，通常需要更精細(xì)的網(wǎng)格等。無網(wǎng)格法則有效地避免了有限元的一些缺點(diǎn)，具有精度高、 計算模型不需劃分網(wǎng)格等特點(diǎn)，但無網(wǎng)格法本身存在諸多缺陷，如其形函數(shù)不具備克羅內(nèi)克爾性質(zhì)，導(dǎo)致不能直接施加邊界條件、計算效率降低等。針對有限元法和無網(wǎng)格技術(shù)的特點(diǎn)，Melenk等[5-6]結(jié)合有限元和無網(wǎng)格技術(shù)，提出了混合有限元-無網(wǎng)格法來分析各類力學(xué)問題。其中，單位分解有限元法是一種常用的混合有限元-無網(wǎng)格方法，該方法的基本原理是在不增加支撐點(diǎn)的前提下，通過增加局部支撐函數(shù)的階次來構(gòu)造高階的全局有限元公式。Zhang等[7-8]提出了一種基于單位分解的有限元-最小二乘點(diǎn)插值法(Finite Element-Least Square Point Interpolation Method，F(xiàn)E-LSPIM)，該方法采用的有限元-無網(wǎng)格四邊形單元將有限元形函數(shù)和最小二乘點(diǎn)插值形函數(shù)相結(jié)合，綜合了有限元法和無網(wǎng)格法各自的優(yōu)點(diǎn)，其形函數(shù)具有克羅內(nèi)克爾性質(zhì)，具有單元兼容性以及高階完備性，并且成功地應(yīng)用于靜力學(xué)和動力學(xué)分析中。隨后姚凌云等[9]將FE-LSPIM推廣到二維聲場的研究之中，并取得了良好的效果。

為了提高結(jié)構(gòu)-聲場耦合分析精度，降低結(jié)構(gòu)-聲場耦合分析中對結(jié)構(gòu)網(wǎng)格尺寸的要求，本文將FE-LSPIM推廣到板結(jié)構(gòu)動力學(xué)和三維聲場的分析中，提出用于板結(jié)構(gòu)-聲場耦合問題分析的FE-LSPIM/FE-LSPIM。運(yùn)用FE-LSPIM/FE-LSPIM分析板結(jié)構(gòu)-聲場耦合問題時，板結(jié)構(gòu)域采用有限元-無網(wǎng)格四邊形單元，聲場域采用有限元-無網(wǎng)格六面體單元。以一個六面體板結(jié)構(gòu)-聲場耦合模型為數(shù)值算例進(jìn)行分析，結(jié)果表明，F(xiàn)E-LSPIM/FE-LSPIM繼承了FE-LSPIM適用性好、精度高的特點(diǎn)。與FEM/ FEM和SFEM/ FEM相比，F(xiàn)E-LSPIM/FE-LSPIM的精度更高，對網(wǎng)格尺寸的質(zhì)量要求更低，能很好地應(yīng)用于板結(jié)構(gòu)- 聲場耦合分析，具有良好的工程應(yīng)用前景。

1 基于FE-LSPIM/FE-LSPIM的板結(jié)構(gòu)動力學(xué)分析

1.1 板單元的有限元-最小二乘點(diǎn)插值

根據(jù)剪切變形理論，如圖1所示，明德林-瑞斯納 (Mindlin-Reissner)板模型的位移分量u,v分別表示為：

u=-wθx(x,y)

v=-wθy(x,y)

(1)

式中:u,v,w分別為板中面x,y,z三個方向的位移，θx和θy分別為xz和yz平面內(nèi)的轉(zhuǎn)角。定義彎曲應(yīng)變?yōu)棣剩羟袘?yīng)變?yōu)棣?，有?/p>

(2)

板結(jié)構(gòu)的橫向剪切剛度本構(gòu)矩陣Ds和彎曲剛度本構(gòu)矩陣Db分別可表示為：

(3)

式中，E為彈性模量；t為板單元厚度；ν為柏松比；μ為剪切修正系數(shù)[1]。

圖1 板單元示意圖

根據(jù)標(biāo)準(zhǔn)的有限元法，對板結(jié)構(gòu)域Ω進(jìn)行離散，四邊形單元數(shù)為Ne，節(jié)點(diǎn)數(shù)為Nd。對于板單元，廣義位移u={θxθyw}T是各自獨(dú)立插值的，其表示為：

u(x,y)=Nplateue

(4)

式中:Nplate為四邊形板單元形函數(shù)，ue為節(jié)點(diǎn)的位移，可以寫成：

(I=1,2,3,4)

(5)

ue={u1(x,y)u2(x,y)u3(x,y)u4(x,y)}T

ui(x,y)={θxiθyiwi}T(i=1,2,3,4)

(6)

節(jié)點(diǎn)位移函數(shù)ui(x,y)由支撐節(jié)點(diǎn)通過LSPIM插值得到。以廣義位移中的撓度對整個求解過程進(jìn)行推導(dǎo)，撓度w(x,y)可表示為：

w(x,y)=NIwe

(7)

向量we表示對應(yīng)的四邊形單元四個節(jié)點(diǎn)的撓度近似函數(shù)wi(x,y)(i=1,2,3,4)，可表示為：

we={w1(x,y)w2(x,y)w3(x,y)w4(x,y)}T

(8)

其中撓度wi(x,y)(i=1,2,3,4)可表示為：

wi(x,y)=Φiwi(i=1,2,3,4)

(9)

wi=[w1w2w3…wn]T

(10)

式中，Φi為節(jié)點(diǎn)i的LSPIM形函數(shù)，它由節(jié)點(diǎn)i的支撐域點(diǎn)通過LSPIM構(gòu)成；wi為支撐節(jié)點(diǎn)的撓度向量函數(shù)；n為節(jié)點(diǎn)i的支撐域點(diǎn)的個數(shù)。

矩陣Φ4×M由Φi(i=1, 2, 3, 4)組裝得到，列數(shù)M等于單元支撐域的節(jié)點(diǎn)數(shù)。節(jié)點(diǎn)支撐域Ω1={1, 2, 3, 4, 5, 6, 7, 8, 9},Ω2={1, 2, 3, 4, 8, 9, 10, 11, 12},Ω3={1, 2, 3, 4, 11, 12, 13, 14, 15} ，Ω4={1, 2, 3, 4, 5, 6, 14, 15, 16}；單元支撐域Ω=Ω1⊕Ω2⊕Ω3⊕Ω4={1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,12, 13, 14, 15, 16}，如圖2所示。

圖2 節(jié)點(diǎn)及單元支撐域定義

將式(10),式(9),式(8)代入式(7)中，得到四邊形單元內(nèi)場點(diǎn)撓度變量的近似形式為：

(11)

式中，Ψ為FE-LSPIM的形函數(shù)矩陣，定義為：

(12)

同理，轉(zhuǎn)動θx和θy可分別表示為：

(13)

1.2 板的動力學(xué)方程

無阻尼板結(jié)構(gòu)動力學(xué)方程的Galerkin弱形式為：

∫ΩδκTDbκdΩ+∫ΩδγTDsdΩ+

(14)

離散后板結(jié)構(gòu)的動力學(xué)方程為：

(15)

式中，K=Kb+Ks為板單元剛度矩陣；Kb為彎曲剛度矩陣，Ks為剪切剛度矩陣，M為板單元質(zhì)量矩陣，F(xiàn)b為體積力列陣。

剛度矩陣K可表示為：

K=Kb+Ks=

∫Ω(Bb)TDbBbdΩ+∫Ω(Bs)TDsBsdΩ

(16)

式中：

(17)

質(zhì)量矩陣M可表示為：

(18)

載荷矢量可表示為：

Fb=∫ΩQTbdΩ

(19)

2 板結(jié)構(gòu)-聲場耦合分析

2.1 聲場有限元-最小二乘點(diǎn)插值模型

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

2p+k2p=0

(20)

式中，p為聲壓；k為波數(shù)，k=ω/c；ω為圓頻率；c為聲速。

根據(jù)伽遼金原理，聲場問題的弱形式可以寫成如下形式：

∫Ωδp·

(21)

式中，qf表示單位體積的附加載荷。

按照FEM方式將聲場域離散為六面體網(wǎng)格，節(jié)點(diǎn)聲壓值近似寫成：

(22)

式中，m為離散單元的節(jié)點(diǎn)個數(shù)，Nf為標(biāo)準(zhǔn)有限元流體單元的形函數(shù)，向量p為六面體單元8個節(jié)點(diǎn)的聲壓近似函數(shù)pi(x,y)(i=1,2,…,8)。節(jié)點(diǎn)聲壓函數(shù)pi(x,y)由支撐節(jié)點(diǎn)通過LSPIM插值得到，其支撐域范圍如圖3所示。其插值過程與板單元的有限元-最小二乘插值相同，在此不再贅述。

圖3 六面體單元節(jié)點(diǎn)i的支撐域點(diǎn)

有限元-最小二乘插值得到的六面體單元內(nèi)的聲壓近似函數(shù)為：

(23)

將式(22)代入(21)中，得到系統(tǒng)的離散方程：

p={p1,p2,…,pn}T

Fs=ρ∫ΩsfΨTΨdΓ

(24)

式中，Kf為聲學(xué)剛度矩陣，Bf為聲學(xué)梯度矩陣，Mf為聲學(xué)質(zhì)量矩陣，p為節(jié)點(diǎn)的聲壓矢量，F(xiàn)s為載荷向量。

2.2 結(jié)構(gòu)-聲場耦合模型

為實(shí)現(xiàn)結(jié)構(gòu)FE-LSPIM模型和聲場FE-LSPIM模型的耦合，如圖4所示，在板結(jié)構(gòu)與聲場的交界面即耦合界面上，應(yīng)滿足位移和壓力連續(xù)的條件。引入界面法向矢量n=nf=-ns，位移連續(xù)條件和壓力連續(xù)條件可表示為：

(25)

圖4 板結(jié)構(gòu)-聲場耦合系統(tǒng)示意圖

根據(jù)式(24)，結(jié)構(gòu)作用在聲場耦合界面Ωsf上的載荷Fs為：

(26)

根據(jù)式(19)，聲場作用在結(jié)構(gòu)耦合界面Ωsf上的載荷Ff為：

(27)

式中Ns為結(jié)構(gòu)域單元的形函數(shù)，Nf是聲場單元的形函數(shù)，這里都采用四邊形等參單元的形函數(shù)。

引入耦合矩陣L：

(28)

從結(jié)構(gòu)域的FE-LSPIM模型和聲場的FE-LSPIM模型，可得到結(jié)構(gòu)-聲場耦合的FE-LSPIM/ FE-LSPIM模型：

(29)

假設(shè)位移和聲壓為整個時域上的諧波，上式可以寫成：

(30)

3 數(shù)值算例

圖5 帶矩形板的六面體結(jié)構(gòu)-聲場耦合模型

圖5所示為六面體附帶矩形板的結(jié)構(gòu)-聲學(xué)耦合模型，尺寸為0.414 m×0.314 m×0.360 m。六面體內(nèi)為空氣域，下表面為矩形平面板結(jié)構(gòu)，尺寸0.414 m×0.314 m，板四邊簡支，其它表面為剛性邊界。矩形平面板材料參數(shù)為：彈性模量E=71 GPa，泊松比ν=0.3，密度ρs=2 700 kg/m3，板結(jié)構(gòu)厚度為0.001 m?？諝庥蚵晫W(xué)參數(shù)為：密度ρf=1.21 kg/m3，聲速c=343 m/s。施加在矩形平面板中心點(diǎn)+Z方向的簡諧激勵力幅值為1 N。

本文計算結(jié)構(gòu)-聲場耦合系統(tǒng)問題均通過Matlab程序?qū)崿F(xiàn)。因?yàn)橛邢拊椒ňW(wǎng)格越密結(jié)果越準(zhǔn)確，所以參考值通過FEM方法計算較密的網(wǎng)格得到。將平板結(jié)構(gòu)域和聲場域耦合平面按不同網(wǎng)格密度劃分成節(jié)點(diǎn)均勻分布的四邊形網(wǎng)格，不同耦合平面網(wǎng)格密度下的結(jié)構(gòu)域和聲域單元數(shù)目與單元類型如表1所示。

表1 六面體模型的網(wǎng)格劃分

首先應(yīng)用FE-LSPIM計算矩形彈性板的特征頻率值，并與SFEM和FEM進(jìn)行對比研究。將矩形彈性板離散成網(wǎng)格密度為8×8的網(wǎng)格模型，用FE-LSPIM、SFEM和FEM分別求得矩形彈性板的前14階固有頻率值，結(jié)果如表2所示。表2中參考值由文獻(xiàn)[1]得到。從表2中可以看出，SFEM所得結(jié)果略優(yōu)于FEM所得結(jié)果；除了第一階特征頻率外，F(xiàn)E-LSPIM所得結(jié)果明顯優(yōu)于SFEM和FEM所得結(jié)果。

表2 應(yīng)用FE-LSPIM、SFEM和FEM計算的板固有頻率值

表3 應(yīng)用FE-LSPIM、FEM計算的空腔固有頻率值

表3所示為用FE-LSPIM計算得到聲腔的前14階非零固有頻率值，耦合平面網(wǎng)格密度為8×8。表中同時給出了用FEM在相同網(wǎng)格模型下計算得到的聲腔的前14階非零固有頻率值以及參考值，參考值由文獻(xiàn)[1]得到。從表3中可以看出，隨著特征頻率階次的增加，由FE-LSPIM和FEM所得結(jié)果的誤差呈變大趨勢，但FE-LSPIM所得結(jié)果明顯比FEM所得結(jié)果更靠近參考值。

表4 應(yīng)用FE-LSPIM/FE-LSPIM、SFEM/FEM和FEM/FEM計算的耦合六面體聲腔固有頻率值

表4所示為用FE-LSPIM/FE-LSPIM計算得到耦合聲腔的前14階非零固有頻率值，耦合平面網(wǎng)格密度為8×8。為了便于對比分析，表中同時給出了用SFEM/FEM、FEM/FEM在相同網(wǎng)格模型下計算得到的耦合聲腔的前14階非零固有頻率值以及參考值，參考值由文獻(xiàn)[1]得到。從表4中可以看出：

(1) 除了第一階特征頻率外，應(yīng)用FE-LSPIM/FE-LSPIM計算得到耦合聲腔的特征頻率比應(yīng)用SFEM/FEM、FEM/FEM計算所得結(jié)果更靠近參考值；

(2) 隨著特征頻率階次的增加，F(xiàn)E-LSPIM/FE-LSPIM、SFEM/FEM和FEM/FEM的計算結(jié)果更加遠(yuǎn)離參考值，計算誤差整體呈變大趨勢，但是在大多數(shù)情況下，F(xiàn)E-LSPIM/FE-LSPIM計算結(jié)果的誤差明顯小于SFEM/FEM和FEM/FEM計算結(jié)果的誤差，相對能取得較高精度的結(jié)果。

為分析單元尺寸對FE-LSPIM/FE-LSPIM耦合聲場計算結(jié)果的影響，將平板結(jié)構(gòu)域和聲場域耦合平面按不同網(wǎng)格密度劃分成節(jié)點(diǎn)均勻分布的四邊形網(wǎng)格，網(wǎng)格密度分別?。?×6、8×8、10×10、12×12、16×16，相應(yīng)的結(jié)構(gòu)域單元數(shù)和聲域單元數(shù)如表1所示。為評價FE-LSPIM/FE-LSPIM的計算效果，同時給出FEM/FEM和SFEM/FEM的聲壓計算結(jié)果進(jìn)行對比研究，參考值通過NASTRAN計算耦合平面網(wǎng)格密度為50×50的模型得到。表5和表6分別表示激振頻率為150Hz及230Hz時應(yīng)用FE-LSPIM/FE-LSPIM、SFEM/FEM和FEM/FEM計算的域點(diǎn)1的聲壓值。

表5 激振頻率為150 Hz時域點(diǎn)1聲壓

表6 激振頻率為230 Hz時域點(diǎn)1聲壓

綜合表5和表6可知：隨著網(wǎng)格密度的增加(即單元尺寸的減少)，F(xiàn)E-LSPIM/FE-LSPIM、SFEM/FEM和FEM/FEM的計算結(jié)果都逼近參考結(jié)果，說明FE-LSPIM/FE-LSPIM、SFEM/FEM和FEM/FEM都是收斂的，且FE-LSPIM/FE-LSPIM比SFEM/FEM和FEM/FEM收斂更快。

為了進(jìn)一步評價FE-LSPIM/FE-LSPIM 分析結(jié)構(gòu)-聲場耦合問題的效果，本文計算了域點(diǎn)1和2的聲壓頻率響應(yīng)，計算頻率范圍為20～320 Hz，并給出了SFEM/ FEM和FEM/ FEM的計算結(jié)果作為對比，參考值通過FEM/ FEM計算耦合平面網(wǎng)格密度為16×16的模型得到。

圖6和圖7分別表示耦合平面網(wǎng)格密度為8×8時三種方法所計算的域點(diǎn)1和2的聲壓頻率響應(yīng)曲線。從圖中可以看出：通過SFEM/FEM、FEM/FEM和FE-LSPIM/FE-LSPIM所得的聲壓頻率響應(yīng)曲線趨勢與參考值相同，由于網(wǎng)格尺寸較大，計算結(jié)果均與參考值相差較大，但與SFEM/FEM和FEM/FEM相比，F(xiàn)E-LSPIM/FE-LSPIM的計算結(jié)果更接近參考值。這表明FE-LSPIM/FE-LSPIM比SFEM/ FEM和FEM/ FEM更精確，在同等粗糙的網(wǎng)格模型下計算結(jié)果更好。

圖8為不同網(wǎng)格密度下用FE-LSPIM/FE-LSPIM、SFEM/FEM和FEM/FEM計算得到的域點(diǎn)1聲壓頻率響應(yīng)曲線，計算頻率范圍為20～320 Hz。FE-LSPIM/FE-LSPIM計算的耦合平面網(wǎng)格密度為8×8； SFEM/FEM和FEM/FEM計算的耦合平面網(wǎng)格密度分別取10×10和12×12。從圖中可以看出，隨著網(wǎng)格密度的增加，計算結(jié)果越準(zhǔn)確，聲壓頻率響應(yīng)曲線整體向左移動靠近參考值。FE-LSPIM/FE-LSPIM在網(wǎng)格密度為8×8時的計算結(jié)果與SFEM/FEM和FEM/FEM在網(wǎng)格密度為12×12時的計算結(jié)果相接近，從而進(jìn)一步說明了對粗糙網(wǎng)格模型，F(xiàn)E-LSPIM/FE-LSPIM比SFEM/ FEM和FEM/ FEM的計算精度更高。

圖6 域點(diǎn)1聲壓頻率響應(yīng)曲線

4 結(jié) 論

本文將FE-LSPIM推廣用于Mindlin-Reissner板結(jié)構(gòu)動力學(xué)和三維聲場分析，提出用于板結(jié)構(gòu)-聲場耦合問題分析的FE-LSPIM/FE-LSPIM，推導(dǎo)了FE-LSPIM/FE-LSPIM分析板結(jié)構(gòu)-聲場耦合問題的計算公式。并以一六面體聲場-結(jié)構(gòu)耦合模型為研究對象進(jìn)行分析，研究結(jié)果表明：FE-LSPIM/FE-LSPIM能很好地應(yīng)用于板結(jié)構(gòu)- 聲場耦合分析中，其計算結(jié)果比SFEM/ FEM和FEM/ FEM收斂快，并且對單元網(wǎng)格尺寸要求比SFEM/ FEM和FEM/ FEM低，在計算較大網(wǎng)格尺寸模型時，能得到比SFEM/ FEM和FEM/ FEM更高的計算精度。因此，用FE-LSPIM/FE-LSPIM分析板結(jié)構(gòu)-聲場耦合問題可以減小計算規(guī)模、節(jié)省計算時間，具有良好的工程應(yīng)用前景。

[1]He Z C, Liu G R, Zhong Z H, et al. A coupled edge-/face-based smoothed finite element method for structural-acoustic problems[J]. Applied Acoustics, 2010, 71: 955-964.

[2]He Z C, Liu G R, Zhong Z H, et al. Coupled analysis of 3D structural-acoustic problems using the edge-based smoothed finite element method/finite element method[J]. Finite Elements in Analysis and Design, 2010, 46: 1114-1121.

[3]Atalla N, Bernhard R J. Review of numerical solution for low-frequency structural-acoustic problems. Applied Acoustics, 1994(43):271-249.

[4]姚凌云,于德介,臧獻(xiàn)國. 殼結(jié)構(gòu)聲場耦合分析的光滑有限元-有限元法[J]. 中國機(jī)械工程，2010，21(15)：115-120.

YAO Ling-yun, YU De-jie, ZANG Xian-guo. Analysis of shell structural-acoustic coupling based on smoothed finite element and finite element method[J]. Chinese Mechanical Engineering, 2010，21(15)：115-120.

[5]Melenk J M, Babu ska I. The partition of unity finite element method：basic theory and applications[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4)：289-314.

[6]Wenterodt C, Estorff O V. Dispersion analysis of the meshfree radial point interpolation method for the Helmholtz equation[J]. International Journal for Numerical Methods in Engineering, 2009,77：1670-1689.

[7]Rajendran S, Zhang B R. A “FE-Meshfree” QUAD4 Element for Free-vibration Analysis[J]. Computer Methods in Applied Mechanicals and Engineering, 2007, 197: 3595-3604.

[8]Rajendran S, Zhang B R. A “FE-Meshfree” QUAD4 Element Based on Partition of Unity[J]. Computer Methods in Applied Mechanicals and Engineering, 2007, 197: 128-147.

[9]姚凌云,于德介,臧獻(xiàn)國. 聲學(xué)數(shù)值計算的有限元-最小二乘點(diǎn)插值法[J]. 中國機(jī)械工程，2010，21(23)：2816-2820.

YAO Ling-yun, YU De-jie, ZANG Xian-guo. A finite element-least square point interpolation method for acoustic numerical computation[J]. Chinese Mechanical Engineering, 2010，21(23)：2816-2820.