點(diǎn)陣夾芯圓柱殼自由振動分析的波有限元方法

韓少燕, 高汝鑫

(1. 西安交通大學(xué)城市學(xué)院 機(jī)械工程學(xué)院, 西安 710018;2. 北京理工大學(xué) 先進(jìn)結(jié)構(gòu)技術(shù)研究院, 北京 100081;3. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實驗室, 遼寧 大連 116024)

0 引 言

點(diǎn)陣夾芯結(jié)構(gòu)被認(rèn)為是最具應(yīng)用前景的下一代輕量化多功能結(jié)構(gòu)形式之一,其有著傳統(tǒng)材料結(jié)構(gòu)所不具備的力學(xué)性能,例如較高的比剛度、比強(qiáng)度、能量吸收能力和性能可設(shè)計等[1]．常見的點(diǎn)陣夾芯結(jié)構(gòu)單胞的拓?fù)錁?gòu)型有金字塔型[2]、四面體型[2]、Kagome型[2]、三維全三角型[3]和八面體型[4]等．點(diǎn)陣夾芯結(jié)構(gòu),如點(diǎn)陣夾芯梁、板和殼結(jié)構(gòu),因其上述的優(yōu)異性能已成功應(yīng)用于航空航天領(lǐng)域[5]．國內(nèi)外學(xué)者針對點(diǎn)陣夾芯結(jié)構(gòu)的力學(xué)和熱學(xué)性能開展了眾多研究[6-7]．根據(jù)點(diǎn)陣夾芯結(jié)構(gòu)的構(gòu)成特點(diǎn),其多用于抗沖吸能結(jié)構(gòu)設(shè)計[8]．點(diǎn)陣夾芯結(jié)構(gòu)應(yīng)用于航空航天、船舶等機(jī)械結(jié)構(gòu)時,通常會遭遇復(fù)雜的振動環(huán)境,所以研究點(diǎn)陣夾芯結(jié)構(gòu)的振動特性是十分必要的．近年來,學(xué)者們針對點(diǎn)陣夾芯結(jié)構(gòu)的線性和非線性振動開展了相關(guān)研究工作[9-11]．

波在周期結(jié)構(gòu)中傳播時滿足Bloch定理[12],點(diǎn)陣夾芯結(jié)構(gòu)通常具有周期性,故波有限元法為周期結(jié)構(gòu)的振動分析提供了一種快速求解方案．然而,波有限元方法往往面臨較多的數(shù)值問題,通常可分為三類[13]： ① 當(dāng)頻率較高時出現(xiàn)的有限元離散誤差和色散誤差[14]; ② 當(dāng)胞元長度與變形波長相當(dāng)時出現(xiàn)的空間離散和周期結(jié)構(gòu)效應(yīng)[12,15]; ③ 傳遞矩陣的條件數(shù)過大,矩陣病態(tài)導(dǎo)致本征值和本征向量無法求解[16]．對于前兩種數(shù)值問題,可以通過減小單元尺寸加密網(wǎng)格來避免;而傳遞矩陣數(shù)值病態(tài)的問題,就只能通過建立高精度的數(shù)值求解方法來解決．為此,學(xué)者們針對傳遞矩陣本征值和本征向量的求解提出了眾多方法：Zhong和Williams[16]應(yīng)用結(jié)構(gòu)力學(xué)和最優(yōu)控制理論的類比關(guān)系,提出了一種直接求解方法,該方法將傳遞矩陣的本征值的求解轉(zhuǎn)化為兩個非病態(tài)矩陣的廣義特征值的求解;傳遞矩陣的本征值問題可以轉(zhuǎn)化為回文二次特征值問題,Huang等[17]針對回文二次特征值問題提出了保結(jié)構(gòu)算法;Wang等[18]進(jìn)一步擴(kuò)展了Zhong和Williams的直接求解法,突破了原方法只能處理單根的限制,使得該方法可以處理多重根問題．但是上述傳遞矩陣特征值的求解方法多用來處理較為簡單的胞元或者自由度規(guī)模較小的傳遞矩陣,當(dāng)矩陣規(guī)模較大時,上述方法會出現(xiàn)累積誤差．鑒于上述原因,目前波有限元方法通常將較為規(guī)則的一維或二維單元作為單胞求解均質(zhì)結(jié)構(gòu)的振動分析問題[19-21],或用來處理胞元構(gòu)型較為簡單的周期(點(diǎn)陣)結(jié)構(gòu)的振動分析問題[22-24]．

