隨歲寒，謝文龍，劉金建

(1.商丘工學(xué)院 機械工程學(xué)院，河南 商丘 476000；2.無錫金元啟信息技術(shù)科技有限公司，江蘇 無錫 214000)

功能梯度材料(functionally graded material, FGM)是一種由兩種或兩種以上結(jié)構(gòu)和性能逐漸變化的材料組成的非均質(zhì)復(fù)合材料，它具有消除應(yīng)力集中、降低殘余應(yīng)力、提高連接強度等普通復(fù)合材料所不具備的優(yōu)良特性。隨著材料科學(xué)的快速發(fā)展，功能梯度材料板在航空航天、機械工程等領(lǐng)域已經(jīng)得到廣泛應(yīng)用，學(xué)術(shù)界關(guān)于該材料板的振動問題也有了很多研究成果。Zhong等[1]研究了功能梯度壓電矩形板的自由振動和強迫振動。Li等[2]基于修正耦合應(yīng)力理論研究了功能梯度壓電微板的靜態(tài)彎曲和自由振動。Li等[3]、王平遠等[4]研究了功能梯度納米板的彎曲、屈曲和振動行為，發(fā)現(xiàn)除材料梯度參數(shù)外，材料的內(nèi)秉特性參數(shù)和非局部參數(shù)對納米板的振動也有重要影響。Jha等[5]提出了一種高階剪切法向變形理論，用于分析功能梯度彈性、矩形簡支板的自由振動。龐有卿等[6]基于高階函數(shù)剪切變形板理論，研究了可移動簡支邊界條件下功能梯度石墨烯增強復(fù)合材料板的自由振動。蹇越傲等[7]基于經(jīng)典板理論,研究了熱載荷作用下功能梯度圓板的非線性動力問題。黃小林等[8]考慮黏彈性地基的影響，基于復(fù)合材料薄板理論研究了石墨烯增強功能梯度矩形薄板的自由振動和動力響應(yīng)。包海軍等[9]研究了功能梯度旋轉(zhuǎn)圓板的熱彈性自由振動。周平等[10]、劉旭等[11]研究了功能梯度軸對稱圓板的自由振動，以及旋轉(zhuǎn)功能梯度納米環(huán)板在熱環(huán)境中的振動頻率。

值得注意的是，現(xiàn)有文獻研究功能梯度板的彎曲及其振動問題，通常考慮板的形狀為矩形、圓形或者圓環(huán)形，以便于建立直角坐標(biāo)或者極坐標(biāo)下的系統(tǒng)控制方程，借助解析法[1,5-6,10]、伽遼金法[8-9]和微分求積法[11]等進行數(shù)值求解。此外，也有學(xué)者利用有限元法研究功能梯度板的彎曲及振動問題，如Deepak等[12]利用ANSYS軟件建模研究了功能梯度方形板的彎曲及自由振動，其建模方法為將板沿厚度方向劃分為若干層，每層作為一個薄板取同一參數(shù)，各層材料參數(shù)不同，采用SHELL板單元離散。這一方式的不足在于，各層內(nèi)部材料參數(shù)恒定且層間材料參數(shù)不連續(xù)，因而求解精度不高且效率低，不適合工程應(yīng)用。

工程中加工及安裝誤差可能造成矩形板的相鄰邊不垂直，使得原設(shè)計的矩形板成為斜板，在這種情況下，上述數(shù)值求解方法將難以完成。等參單元具有較高次的位移模式，能更好地反映結(jié)構(gòu)的復(fù)雜應(yīng)力分布狀態(tài)，同時也能很好地適應(yīng)曲線邊界和非正交的直線邊界。因此，本研究采用8節(jié)點四邊形等參單元離散求解域，利用虛功原理推導(dǎo)功能梯度斜板的質(zhì)量矩陣和剛度矩陣。將本方法所得結(jié)果與ANSYS軟件計算結(jié)果進行了對比，隨后重點分析了梯度指數(shù)和斜角板固有振動的影響，驗證了本方法精度更高。

1 數(shù)學(xué)模型

考慮邊長分別為a和b、厚度為h、斜角為α的功能梯度斜板，建立如圖1所示坐標(biāo)系。假設(shè)板的兩表面分別為陶瓷和金屬，則其彈性模量E和密度ρ可分別表示為

圖1 四邊簡支功能梯度斜板

(1)

采用8節(jié)點四邊形等參離散求解域，則單元廣義坐標(biāo)如下：

jdffdj7=[w1φx1φy1…w8φx8φy8]T。

(2)

將局部坐標(biāo)系下的斜板轉(zhuǎn)化為自然坐標(biāo)系下的方板(圖2)：

圖2 由局部坐標(biāo)系到自然坐標(biāo)系的幾何形狀轉(zhuǎn)換

x=N1x1+N2x2+N3x3+N4x4+N5x5+N6x6+N7x7+N8x8，

