一種基于譜元法的層合板SEA參數計算方法

袁金龍，陳海波，王 昆，操小龍

(1.中國科學院 材料力學行為和設計重點實驗室，中國科學技術大學近代力學系，合肥 230027；2.北京機電工程研究所，北京 100074)

復合材料具有強度高、重量輕、抗疲勞等特點，已經廣泛應用于航空航天等領域。飛行器在飛行過程中承受氣動熱、寬頻噪聲、振動等復雜惡劣環境，這些環境因素對其性能和可靠性有著重要的影響。而針對構件高頻域的動態響應，有限元和邊界元需要用很小的網格尺寸來描述結構特征長度，使得數值分析的計算量和存儲量變得非常巨大[1]，甚至不能實現。統計能量分析方法將復雜系統劃分為不同的模態群，并從統計意義上把大系統分解為若干個便于分析的獨立子系統，對共振模態響應進行平均，從而得到光滑準確的能量值[2]。

統計能量分析方法的預示精度在很大程度上取決于統計能量分析參數的準確性。對于三明治板這類夾芯層合板結構的研究，國內外學者已經做了很多工作。 Renji等[3]推導了考慮芯層橫向剪切剛度影響的夾層板四階振動控制方程，并由控制方程在波數空間積分得到了正交各向異性復合材料夾層板的模態密度公式，最后通過點導納方法測得的模態密度試驗值對理論公式進行了驗證。Zhou等[4]由夾層板六階振動控制方程推導了模態密度的計算公式，計算了蜂窩夾層板的傳聲損失。然而該研究只針對各向同性面板的夾層板，忽略了芯層面內剛度對夾層板模態密度的影響。Han等[5]考慮面板彎曲剛度和芯層面內剛度, 提出了一種改進的經典夾層板理論,能更準確地計算三明治板的模態密度。

對于一般的層合板結構，三明治板理論不再適用，這時需要采用層合板理論或數值計算方法。經典層合板理論(CLPT)是基于Kirchhoff薄板假設,認為板變形后的中面法線依然垂直于板中面，適用于薄板結構低頻分析[6]；一階剪切板理論(FSDT)考慮橫向剪切和轉動慣量影響，認為板變形后其中面法線仍為直線但不與中面垂直，適用于更高頻率的分析[7]；而更高階的剪切板理論控制方程復雜，波數求解困難。Chronopoulos等[8]利用波有限元方法計算了正交各向異性復合材料層合板的頻散特性，并利用波數和群速度數值求解了層合板的模態密度，但是計算量巨大。

譜元法(SFEM)是一種典型的半解析方法，是在有限元方法的基礎上提出的，它綜合了有限元方法處理復雜結構問題的靈活性以及收斂速度快和計算精度高的特點，對結構模擬波傳播特性十分有效，是提取色散曲線的一種重要方法。譜元法假設平面簡諧波沿著波導軸方向傳播，并基于有限元方法對波導截面內的位移場進行離散，利用類似有限元的方法進行單元剛度矩陣和單元質量矩陣的組裝，從而獲得一個標準的波數特征值方程。Zak等[9]將裂紋近似為剛度線性變化的彈簧，采用二維譜元法研究了含裂紋的各向同性板中的彈性波傳播特性，并進而提出了一種損傷識別的方法。Peng等[10]采用三維譜元法模擬彈性波在各向同性板中的傳播特性，并與采用二維譜元法的模擬結果進行對比，通過控制加載方式得到了蘭姆波的不同模態。Kudela 等[11]基于Mindlin 板理論采用二維譜單元法研究了彈性波在復合材料層合板中的傳播，并研究了鋪層角以及纖維體積含量對層合板彈性波傳播特性的影響。Datta等[12]采用譜元法研究了彈性波在對稱交錯層疊層合板結構的傳播特性。Shorter[13]采用一維線性譜單元對線彈性各向同性層合板進行波導分析,推導了層合板總內損耗因子的計算公式。徐超等[14]推導了一種任意形狀的時域譜單元，研究了高頻導波在二維功能梯度材料結構中的傳播特性。