本文基于波傳播理論和有限元方法,針對點(diǎn)陣夾芯圓柱殼的自由振動分析,發(fā)展了波有限元方法．首先,根據(jù)自由波在圓柱殼軸向和周向的傳播規(guī)律,基于波有限元方法建立了胞元的振動控制方程．該振動控制方程雖然建立在胞元自由度上,但可描述整個點(diǎn)陣夾芯圓柱殼的自由振動行為．其次,對于胞元振動控制方程,基于Neumann級數(shù)推導(dǎo)了約束剛度矩陣的逆矩陣的顯式表達(dá)式,不僅提高了求逆效率,而且將振動控制方程中的頻率分離出來,從而將點(diǎn)陣夾芯圓柱殼的固有頻率求解轉(zhuǎn)化為胞元邊界自由度規(guī)模的二次特征值問題．最后,根據(jù)駐波形成的條件,給出了圓柱殼軸向和周向的波傳播參數(shù),進(jìn)而求解二次特征值問題即可得到點(diǎn)陣夾芯圓柱殼的固有頻率．相比于全尺寸有限元模型,本文的方法可顯著降低計算模型的自由度規(guī)模,從而顯著降低計算成本,提高計算效率,且沒有數(shù)值病態(tài)的問題．另外,推導(dǎo)過程不指定胞元的結(jié)構(gòu)形式,且其剛度矩陣和質(zhì)量矩陣可以方便地從任意的商業(yè)有限元軟件中導(dǎo)出,便于工程應(yīng)用．

1 胞元的波有限元控制方程

考慮點(diǎn)陣夾芯圓柱殼的自由振動,其中一個胞元的有限元模型如圖1所示,該胞元的有限元振動方程可以寫為

Dq=f,

(1)

其中

D=K-ω2M

(2)

為胞元的動剛度矩陣,q為胞元的節(jié)點(diǎn)自由度向量,f為胞元的節(jié)點(diǎn)力向量,K和M分別為胞元的剛度矩陣和質(zhì)量矩陣,ω為點(diǎn)陣夾芯圓柱殼自由振動的圓頻率．

根據(jù)圖1對節(jié)點(diǎn)自由度的劃分,可以將式(1)寫為分塊的形式：

(3)

其中

(4)

圖1 點(diǎn)陣夾芯圓柱殼胞元的有限元模型

胞元內(nèi)部節(jié)點(diǎn)作用力為零,即fi=0,展開式(3)第2行并移項可得

(5)

將式(5)代入到式(3)第1行可得

(6)

根據(jù)Bloch-Floquet定理,胞元相鄰邊界上的節(jié)點(diǎn)自由度向量具有以下關(guān)系：

q2=λxq1,qR=λxqL,q4=λθq1,qt=λθqb,q3=λθλxq1,

(7)

其中λx和λθ分別為胞元在軸向和周向的波傳播參數(shù)．

根據(jù)式(7)和(4),可以將qe表示為

qe=ΛQ,

(8)

其中

(9)

(10)

將式(8)代入式(6)得到

(11)

根據(jù)胞元力的平衡條件[25]有

(12)

(13)

(14)

將式(12)—(14)寫為矩陣形式有

ΛHfe=0．

(15)

將式(11)代入式(15)有

(16)

其中

(17)

式(16)即為點(diǎn)陣夾芯圓柱殼胞元的波有限元控制方程．

