文潤(rùn)旭， 袁麗蕓*， 劉嘉明， 劉麗萍

(1.廣西科技大學(xué)機(jī)械與交通工程學(xué)院， 柳州 545006； 2.廣西汽車(chē)零部件與整車(chē)技術(shù)重點(diǎn)實(shí)驗(yàn)室， 柳州 545006)

彈性空腔結(jié)構(gòu)應(yīng)用在工程中的各個(gè)領(lǐng)域，如船舶的船艙、飛機(jī)的機(jī)艙等，可以進(jìn)一步簡(jiǎn)化為彈性空腔的模型。對(duì)于彈性空腔的振動(dòng)與噪聲問(wèn)題一直是迫切需要研究人員解決的，提高彈性空腔聲振特性性能是科技前瞻性的研究。有限元法是現(xiàn)如今工程分析中最廣泛的數(shù)值方法，并隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展有限元法受到科技工作者的高度重視并得以應(yīng)用，它具有靈活性、高效性、邊界適應(yīng)性好等特點(diǎn)，并且能夠適應(yīng)于多物理場(chǎng)的數(shù)值分析，如聲固耦合、磁固耦合等。有限元數(shù)值方法的主要思想是通過(guò)把問(wèn)題域離散為若干單元，并將它們的邊界上的結(jié)點(diǎn)相互連接，用近似函數(shù)表示求解域的未知變量，使用變分原理建立求解未知的代數(shù)方程。其中，二維板結(jié)構(gòu)的三角形單元，由于邊界適應(yīng)性強(qiáng)，位移函數(shù)構(gòu)造簡(jiǎn)單，計(jì)算規(guī)模小，得到了廣泛的應(yīng)用。顧太平等[1]基于直角坐標(biāo)系的三角形離散網(wǎng)格，采用面積坐標(biāo)構(gòu)造了薄板三角形單元內(nèi)部任意一點(diǎn)的撓度位移插值函數(shù)。薛俊好等[2]在文獻(xiàn)[1]的基礎(chǔ)上，進(jìn)一步給出了該三角形板單元的應(yīng)變矩陣計(jì)算公式。然而，上述兩篇文獻(xiàn)均缺少進(jìn)一步的算例分析。由于插值型函數(shù)直接由原始網(wǎng)格生成，當(dāng)單元網(wǎng)格畸變時(shí)，其精度將難以保證，角形板在進(jìn)行彈性板的動(dòng)力學(xué)問(wèn)題分析中，具有較好的邊界適應(yīng)性，然而需要進(jìn)一步處理單元畸變引起的精度較差問(wèn)題。

針對(duì)三維空腔問(wèn)題，常采用的方法包括有限元法、邊界元法、光滑有限元法等，應(yīng)用廣泛的三維單元包括四面體單元和六面體單元。閔祥斗等[3]利用有限元方法對(duì)高速列車(chē)車(chē)室進(jìn)行聲固耦合建模，經(jīng)過(guò)Ansys軟件仿真預(yù)測(cè)中低段噪聲。楊年炯等[4]針對(duì)駕駛室轟鳴聲問(wèn)題使用有限元方法對(duì)引氣管仿真繼而優(yōu)化，其中在引氣管不平整面采用三角形網(wǎng)格劃分，避免了網(wǎng)格畸變產(chǎn)生的計(jì)算誤差。Du等[5]對(duì)彈性板為邊界的三維空腔進(jìn)行了聲振耦合分析，并且使用Rayleigh-Ritz法，基于結(jié)構(gòu)-聲學(xué)耦合系統(tǒng)的能量表達(dá)式，確定了其低階模態(tài)頻率。白潔[6]采用Galerkin法對(duì)三維聲場(chǎng)進(jìn)行離散，生成四面體單元，車(chē)身采用四邊形單元，構(gòu)造了車(chē)身三維彈性空腔的有限元耦合模型，采用Virtual Lab獲得了耦合系統(tǒng)的固有頻率。結(jié)果表明在駕駛艙不規(guī)則區(qū)域邊界采用四邊形網(wǎng)格離散，會(huì)產(chǎn)生奇異點(diǎn)，從而影響計(jì)算精度，甚至發(fā)生不收斂現(xiàn)象。姚昊萍等[7]構(gòu)建了由2塊彈性板和4塊剛性板構(gòu)成的聲-固耦合空腔模型，利用了瑞利-里茨方法求解耦合系統(tǒng)的特征值，使用紐曼邊界條件，分析了空腔內(nèi)的聲壓響應(yīng)，通過(guò)算例結(jié)果驗(yàn)證了此模型解析解與實(shí)驗(yàn)值一致性。戴騰飛[8]應(yīng)用光滑有限元法對(duì)柴油機(jī)缸蓋內(nèi)外聲學(xué)建模，并驗(yàn)證此方法的可行性。需要指出的是，上述研究中聲場(chǎng)區(qū)域均采用有限元法構(gòu)造聲學(xué)剛度矩陣，由于有限元本身的局限性，其剛度矩陣偏硬，因此，計(jì)算精度有待進(jìn)一步的提高。而光滑有限元法在處理聲學(xué)問(wèn)題上具有對(duì)剛度“軟化”的作用，逐漸得到了學(xué)者們的青睞。郭小斌等[9]總結(jié)了光滑有限元法的類(lèi)型和優(yōu)點(diǎn)，此方法能夠軟化剛度矩陣、處理體積鎖定問(wèn)題，應(yīng)用也極其廣泛，其中在聲學(xué)的應(yīng)用中具有較高的收斂性和精度。Li等[10]使用光滑有限元法對(duì)結(jié)構(gòu)-聲學(xué)耦合問(wèn)題分析，由于板和彈性空腔都是采用光滑有限元法進(jìn)行模型的建立，所以具有對(duì)有限元“過(guò)硬”情況軟化的效果，因此比有限元法精度更高，但其彈性板所采用的三角形單元面內(nèi)位移采用了一階剪切變形理論，與橫向位移相關(guān)，仍具有一定的局限性。張智瑯等[11]基于光滑有限元-解耦有限粒子法(smooth finite element method, decoupled finite particle method, SFEM-DFPM)模擬流固耦合問(wèn)題，能夠解決傳統(tǒng)有限元“過(guò)剛”的問(wèn)題，此耦合方法能有效處理自由液面的流固耦合問(wèn)題。

