陳秀濤 卿光輝

(中國(guó)民航大學(xué)航空工程學(xué)院，天津300300)

通常情況下，帶有加強(qiáng)筋的板殼結(jié)構(gòu)具有剛度大、耗材少和重量輕等特點(diǎn)，因此被廣泛應(yīng)用于各類實(shí)際工程中。但由于加筋結(jié)構(gòu)的振動(dòng)行為對(duì)整體結(jié)構(gòu)的穩(wěn)定性影響十分顯著，所以有關(guān)加筋結(jié)構(gòu)的振動(dòng)特性分析一直是國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)問題。

Mukherjee 等[1]和 Mukhopadhyay 等[2]對(duì)加筋結(jié)構(gòu)振動(dòng)分析的早期文獻(xiàn)進(jìn)行了綜述：由于加筋結(jié)構(gòu)的復(fù)雜性，基于位移有限元模型的數(shù)值分析方法最為普遍。Trkmen 等[3]提出了暴露于沖擊波的加筋板的位移有限元分析模型，并進(jìn)行了實(shí)驗(yàn)研究。Zhao等[4]使用能量方法，研究了簡(jiǎn)單支撐的交叉加筋多層圓柱殼的自由振動(dòng)。Rikards 等[5]發(fā)展了一種精度較高的三角形單元，研究了加筋復(fù)合材料層合殼的自由振動(dòng)。Guo 等[6]對(duì)加筋層合板結(jié)構(gòu)提出了分層有限元模型，并對(duì)這類結(jié)構(gòu)進(jìn)行了屈曲分析。游翔宇等[7]基于光滑有限元法對(duì)加筋板的自由振動(dòng)進(jìn)行了研究。

值得說明的是，一般位移有限元法是采用有限的節(jié)點(diǎn)代替了真實(shí)結(jié)構(gòu)中無窮多的點(diǎn)，結(jié)果導(dǎo)致有限元模型比真實(shí)模型的剛度大，故而固有頻率結(jié)果大于真實(shí)解。從理論上講，由于合理地引入非協(xié)調(diào)位移項(xiàng)，使其近似多項(xiàng)式的階次完備，致使非協(xié)調(diào)元的柔性增加。系統(tǒng)控制方程的柔性也相應(yīng)地增加，因而采用非協(xié)調(diào)元的模型可以得到精度更高的數(shù)值結(jié)果。然而，即使采用收斂性較好的非協(xié)調(diào)位移元分析加筋板的振動(dòng)特性問題，其相應(yīng)的有限元模型的剛度也是偏硬的，因而在有限元網(wǎng)格較稀疏的情況下，固有頻率結(jié)果通常是大于真實(shí)解的。

近些年來，一些研究人員應(yīng)用Hamilton 正則方程半解析法分析加筋層合結(jié)構(gòu)[8]或厚度不連續(xù)結(jié)構(gòu)的振動(dòng)特性問題[9]。Hamilton 正則方程半解析法分析板殼結(jié)構(gòu)的顯著優(yōu)勢(shì)是允許材料和幾何特性沿厚度方向改變，并可以方便地分析層合板殼結(jié)構(gòu)類的振動(dòng)問題。通常情況下，基于Hamilton 正則方程半解析法的模型，數(shù)值結(jié)果相對(duì)比較準(zhǔn)確。但是，該方法有一個(gè)明顯的缺點(diǎn)，即由于層合板殼類結(jié)構(gòu)的每層在面域內(nèi)對(duì)位移變量和面外應(yīng)力變量進(jìn)行了有限元離散，因而每一層的單元組裝后導(dǎo)致規(guī)模較大的數(shù)值矩陣。在形成結(jié)構(gòu)整體的控制方程前需要對(duì)每一層的數(shù)值矩陣分別進(jìn)行關(guān)于厚度變量的指數(shù)矩陣運(yùn)算，這一過程的數(shù)值運(yùn)算成本高。另一方面，完成了每一層矩陣指數(shù)運(yùn)算后，為保證不同材料層的層間位移變量和應(yīng)力變量的連續(xù)性，又要對(duì)各層的控制方程進(jìn)行乘積運(yùn)算，這一計(jì)算步驟也需要一定的時(shí)間?？傊琀amilton 正則方程半解析法面對(duì)工程中龐大的結(jié)構(gòu)問題，單元離散的數(shù)量大，數(shù)值運(yùn)算相對(duì)于常規(guī)的有限元位移法來說，矩陣指數(shù)運(yùn)算量大，時(shí)間消耗大。