與傳統的有限元方法相比，譜元方法只需要對波導截面位移場進行離散，降低了模型的自由度，從而極大地提高了計算效率。而且由于在波的傳播方向采用精確解，從而避免了有限元離散時產生的數值誤差，進一步提高了計算結果的準確性。但將該法應用于統計能量參數獲取的研究還未見報道。本文采用二次譜單元對一般層合板厚度方向進行離散，建立波數的特征根方程，并在波數求解時采用皮爾遜相關系數對其對應的模態進行分類，進而計算一般層合板的統計能量分析參數。通過對薄板、三明治夾芯板以及厚板的數值模擬分析，驗證了該方法的可行性和有效性，這對一般層合板的統計能量分析具有工程指導意義。

1 譜元法(SFEM)

1.1 單元位移場

平面波在平板中傳播時，高頻段波傳播對板面內尺寸如長、寬不敏感，因此在厚度方向可采用一維單元離散,由波的傳播特性，單元位移場可表示為：

(1)

式中:ω為圓頻率，k為波數，φ為方位角，u(z)、v(z)和w(z)為厚度方向z坐標處波傳播的振幅。長度為Li的單元中，位移振幅通過一維形函數插值可表示為：

(2)

圖1 一維網格

p=u,v,w

(3)

N(z)為單元的形函數，文獻[13]采用的是二節點單元，為了進一步提高計算精度和效率，本文采用三節點二次單元：

(4)

1.2 單元質量矩陣與剛度矩陣

單元的應變矩陣為：

(5)

其中：

(6)

將式(1)和式(2)代入式(5)中得

(7)

其中：

B1=LzN,z,B2=cos(φ)LxN+cos(φ)LyN

(8)

N,z為形函數矩陣對坐標z的偏導，單元內的應力為：

(9)

式中:C為模量矩陣，S為柔量矩陣，對于正交各向異性材料的獨立彈性常數為9個，即三個拉壓彈性模量、三個剪切模量和三個主泊松比：

在層合板坐標系下，模量矩陣為：

Cx-y=T(θ)-1CRT(θ)R-1

(10)

式中：m=cos(θ)，n=sin(θ)，θ為材料主軸方向與參考方向的夾角。層合板動能為：

(11)

層合板勢能為：

(12)

H代表共軛轉置。將式(7)代入式(12)可得：

U=

(13)

由哈密頓變分原理可得

(14)

其中：

將所有單元的剛度矩陣和質量矩陣組裝起來即可獲得總體的剛度矩陣和質量矩陣：

[K1(φ)+K2(φ)k+K3(φ)k2-ω2M]q0=0

(15)

2 統計能量分析

2.1 模態密度求解

當方位角和頻率給定時，方程式(15)就變成k的二次特征根問題，二次特征根問題又可以進一步轉化成下列線性特征值問題[13]：

(16)

其中：

q1=kq0;C2=K3;

C1=K2;C0=K1-ω2M

(17)

特征值的實部和虛部的符號代表給定波的傳播方向，特征向量代表給定頻率和方位角的橫截面波的形狀。在低頻段，振動波主要為彎曲波、剪切波和拉伸波，當頻率繼續升高時將出現高階振動波，因此需要對波的傳播類型進行分類。層合板理論和三明治板理論都是通過求解高次方程來獲得波數，存在波數分類困難的問題。譜元法在求解特征值的同時得到了特征向量，從而可用它來對波數進行分類。文獻[8]采用特征向量余弦相似度的平方進行波數分類。余弦相似度指通過計算兩個向量的夾角余弦值來評估它們的相似度，它適用于數據稀疏的情況，但易受向量平移的影響。而皮爾遜相關系數是利用兩個變量之間的協方差和標準差的商來衡量兩個變量的相關性，相比于余弦相似度，它對數據進行了中心化，具有平移不變性和尺度不變性的特點，且不受維度缺失的影響[15]。

根據模態相似原則，本文采用皮爾遜相關系數的絕對值來進行振型分類：

(18)

式中:Qm和Qn為相鄰方位角或者相鄰頻率點的特征向量，E表示取均值。對振型相似的模態，P的絕對值接近1；在實際分類時，先取較大值，再取較小值，直到篩選出最相似的模態。利用方程(16)和(18)即可獲得各階模態對應的波數。

根據統計能量分析方法中模態密度的定義，在波數空間中，圓頻率ω下的模態數N(ω)可由圓頻率ω曲線包含的面積獲得：

(19)

式中:a和b為矩形板的長和寬。由統計能量法模態密度的定義有：

(20)

2.2 結構與聲腔的耦合

耦合損耗因子是子系統間耦合作用大小的一種度量。結構與聲腔相互作用時結構的損耗因子ηrad與輻射阻抗Rrad有關：