2 約束動剛度矩陣的顯式求逆

(18)

利用Neumann級數(shù)[26]將式(18)近似為

(19)

將式(19)代入到式(17),可得

工藝過程：將漢生膠分散在甘油里面，再倒入適量燒開的去離子水，攪拌均勻，加入B相剩余組分，加熱攪拌至 85 ℃；混合A相并加熱至85 ℃；混合A,B兩相，均質(zhì)5 min,并保溫攪拌30 min；冷卻至60 ℃加入C相（用少量水分散），繼續(xù)冷卻攪拌至45 ℃以下加入D相組分，繼續(xù)攪拌冷卻至室溫。

(20)

其中

(21)

(22)

(23)

將式(20)代入到式(16),可得

(μ2A+μB+C)Q=0,

(24)

其中

(25)

(26)

(27)

μ=ω2．

(28)

取z=[μQTQT]T,則式(24)變換為

(29)

從式(29)容易看出,其中的系數(shù)矩陣是軸向波傳播參數(shù)λx和周向波傳播參數(shù)λθ的函數(shù),均與圓頻率ω?zé)o關(guān)．對于給定的軸向波傳播參數(shù)λx和周向波傳播參數(shù)λθ,點(diǎn)陣夾芯圓柱殼的頻率求解轉(zhuǎn)換為求解式(29)的廣義特征值問題．另外,從式(29)還可以看出,本文方法的求解自由度規(guī)模為通常不大于一個胞元的自由度,不需要對整個點(diǎn)陣夾芯圓柱殼進(jìn)行建模,且對胞元的具體構(gòu)型沒有特別要求,即胞元的結(jié)構(gòu)矩陣可由現(xiàn)有的商業(yè)軟件導(dǎo)出．

3 點(diǎn)陣夾芯圓柱殼中的波傳播分析

根據(jù)結(jié)構(gòu)中自由波的傳播與自振模式的對應(yīng)關(guān)系：當(dāng)波傳播參數(shù)取適當(dāng)?shù)闹?使得波在結(jié)構(gòu)中反射疊加后形成駐波時,結(jié)構(gòu)達(dá)到了某一階自振模式[27-28]．圓柱殼自由振動涉及周向和軸向兩個方向的波傳播．首先考慮周向,考慮周期性邊界條件,容易知道圓柱殼周向的波數(shù)為整數(shù),則胞元在周向的波傳播參數(shù)為

λθ=einθcell,n=0,±1,±2,…,

(30)

其中θcell為胞元在周向掠過的弧度．

其次考慮胞元軸向波傳播參數(shù)的求解,本文利用自由波在梁中的傳播分析來求解圓柱殼軸向的波傳播參數(shù)．根據(jù)駐波形成的條件,波在軸向傳播一個周期需要滿足以下條件：

2kL+(φL+φR)=2mπ,m=0,±1,±2,…,

(31)

其中L為圓柱殼的長度,k為軸向半波數(shù),φL為入射波和反射波在圓柱殼左端邊界上的相位差,φR為入射波和反射波在圓柱殼右端邊界上的相位差．對于簡支邊界條件,入射波和反射波的相位差為-π,對于固支或自由邊界條件,入射波和反射波的相位差為-π/2[27-28]．

將式(31)移項可得軸向半波數(shù),進(jìn)而軸向的波傳播參數(shù)為

(32)

其中Lcell為胞元沿軸向的長度．