傳統(tǒng)的混合有限元法有很多優(yōu)點(diǎn)，例如，通常只需C0連續(xù)的插值多項(xiàng)式來近似表達(dá)待求的場(chǎng)量，位移和應(yīng)力計(jì)算結(jié)果精度高。最近，Qing 等[10]結(jié)合最小勢(shì)能原理和H–R 變分原理建立了分析靜力學(xué)問題的非協(xié)調(diào)廣義混合元。文獻(xiàn)[11] 擴(kuò)展了這種非協(xié)調(diào)廣義混合元的應(yīng)用范疇。

本文在文獻(xiàn) [10-11] 的基礎(chǔ)上，根據(jù) Hamilton原理，首先建立考慮結(jié)構(gòu)振動(dòng)特性的廣義混合變分原理，然后建立了非協(xié)調(diào)廣義混合元的動(dòng)力學(xué)模型，用于加筋板的振動(dòng)特性分析。

1 位移元的振動(dòng)特征方程

通常情況下，以位移元為基礎(chǔ)的無阻尼自由振動(dòng)方程[12]為

式中，M表示結(jié)構(gòu)的質(zhì)量矩陣，為單元的一致質(zhì)量矩陣，N表示形函數(shù)矩陣,ρ表示材料的密度，V表示連續(xù)體的體積，K表示結(jié)構(gòu)的剛度矩陣，是單元的剛度矩陣，B為應(yīng)變矩陣，D是材料的彈性矩陣。

設(shè)q(t) 為待求的解，考慮簡(jiǎn)諧運(yùn)動(dòng)問題時(shí)，位移解的表達(dá)形式為

式中，Φ表示固有振型，ω表示對(duì)應(yīng)固有振型的頻率。

將式 (2) 代入式 (1) 可得到無阻尼自由振動(dòng)的特征方程為

對(duì)式(3) 的求解方法較多，例如，逆迭代法、子空間迭代法和行列式搜索法等。

2 廣義混合元的振動(dòng)特征方程

根據(jù)廣義混合變分原理[13]，可設(shè)結(jié)構(gòu)的混合能(包括應(yīng)變能和余能) 為

式中，σ表示應(yīng)力向量，C表示材料的剛度系數(shù)矩陣，?表示微分算子，u表示位移向量。

處于自由振動(dòng)狀態(tài)下的結(jié)構(gòu)的動(dòng)能為

式中，ω表示固有頻率，ρ表示材料密度。

考慮動(dòng)力學(xué)問題的Hamilton 變分原理為

因此含有固有頻率ω參數(shù)的廣義混合變分原理可表示為

式中，參數(shù)α的取值范圍為 0α1。

參考文獻(xiàn)[10-11] 的推導(dǎo)過程，由式(7) 可得到六面體非協(xié)調(diào)廣義混合單元的自由振動(dòng)特征方程為

式中，h為單元柔度矩陣，g為單元杠桿矩陣，pe表示單元節(jié)點(diǎn)的應(yīng)力參數(shù)，Φe表示單元固有振型參數(shù)，

各式中的形函數(shù)矩陣N和非協(xié)調(diào)項(xiàng)形函數(shù)矩陣Nr表達(dá)形式可參見文獻(xiàn)[10-11]。

由式(8) 可導(dǎo)出

將式 (9) 代入式 (10) 中有

其中，剛度矩陣κ=αk+(1?α)gTh?1g。

式(11)表明六面體非協(xié)調(diào)廣義混合元的剛度矩陣κ與柔度矩陣g相關(guān)，參數(shù)α的取值范圍為0～1, 并且不同的取值可得到不同的剛度矩陣。若式 (7) 中α= 1，則式 (11) 中的κ與位移元下的剛度矩陣一致，同為基于最小勢(shì)能原理的非協(xié)調(diào)位移元的剛度矩陣。α1 時(shí)，其不同的取值可以調(diào)節(jié)剛度矩陣柔性。文獻(xiàn)[14] 推導(dǎo)含參數(shù)α的廣義混合元模型時(shí)，小變形線性彈性問題是最基本的假設(shè)。另一方面，假設(shè)含參數(shù)α的廣義混合元模型的應(yīng)變能與精確的應(yīng)變能相等。基于這一假設(shè)，文獻(xiàn) [14] 導(dǎo)出了參數(shù)α= 0.75 的最優(yōu)值。因此，就小變形線性彈性問題而言，廣義混合元模型中的α= 0.75 對(duì)計(jì)算結(jié)果的最主要影響是數(shù)值解更接近問題的真實(shí)解。