y=N1y1+N2y2+N3y3+N4y4+N5y5+N6y6+N7y7+N8y8，

(3)

式中：Ni為形函數(shù)，是自然坐標(biāo)的函數(shù)。根據(jù)等參元概念，構(gòu)造單元中面上任一點的廣義位移模式：

(4)

式中：N=[N10 0N20 0N30 0N40 0N50 0N60 0N70 0N80 0]；

P=[0N10 0N20 0N30 0N40 0N50 0N60 0N70 0N80]；

Q=[0 0N10 0N20 0N30 0N40 0N50 0N60 0N70 0N8]。

按照偏微分規(guī)則，形函數(shù)關(guān)于自然坐標(biāo)的偏導(dǎo)數(shù)可表示為

(5)

(6)

將式(5)和式(6)整理成矩陣形式

(7)

(8)

根據(jù)Reissner-Mindlin板理論，位移場為

u(x，y，z)=-zφx，

(9)

v(x，y，z)=-zφy，

(10)

w(x，y，z)=-w。

(11)

根據(jù)位移場(9)至(11)，得到慣性力做功的變分為

(12)

ε=Bvppzvvn，

(13)

(14)

式中：κ為切變模量修正系數(shù)。將方程(14)簡化為

σ=Dε。

(15)

結(jié)合式(13)和式(15)，得到應(yīng)變能變分為

(16)

將式(12)和式(16)代入如下虛功原理表達式：

δU=δW，

(17)

進而得到系統(tǒng)有限元平衡方程

(18)

需要指出的是，質(zhì)量矩陣M和剛度矩陣K需要結(jié)合式(8)利用高斯積分進行計算，并且這兩個矩陣與斜角α和材料梯度指數(shù)k相關(guān)。

2 數(shù)值算例

采用四邊簡支的邊界條件，求解式(18)可得功能梯度斜板的各階固有頻率，計算所得物理參數(shù)見表1。

表1 物理參數(shù)

圖3給出了ANSYS軟件所得前4階固有頻率及其振型，計算參數(shù)為梯度指數(shù)k=0和斜角α=30°。由式(1)可知，當(dāng)梯度指數(shù)k=0時代表板為均質(zhì)陶瓷板；當(dāng)k=+∞時代表板為均質(zhì)金屬板，此時各階振型與均質(zhì)陶瓷板是類似的。為驗證方法的有效性，利用本方法分別計算了k=0和k=+∞兩種條件下的前4階固有頻率，并與ANSYS軟件計算結(jié)果進行了對比。表2和表3分別為α=30°和α=15°兩種情況下的對比數(shù)據(jù)。從表2可看出，對于均質(zhì)金屬板或者陶瓷板，各階固有頻率誤差百分比相同，隨著階次的升高，誤差略微增大，最大值為第4階時的2.53%。對比表2和表3，當(dāng)斜角α變小時，各階固有頻率的誤差都有較大幅度降低，可見本方法對于小斜角板具備更高精度。

圖3 ANSYS計算前4階固有頻率及其振型(k=0，α=30°)

表2 前4階固有頻率對比(α=30°)

表3 前4階固有頻率對比(α=15°)

圖4和圖5給出了α=30°和α=15°兩種條件下前4階固有頻率隨梯度指數(shù)的變化規(guī)律。隨著梯度指數(shù)的增大，材料中的陶瓷含量逐漸降低，由于陶瓷的彈性模量約為金屬的兩倍，而密度則約為金屬的一半，因而隨著梯度指數(shù)的增大，各階固有頻率逐漸降低且降低幅度逐漸放緩。同時可以看到，在任意梯度指數(shù)下，α=30°時相對于α=15°時各階固有頻率更高。這是由于斜角越大，則板的面積越小、相對剛度越大，從而小面積對應(yīng)著更高的固有頻率。

圖4 前4階固有頻率與梯度指數(shù)的關(guān)系(α=30°)

圖5 前4階固有頻率與梯度指數(shù)的關(guān)系(α=15°)

3 結(jié)語

本研究主要分析了功能梯度斜板的自由振動問題，采用8節(jié)點四邊形等參單元離散求解域，利用虛功原理得到了單元質(zhì)量矩陣和剛度矩陣，求解得到了系統(tǒng)前4階固有頻率，討論了功能梯度板的斜角和梯度指數(shù)對固有頻率的影響。結(jié)果表明，隨著梯度指數(shù)的增大，陶瓷含量降低而金屬含量升高，從而各階固有頻率逐漸降低且降低幅度逐漸放緩。斜角增大使板的剛度增大，進而各階固有頻率增大。均質(zhì)材料條件下將本方法所得結(jié)果與ANSYS軟件計算結(jié)果進行對比，顯示出了良好的一致性。誤差對比結(jié)果也表明，本方法對于小斜角板精度更高。