綜上所述，目前對(duì)彈性空腔的研究多采用有限元法或光滑有限元法進(jìn)行，其中聲場(chǎng)區(qū)域采用光滑技術(shù)處理可以將聲剛度矩陣進(jìn)行“軟化”，從而提高計(jì)算精度；而彈性板采用4節(jié)點(diǎn)3自由度單元或3節(jié)點(diǎn)3自由的板單元進(jìn)行離散，面內(nèi)位移與橫向位移相關(guān)，計(jì)算精度有待進(jìn)一步的提高。最近，劉嘉明等[12]用邊光滑有限元進(jìn)行二維約束層阻尼( constrained layer damping ，CLD)全敷設(shè)矩形腔的聲振分析，并與有限元法比對(duì)，具有較好的計(jì)算精度，其研究思路可進(jìn)一步拓展到三維情形?；诖?，擬采用3節(jié)點(diǎn)5自由度的三角形單元對(duì)薄板進(jìn)行離散處理，以便得到邊界適應(yīng)性好、精度較高的薄板有限元模型；再對(duì)彈性空腔的聲剛度矩陣進(jìn)行光滑處理，以軟化聲場(chǎng)剛度矩陣；然后考慮彈性板和空腔的耦合條件，從而得到三維彈性空腔的一種新的聲振耦合模型，并進(jìn)一步對(duì)該類(lèi)模型的聲振特性進(jìn)行數(shù)值探討。

1 三維彈性空腔聲固耦合模型

如圖1所示，三維彈性空腔是由聲腔域和柔性板組成，其中頂部采用柔性板，其他五個(gè)壁是剛性的，對(duì)于此類(lèi)彈性空腔使用一種新方式組合建模，分別對(duì)板進(jìn)行三角形單元有限元法建模、聲腔域進(jìn)行光滑有限元建模。下面分別闡述板有限元模型、聲腔域光滑有限元模型，以及三維彈性空腔耦合模型。

圖1 三維彈性空腔耦合模型Fig.1 Three-dimensional elastic cavity coupling model

1.1 薄板有限元模型的建立

針對(duì)上述板單元，引入薄板的中面位移向量，即

Δ=[u,v,w,θx,θy]T

(1)

根據(jù)文獻(xiàn)[1]，采用邊界適應(yīng)性較好三角形單元對(duì)板結(jié)構(gòu)進(jìn)行離散。引入面積坐標(biāo)L1、L2、L3，則單元內(nèi)任意點(diǎn)(xi,yi)的撓度可寫(xiě)為

(2)

式(2)中：形函數(shù)具體表達(dá)為

(3)

(4)

(5)

同時(shí)，考慮三節(jié)點(diǎn)三角形單元任意點(diǎn)(x,y)獨(dú)立于橫向位移的面內(nèi)位移，即