對(duì)式(11)求和可得整體有限元模型的振動(dòng)特性方程為

關(guān)于非協(xié)調(diào)廣義混合元模型施加邊界條件的方法：最終的振動(dòng)特性方程 (12) 中只含有位移變量。因此方程(12) 只需考慮位移邊界條件的引入，其引入方法與位移元法中的方法相同。

3 實(shí)例分析

以下數(shù)值分析均采用式(12)對(duì)算例進(jìn)行數(shù)值分析，其中參數(shù)α=0.75。此外，在關(guān)于加筋板的算例中，為了提高計(jì)算效率，對(duì)式(12) 中的單元質(zhì)量矩陣采用集中質(zhì)量矩陣的形式。

3.1 層合板的自由振動(dòng)特性分析

如圖 1 所示的三層板，每層均為正交異性材料，坐標(biāo)軸沿彈性主方向，上、下兩層材料相同。令分別表示第一和第二層的C11值，ρ1和ρ2分別是第一層和第二層材料的密度。當(dāng)δ= 1，γ= 1 時(shí)，該板即為一單層正交異性板，板長(zhǎng)度a和寬度b，a=b，厚度h，h/a=0.1，h1/h=h3/h=0.1，h2/h=0.8。邊界條件如圖2 所示，為四邊簡(jiǎn)支，即當(dāng)x=0 和x=a時(shí)，u=w=0，當(dāng)y=0 和y=b時(shí)，v=w=0。剛度系數(shù)比為

圖1 層合板尺寸圖

圖2 四邊簡(jiǎn)支約束

以下通過改變網(wǎng)格密度來討論非協(xié)調(diào)廣義混合元的收斂性。

在γ= 1 和δ= 1 的模型下，在x方向和y方向劃分同等單元(從2 個(gè)單元遞增到12 個(gè)單元)，在厚度方向z劃分10 個(gè)單元，保持不變。

第一階固有頻率的收斂曲線如圖3 所示。從圖中可以明顯看出，隨著網(wǎng)格密度的增加，通過非協(xié)調(diào)位移元和非協(xié)調(diào)廣義混合元計(jì)算的頻率結(jié)果都逐漸趨于平穩(wěn)。當(dāng)網(wǎng)格密度為8×8×10 時(shí)，兩種方法的結(jié)果都逼近精確解，但非協(xié)調(diào)廣義混合元的結(jié)果更加精確。

圖3 不同網(wǎng)格密度下的1 階固有頻率

以下根據(jù)δ，γ的不同取值，分別求其一階固有頻率。網(wǎng)格劃分情況如圖 4 所示。兩種有限元模型的固有頻率結(jié)果見表 1，且兩種模型的 1 階振型圖一致，如圖5 所示。

比較表1 中的數(shù)據(jù)，不難看出，非協(xié)調(diào)位移元的結(jié)果高于文獻(xiàn)解。這是因?yàn)槲灰圃Ｐ偷膭偠认鄬?duì)于實(shí)際的模型偏硬。在網(wǎng)格模型相同的情況下，因?yàn)榉菂f(xié)調(diào)廣義混合元模型的剛度更接近真實(shí)結(jié)構(gòu)，所以本文的非協(xié)調(diào)廣義混合元的結(jié)果優(yōu)于非協(xié)調(diào)位移元。這一結(jié)論與文獻(xiàn)[14] 中的相關(guān)論述相同。

圖4 層合板網(wǎng)格圖

√/C(2)11 )本文方法0.048 06 1.33 0.057 65 1.06 150.077 41相對(duì)誤差/% — 1.51 0.35

圖5 1 階振型圖

3.2 單加強(qiáng)筋偏心加筋板的自由振動(dòng)特性分析

帶有單加強(qiáng)筋的偏心加筋板，尺寸圖如圖6 所示，材料參數(shù)為：彈性模量E= 68.7 GPa，泊松比ν= 0.29，密度ρ= 2823 kg/m3。幾何參數(shù)：a=b= 0.203 2 m，t= 0.001 371 6 m，w=0.006 35 m，h=0.012 7 m，邊界條件為四邊固支，如圖 7 所示，當(dāng)x=0 和x=a時(shí)，u=v=w=0，當(dāng)y=0 和y=a時(shí)，u=v=w=0。對(duì)該結(jié)構(gòu)的振動(dòng)特性進(jìn)行分析。