(21)

式中:Ms為板結構的總質量,對于四邊簡支板的輻射阻抗可由下式獲得[15]:

(22)

式中:As為板結構的面積，fc為板結構的臨界頻率，c為空氣中聲音的傳播速度，ρa為空氣密度。P為板邊緣總長度。

2.3 內損耗因子

子系統的動力學特性與系統的剛度、質量和阻尼有關，阻尼在頻響曲線共振區起主導作用，因此阻尼特性對研究結構的動力學響應和聲傳遞特性具有重要意義。在統計能量分析方法中，一般用內損耗因子來表示系統的阻尼損耗特性。內損因子可定義為：

(23)

式中:Pdiss為時間平均的功率耗散。總功率損耗因子可以由各層的應變能獲得：

(24)

式中:n是層合板的總層數，將式(24)代入(23)得：

(25)

式中:Ki是層合板第i層的組裝剛度，即第i層對總體剛度矩陣的貢獻量。考慮到高頻振動時結構動能與勢能的互等關系，將第k階傳播波的振型向量代入式(5)，該階波的內損耗因子可表示為：

(26)

式中:Pk是第k階傳播波的特征向量。由式(26)可知，總內損耗因子與各層的勢能有關，單層的勢能越大，其對總內損耗因子的貢獻越大。

2.4 傳聲損失

(27)

圖2 傳聲模型

其中

(28)

式中：Vi表示聲源室和接收室的體積，ηi為內損耗因子，聲場的內損耗因子可以由混響時間T3獲得，ni為模態密度，聲場的模態密度可由下式獲得：

(29)

式中:Ai為聲場的總表面積，Si聲場的邊緣的總長度。

3 數值算例

3.1 各向同性薄板

對于厚度為h的各向同性薄板，其彎曲、剪切和拉伸的波數理論值可由下式給出：

(30)

以厚度為1 mm的均質鋁板A為研究對象，鋁板的密度為2 700 kg/m3,彈性模量為71 GPa，泊松比為0.329 6。圖3為采用不同譜單元SFEM計算結果與式(30)理論值相比的誤差。由圖3可知，二次單元只需一個單元就能獲得較準確的結果，而線性單元用了4個還未達到相應的精度，由此可見二次單元比線性單元的效率更高。圖4為分別采用式(30)與采用二次單元的SFEM方法的波數計算結果對比，雖然SFEM只用了1個單元，二者幾乎一致。本文后面的算例均采用二次單元計算。由圖5可知，對于各向同性材料，波數的分布與方位角無關。

圖3 板A彎曲波數SFEM與理論值誤差(方位角為0°)

圖4 板A波數(方位角為0°)

圖5 板A波數與方位角的關系(SFEM)

3.2 正交各向異性薄板

正交各向異性薄板B為碳-碳(Carbon)復合材料板，材料參數如表1所示。由圖6可知，當厚度比較小時，SFEM的波數計算結果與層合板理論CLPT和FSDT幾乎一致。對于正交各向異性材料，波數的分布與方位角有關，如圖7所示，此時板B在靠近x軸時波數較大。

表1 材料力學參數

圖6 板B波數(方位角為0°)

圖7 板B波數與方位角的關系(SFEM)

至此，通過各向同性板和正交各向異性板的波數計算結果對比，驗證了SFEM計算波數的有效性。

3.4 三明治板

三明治板是一類特殊的層合板結構，一般由兩塊較薄、模量較大的面板和輕質、模量較小的厚夾芯層構成。以文獻[4]中的蜂窩板C為研究對象，板C為夾芯層為6.35 mm的蜂窩板結構，面板和夾芯材料的力學性質，如表1所示。

采用SFEM分析板C時，面板劃分1個單元，夾芯層劃分5個單元，其他參數與文獻[4]一致。計算出波數后，模態密度可由式(20)計算。如圖8所示，在較低頻段,SFEM與文獻[4]和文獻[5]的三明治板理論計算結果差距較小，都比實驗值略小；在較高頻段，SFEM計算結果相比于文獻[4]和文獻[5]與實驗值吻合更好。由SFEM獲得的臨界頻率為1 680 Hz，通過式(22)計算輻射阻抗如圖9所示，相比于文獻[4]中的數值解，由SFEM計算的輻射阻抗與實驗值一致性更好。如圖10所示，由SFEM計算的傳聲損失在低頻段和高頻段與實驗值都比較接近，在臨界頻率處二者有一些偏差，總體結果要比文獻[4]中的好。通過三明治板的模態密、輻射效率和傳聲損失的計算結果對比，進一步驗證了本文所提SFEM計算SEA參數的準確性。