(6)

由式(2)、式(6)引入形函數(shù)矩陣Ns，即

(7)

為了盡量減少網(wǎng)格畸變所帶來(lái)的計(jì)算誤差，將采用等參單元進(jìn)行坐標(biāo)映射，將物理坐標(biāo)映射到自然坐標(biāo)系中進(jìn)行計(jì)算。此時(shí)，可將單元內(nèi)部任意點(diǎn)的坐標(biāo)用節(jié)點(diǎn)坐標(biāo)和形函數(shù)表達(dá)為

L3=1-ξ-η

(8)

用變分法可得板結(jié)構(gòu)的動(dòng)力平衡方程為

(9)

式(9)中：us為板節(jié)點(diǎn)位移向量；Ms為板質(zhì)量矩陣；Ks為板的剛度矩陣；Fs為結(jié)構(gòu)的節(jié)點(diǎn)等效外部激勵(lì)向量；Ffs為流體在流固耦合邊界作用在結(jié)構(gòu)上的節(jié)點(diǎn)等效聲壓向量。其表達(dá)式為

(10)

(11)

(12)

(13)

1.2 聲腔域光滑有限元模型的建立

通過(guò)標(biāo)準(zhǔn)的伽遼金方法，得到聲腔域Helmholtz方程的加權(quán)殘值弱形式，即

(14)

式(14)中：kω為聲場(chǎng)波數(shù)；ρc為聲場(chǎng)介質(zhì)密度；ΓN為邊界條件。

(15)

采用光滑梯度技術(shù)，聲學(xué)問(wèn)題的光滑聲壓梯度可表示為

(16)

式(16)中：Bfi(x,y)表示光滑梯度矩陣；M為局部光滑域的節(jié)點(diǎn)總數(shù)。

將式(15)和式(16)代入式(14)中，運(yùn)用變分原理，并通過(guò)單元集成分析，可以得到聲學(xué)域光滑有限元?jiǎng)恿W(xué)方程為

(17)

式(17)中：Kf表示聲學(xué)光滑剛度矩陣；Mf表示聲學(xué)質(zhì)量矩陣；Ff為空腔在邊界處的振動(dòng)所引起的激勵(lì)力向量，其表達(dá)式[12]為

(18)

(19)

(20)

1.3 三維彈性空腔聲固耦合方程的建立

在聲腔域板的耦合處，板的振動(dòng)法向速度會(huì)對(duì)聲腔產(chǎn)生耦合力，并結(jié)合式(20)，該耦合作用力可以表示為

(21)

結(jié)合式(9)、式(17)可得彈性空腔耦合方程為

(22)

對(duì)于自由振動(dòng)問(wèn)題，式(22)的右邊激勵(lì)力向量為0，代入位移和聲壓邊界條件(一般為0)后，不難獲得該類(lèi)問(wèn)題的特征方程為

|K-λM|=0

(23)

式(23)中：K和M是代入邊界條件后，劃掉邊界向量為0時(shí)對(duì)應(yīng)的各行和各列元素后，整理所得的剛度矩陣和質(zhì)量矩陣。進(jìn)一步求解得可得耦合系統(tǒng)的固有頻率和振型。

2 數(shù)值算例

2.1 柔性板的模態(tài)分析

為了驗(yàn)證文中三角形板單元有限元法的準(zhǔn)確性及良好的適應(yīng)性，此算例對(duì)矩形薄板進(jìn)行模態(tài)分析，幾何和材料參數(shù)取自參考文獻(xiàn)[10]，該矩形薄板為柔性鋁板，長(zhǎng)a=0.4 m、寬b=0.32 m、厚h=0.001 m、彈性模量E=71 GPa、泊松比υ=0.3、密度ρ=2.7×103kg/m3，約束條件采用四邊簡(jiǎn)支。

用本文數(shù)值解與文獻(xiàn)[10]中的有限元方法(finite element mothod，F(xiàn)EM)、光滑有限元方法(smoothed finite element mothod，SFEM)數(shù)值進(jìn)行固有頻率及誤差大小的對(duì)比，結(jié)果見(jiàn)表1所示。圖2、圖3則給出了ANSYS和本文所用方法所得簡(jiǎn)支板的前兩階振型圖。其中本文方法所得振型圖根據(jù)MATLAB編程所得模態(tài)振型矩陣，代入Tecplot軟件繪制所得。