圖6 單加強(qiáng)筋偏心加筋板尺寸圖

將板部分劃分成 27×27×1 個(gè)六面體單元，加強(qiáng)筋部分劃分成1×27×8 個(gè)六面體單元，網(wǎng)格劃分情況如圖8 所示。兩種有限元法的前4 階固有頻率結(jié)果見表2。前4 階固有振型見圖9，其中非協(xié)調(diào)位移元的結(jié)果由有限元商業(yè)軟件Abaqus 得到。

圖7 四邊固支約束

圖8 單加強(qiáng)筋偏心加筋板網(wǎng)格劃分圖

表2 單加強(qiáng)筋偏心加筋板前 4 階固有頻率(單位：Hz)

由表2 中數(shù)據(jù)可以看出，非協(xié)調(diào)廣義混合元模型的固有頻率值低于非協(xié)調(diào)位移元模型的結(jié)果數(shù)值，距離文獻(xiàn)解更加接近，結(jié)果精度更高。

比較圖9 中兩種方法的前4 階固有振型不難看出，非協(xié)調(diào)廣義混合元模型和非協(xié)調(diào)位移元模型的固有振型幾乎是一致的，這進(jìn)一步驗(yàn)證了本文方法的正確性。

3.3 雙加強(qiáng)筋偏心加筋板的自由振動(dòng)特性分析

雙加強(qiáng)筋偏心加筋層合板[8]，尺寸如圖10 所示，尺寸參數(shù)：a= 1.2 m，l= 0.4 m，t= 0.01 m，t1=0.02 m，w= 0.04 m，h= 0.08 m。板部分總共三層，每層均為正交異性材料，上下兩個(gè)面層材料相同，假設(shè)面層的密度ρ1=1600 kg/m3，芯層的密度ρ1/ρ2=2，加強(qiáng)筋的材料參數(shù)與面層的一致。邊界條件為四邊固支，如圖11 所示，當(dāng)x=0 和x=a時(shí)，u=v=w= 0，當(dāng)y=0和y=a時(shí)，u=v=w= 0。對(duì)該結(jié)構(gòu)進(jìn)行振動(dòng)分析。

圖9 單加強(qiáng)筋偏心加筋板的前4 階固有振型

圖10 雙加強(qiáng)筋偏心加筋層合板

圖11 四邊固支約束

剛度系數(shù)比分別為

將平板網(wǎng)格劃分成11×11×4，加強(qiáng)筋網(wǎng)格劃分成1×14×4，網(wǎng)格劃分情況如圖12 所示。

圖12 雙加強(qiáng)筋偏心加筋板網(wǎng)格劃分圖

由表 3 可知，對(duì)于材料屬性組成比較復(fù)雜的雙加強(qiáng)筋偏心加筋層合板，非協(xié)調(diào)廣義混合元的固有頻率值明顯低于非協(xié)調(diào)位移元模型的結(jié)果值，并且相應(yīng)的前4 階固有振型如圖13 所示。

表3 雙加強(qiáng)筋偏心加筋層合板前4 階固有頻率(單位：Hz)

4 結(jié)論

本文在非協(xié)調(diào)廣義混合元理論的基礎(chǔ)上，建立了關(guān)于應(yīng)力和位移兩類變量的自由振動(dòng)特征方程。

非協(xié)調(diào)廣義混合元的自由振動(dòng)特征方程的明顯特點(diǎn)是在消去應(yīng)力變量時(shí)，柔度矩陣可對(duì)剛度矩陣中的系數(shù)進(jìn)行調(diào)節(jié)，起到增加結(jié)構(gòu)柔性的作用。

實(shí)例分析表明，基于非協(xié)調(diào)廣義混合元模型的結(jié)果的精度高于非協(xié)調(diào)位移有限元的精度。從理論上講，造成差別的主要原因是非協(xié)調(diào)位移有限元模型的剛度偏硬，會(huì)導(dǎo)致固有頻率結(jié)果偏高，而非協(xié)調(diào)廣義混合元模型的剛度更加接近真實(shí)情況，故而固有頻率結(jié)果更接近真實(shí)解。

圖13 雙加強(qiáng)筋偏心加筋層合板前4 階固有振型