需要注意的是,根據(jù)文獻(xiàn)[27-28],自由波在梁中傳播的通解可以由4個分量疊加而成,分別為e-ikx,eikx,e-kx,ekx,其中k的正負(fù)號分別代表著波的傳播方向．可以看出,前兩個分量為傳播波分量,后兩個分量為近場波(near field wave),也稱倏逝波(evanescent field wave)．隨著傳播距離的增加,近場波呈指數(shù)衰減．考慮不同的邊界條件,在兩端簡支的邊界條件下,近場波為零[27-28];而在其他邊界條件,如固支、自由邊界條件等情況下,近場波不為零．然而,需要指出的是,由于近場波衰減很快,且頻率越高,衰減越快,所以當(dāng)頻率較高、波數(shù)較大時,近場波的影響可以忽略不計．鑒于上述分析,本文方法在推導(dǎo)過程中沒有考慮近場波的影響,故本文方法在兩端簡支邊界條件下,具有較高的計算精度,而在其他非簡支邊界條件下,精度受近場波的影響：對于較低階的自振頻率,本文方法的計算精度較差;對于高階自振頻率的求解,波數(shù)的增加使得近場波影響變小,計算精度會越來越高．

最后,點(diǎn)陣夾芯圓柱殼的模態(tài)可以通過入射波和反射波疊加得到[27-28],例如,對于兩端簡支的點(diǎn)陣夾芯圓柱殼,其模態(tài)φ可以表示為

(33)

綜上,本文方法求解點(diǎn)陣夾芯圓柱殼自振頻率和模態(tài)的步驟如下：

① 根據(jù)式(30)和(32),求得點(diǎn)陣夾芯圓柱殼在軸向和周向的波傳播參數(shù)λx和λθ;

② 將兩個方向的波傳播參數(shù)代入式(9),并結(jié)合式(25)—(27)可得式(29)中的系數(shù)矩陣A,B和C;

③ 求解式(29)的廣義特征值問題即可得到廣義特征值μ和廣義特征向量z;

④ 結(jié)合廣義特征值μ和式(28)即可得到自振頻率ω,結(jié)合廣義特征向量z、式(8)和(10)即可得到qe; ⑤ 將qe代入式(5)可以得到qi,進(jìn)一步利用式(33)即可得到點(diǎn)陣夾芯圓柱殼的模態(tài)．

4 數(shù) 值 算 例

本節(jié)利用一個點(diǎn)陣夾芯圓柱殼來驗證本文方法的正確性和高效性,其周向有300個胞元,軸向有250個胞元,如圖2(a)所示．胞元為金字塔型,如圖2(b)所示,其前后面板分別為半徑為0.99 m和1.01 m的圓柱面,軸向長度均為20 mm,周向掠過的角度均為1.2°,厚度均為2 mm,四根連接梁的截面為圓形,連接梁半徑為1 mm．胞元材料的質(zhì)量密度為2 711 kg/m3,彈性模量為68.98 GPa,Poisson比為0.33．點(diǎn)陣夾芯圓柱殼兩端的邊界條件依次考慮為簡支-簡支(SS)、固支-固支(CC)、自由-自由(FF)、固支-簡支(CS)和固支-自由(CF),分別利用本文方法和有限元方法對其進(jìn)行自由振動分析．

(a) 點(diǎn)陣夾芯圓柱殼 (b) 金字塔型胞元 (a) A lattice core sandwich cylindrical shell(b) A pyramidal lattice truss core element圖2 點(diǎn)陣夾芯圓柱殼及其胞元示意圖

4．1 兩端簡支

首先考慮兩端簡支邊界條件,胞元的上下面板均使用4個相同大小的殼單元建模,芯子中的4根連桿均使用3個梁單元建模,胞元有限元模型共有156個自由度;全尺寸有限元模型由胞元有限元模型在兩個方向上陣列得到,共有7 207 200個自由度．

本文方法中使用的胞元剛度矩陣和質(zhì)量矩陣從Nastran中導(dǎo)出,并使用MATLAB編寫計算程序,使用CPU為8核*3.60 GHz,內(nèi)存為32 GB的桌面級電腦,單核計算一個工況所花費(fèi)時間大概為1 s．全尺寸有限元模型使用Nastran建模,計算在CPU為80核*2.3 GHz、內(nèi)存為512 GB的高性能工作站上進(jìn)行,使用18核計算一個工況(前100階模態(tài))所花費(fèi)時間大概為40 min．需要指出的是,計算時間除了受計算規(guī)模、硬件配置影響以外,還受其他因素的影響,故此處兩者的計算時間為多次計算時間取的平均值．以上的計算時間對比驗證了本文方法的高效性．