從表1可知，在低階時(shí)，本文數(shù)值算法的計(jì)算結(jié)果誤差稍高于SFEM計(jì)算結(jié)果，但在中高階時(shí)，其誤差低于SFEM計(jì)算結(jié)果；本文方法整體上計(jì)算結(jié)果誤差小于文獻(xiàn)[10]FEM和ANSYS的計(jì)算結(jié)果；此外，圖2、圖3給出了ANSYS和本文所用方法所得簡(jiǎn)支板的前兩階振型對(duì)比圖，極度相似。綜上，本文數(shù)值算法結(jié)果接近SFEM結(jié)果，并且本文的計(jì)算精度比傳統(tǒng)有限元精度高。因此，使用本文數(shù)值算法對(duì)板進(jìn)行模態(tài)分析是精確的，且計(jì)算規(guī)模小。

圖2 矩形薄板第一階模態(tài)振型Fig.2 The first mode of the rectangular plate

圖3 矩形薄板第二階模態(tài)振型Fig.3 The second mode of the rectangular plate

表1 四邊簡(jiǎn)支矩形薄板前十階固有頻率對(duì)比Table 1 Comparison of the first ten natural frequencies of the triangular thin plate

2.2 三維彈性空腔聲振特性分析

為了評(píng)估本文方法對(duì)三角形單元的薄板理論推導(dǎo)的合理性，進(jìn)一步考查本文中使用的柔性板有限元建模、聲腔域光滑有限元建模，再進(jìn)行新耦合方式的合理性，及三角形單元與四面體單元的適應(yīng)性效果，對(duì)此分析三維彈性空腔聲固耦合的聲振性能。

本例中，為了跟文獻(xiàn)[10]作對(duì)比，封閉矩形空腔的尺寸參數(shù)選為：長(zhǎng)為0.4 m、寬為0.32 m、高為0.36 m。板材為鋁制，同算例1中柔性薄板的幾何參數(shù)和物理參數(shù)一致。柔性薄板由三角形網(wǎng)格建模，該簡(jiǎn)支板的邊界條件位于邊緣?？涨粌?nèi)的三維空氣是由四面體網(wǎng)格模擬，其空氣的物理參數(shù)：密度ρ=1.225 kg/m3、聲速c=340 m/s，計(jì)算得出了由柔性板組合三彈性空腔前十階固有頻率結(jié)果見(jiàn)表2所示。其中，本文模型的網(wǎng)格尺寸為0.04 m,ANSYS有限元網(wǎng)格尺寸為0.02 m。

從表2可知，三維彈性空腔模型固有頻率的對(duì)比中SFEM的計(jì)算結(jié)果誤差小于ANSYS的計(jì)算結(jié)果；本文數(shù)值算法比SFEM計(jì)算精度更高，并且更接近參考解，在中高階模態(tài)中計(jì)算結(jié)果較明顯；進(jìn)一步分析圖4、圖5可得，ANSYS和本文所用方法所得三維彈性空腔前兩階振型基本相同。

圖4 三維彈性空腔模型的一階模態(tài)振型Fig.4 The first mode of the three-dimensional elastic cavity model

圖5 三維彈性空腔模型的二階模態(tài)振型Fig.5 The second mode of the three-dimensional elastic cavity model

表2 三維彈性空腔模型前十階固有頻率對(duì)比Table 2 Comparison of the first ten natural frequencies of the three-dimensional elastic cavity model

綜上，驗(yàn)證了FEM與SFEM耦合算法是正確的，并且計(jì)算效果較好，說(shuō)明了在保證計(jì)算規(guī)模較小的前提下，本文方法仍具有較好的計(jì)算精度。

3 結(jié)論

本文中三維彈性空腔聲固耦合新模型是分別采用三角形單元對(duì)板進(jìn)行有限元建模、光滑有限元法對(duì)聲腔域進(jìn)行光滑處理，再考慮二者的邊界耦合條件進(jìn)一步建立的。通過(guò)與文獻(xiàn)[10]和ANSYS仿真結(jié)果對(duì)比可得出以下結(jié)論。

(1)采用3節(jié)點(diǎn)5自由度的三角形單元對(duì)薄板進(jìn)行離散處理計(jì)算規(guī)模小，精度足夠。

(2)基于FEM與SFEM耦合算法的聲固耦合模型，計(jì)算精度高，規(guī)模也相對(duì)較小。

(3)由于采用的三角形板單元具有良好的邊界適應(yīng)性和較小的計(jì)算規(guī)模，聲腔域的光滑處理也軟化了聲場(chǎng)域的剛度矩陣，可保證較好的計(jì)算精度，本文中所建立的聲振耦合模型可進(jìn)一步地應(yīng)用于實(shí)際工程中形狀更復(fù)雜彈性空腔結(jié)構(gòu)的聲振特性分析中。