圖8 板C的模態密度

圖9 板C的輻射阻抗

圖10 板C的傳聲損失

3.3 一般層合板

層合板D、E和F為五層厚度相同的層合板結構，長和寬均為1 m，每層的材料均為Carbon，鋪層角度為[0/90/0/90/0]。單層的厚度分別為0.1 mm、0.5 mm和1 mm。分別采用CLPT、FSDT和SFEM計算層合板的模態密度，如圖11～13所示，CLPT和FSDT相對于SFEM的偏差如圖14所示。CLPT隨著厚度和頻率升高，模態密度計算偏差越大，FSDT在厚度小時與SFEM吻合得較好，當層合板厚度變大時，偏差增大。CLPT、FSDT和SFEM的計算時間對比如表2所示，FSDT的計算時間比CLPT略大，而SFEM的計算時間比CLPT多9倍左右。影響SFEM計算效率的主要因素是總體剛度矩陣和質量矩陣的維度以及波形分類引起的額外計算量，較細的單元劃分會增加矩陣的維度但可以提高計算精度，因此與普通有限元一樣，實際計算中要尋求符合計算精度的較稀疏網格進行計算。另外，采用SFEM計算SEA參數的原理與層合板理論計算SEA參數是一樣的，都是基于波數的方法，因此其計算時間在求解過程中并不會頻率升高顯著增加。下一步我們將對如何提高SFEM的計算效率做專題研究。

圖11 板D的模態密度

圖12 板E的模態密度

圖13 板F的模態密度

圖14 層合板理論相對于SFEM的偏差

表2 計算時間對比

三明治理論雖然也能較為準確的計算模態密度，但是它僅適用于特殊夾芯結構；層合板理論在厚度和頻率增大時，計算偏差會變大；三明治板理論和層合板理論都是求解高次方程，難于進行模態分類，而SFEM在求解特征值問題時能利用皮爾遜相關系數準確地進行模態分類，適用于一般層合板結構的分析，并且對網格要求低，計算量小，精確度高。對于厚板，還可以適當增加單元數目來提高SFEM計算結果的準確性。

3.4 汽車玻璃層合板

板G是一種典型汽車玻璃(Glass)層合板，由兩塊厚度為4 mm的玻璃板和0.75 mm厚的PVB材料組成，其如圖15所示[13]。

圖15 板G結構示意圖

Glass的材料為各向同性材料，楊氏模量為62 GPa，泊松比為0.24，密度為2 300 kg/m3；PVB為黏彈性材料，泊松比為0.499，其剪切模量與頻率有關，其關系如圖2所示，PVB材料的內損耗因子也是與頻率有關，其關系如圖17所示，Glass的內損耗因子取0[13]。

圖16 PVB材料剪切模量

板G的三層各劃分為一個單元，利用SFEM計算獲得各模態的波數和特征向量，根據式(26)獲得板G的總內損耗因子如圖18所示，基于SFEM的總內損耗因子計算結果與文獻[13]和VA one軟件的結果幾乎是一致的，驗證了SFEM計算總內損耗因子的準確性。彎曲模態的內損耗因子隨頻率的升高先增大后減小，量級為10-1左右；而剪切模態和拉伸模態的損耗因子持續增大，量級為10-5～10-4左右。黏彈性層PVB對彎曲模態損耗因子影響較大，而對剪切模態和拉伸模態的損耗因子影響較小。

圖17 PVB材料內損耗因子

圖18 板G內損耗因子

4 總 結

本文基于二次譜單元建立層合板的特征方程，在波數求解時，采用了皮爾遜相關系數來衡量模態的相似度，根據相似度準確地進行波數分類，進而計算層合板的統計能量分析參數。通過數值模擬驗證了SFEM方法的通用性與有效性。本文的SFEM在處理復雜鋪層、任意厚度和材料分布的層合平板方面，比傳統層合板理論精度更高、適用范圍更廣，且其計算量遠遠小于傳統有限元方法。但要解決曲板、加筋板等復雜層合板結構的建模問題，需要發展高維的譜單元。本文基于SFEM的SEA參數計算方法對飛行器層合板結構的高頻建模有較重要的工程實際意義。