表1給出了本文方法和全尺寸有限元模型計算得到的頻率結(jié)果,由于本文重點(diǎn)關(guān)注圓柱殼的呼吸模態(tài),所以表中沒有給出周向波數(shù)為0或1對應(yīng)的固有頻率．可以看出,兩端簡支邊界條件下,本文方法計算的點(diǎn)陣夾芯圓柱殼的固有頻率與全尺寸有限元模型的結(jié)果吻合很好,驗證了本文方法的有效性．

表1 兩端簡支邊界條件下(SS)點(diǎn)陣夾芯圓柱殼的自由振動頻率對比

圖3給出了本文方法和有限元方法計算得到的點(diǎn)陣夾芯圓柱殼的模態(tài)形狀(m=4,n=3)對比,可以看出兩種方法計算的模態(tài)形狀吻合較好．

4．2 其他邊界條件

本小節(jié)考慮了其他邊界條件下的點(diǎn)陣夾芯圓柱殼的自由振動問題,依次為固支-固支(CC)、自由-自由(FF)、固支-簡支(CS)和固支-自由(CF)．表2—5給出了上述4種邊界條件下,本文方法和有限元方法計算得到的點(diǎn)陣夾芯圓柱殼的固有頻率對比．從表中數(shù)據(jù)可以看出,在上述4種邊界條件下,本文方法的計算結(jié)果與有限元方法的計算結(jié)果誤差比在兩端簡支邊界條件下的結(jié)果誤差大,這是由上述4種邊界條件下近場波的影響造成的．隨著波數(shù)的增加近場波迅速衰減,其對結(jié)果的影響越來越小,表2—5中的誤差數(shù)據(jù)驗證了這一點(diǎn)．另外,從表中還可以看出,雖然本文方法在上述4種邊界條件下的計算結(jié)果與有限元計算結(jié)果存在誤差,但是誤差大小整體在可接受范圍內(nèi)．

(a) 有限元方法 (b) 本文方法(a) The finite element method (b) The present method圖3 本文方法和有限元方法計算的模態(tài)對比(m=4, n=3)

表2 兩端固支邊界條件下(CC)點(diǎn)陣夾芯圓柱殼的自由振動頻率對比

表3 兩端自由邊界條件下(FF)點(diǎn)陣夾芯圓柱殼的自由振動頻率對比

表4 一端固支一端簡支邊界條件下(CS)點(diǎn)陣夾芯圓柱殼的自由振動頻率對比

表5 一端固支一端自由邊界條件下(CF)點(diǎn)陣夾芯圓柱殼的自由振動頻率對比

5 結(jié) 論

本文針對點(diǎn)陣夾芯圓柱殼的自由振動分析發(fā)展了波有限元方法．首先,基于二維波有限元方法建立了點(diǎn)陣夾芯圓柱殼單個胞元的控制方程,該控制方程的自由度規(guī)模顯著小于全尺寸有限元模型;然后,利用Neumann級數(shù)推導(dǎo)了約束剛度矩陣的顯式表達(dá)式,使得控制方程中的頻率分離出來,從而將固有頻率求解轉(zhuǎn)化為胞元的二次特征值問題;最后,考慮駐波的形成條件求得圓柱殼軸向和周向的參數(shù),得到了點(diǎn)陣夾芯圓柱殼的固有頻率．?dāng)?shù)值算例表明：本文方法對兩端簡支點(diǎn)陣夾芯圓柱殼的自由振動具有很高的預(yù)測精度,而對于其他邊界條件,受近場波的影響預(yù)測精度略差,然而隨著波數(shù)的增加,近場波快速衰減,本文方法的預(yù)測精度會越來越高．

致謝本文作者衷心感謝大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實驗室開放課題(GZ21110)對本文的資助．