基于能量有限元法的功能梯度梁振動分析

王 迪, 朱 翔，2, 李天勻，2, 衡 星, 高 雙

(1.華中科技大學 船舶與海洋工程學院，武漢 430074；2.船舶與海洋水動力湖北省重點實驗室，武漢 430074)

功能梯度材料(Functionally Graded Material，FGM)[1]通常是由兩種或兩種以上不同性能的材料組成，材料的基本要素沿結構的某一方向連續變化。功能梯度材料最早是由日本學者提出，由于其優良的結構性能和重要的應用價值，引起人們廣泛的研究和關注。

梁結構是工程中最為廣泛使用的基本結構，自功能梯度材料提出以來，許多學者對功能梯度梁結構的振動特性進行了研究。鄧昊等[2]基于動力學的相關理論，通過對狀態空間變量的替換，研究了功能梯度Timoshenko梁的理論模型。Yang等[3]利用經典的歐拉梁理論，對不同邊界條件的功能梯度梁進行了自由振動分析，研究了功能梯度梁的振動特性。Simsek[4]將梁的經典理論、一階剪切和高階剪切理論與拉格朗日方程結合起來，分析了在不同邊界條件下功能梯度材料梁的固有頻率和振動特性。Alshorbagy等[5]應用虛功原理，推導了材料性質沿厚度方向冪函數變化的功能梯度材料梁自由振動的平衡方程并采用有限元數值方法進行求解對比。Ravikiran等[6]基于高階剪切理論，采用有限元法得到陶瓷-金屬功能梯度材料梁的平衡方程，應用剛度矩陣法求得了功能梯度材料梁的撓度解。

在研究結構振動的分析方法中，能量有限元法(Energy Finite Element Method, EFEM)是近些年發展起來的一種新的數值方法，其主要目的是解決結構中高頻的振動問題[7]。該方法將能量密度作為動力學方程的基本變量，通過有限元的思路將不同結構離散成單元來表示結構的特性，通過求解能量密度進而得到結構的振動特性。孫麗萍等[8]以能量密度為變量,推導了薄板結構受激勵產生彎曲振動時的能量密度控制方程,并得到了控制方程的有限元解。蔡忠云等[9]以能量有限元法建立了復合材料層合梁受激勵時的控制方程, 并研究了復合材料梁的橫向振動問題。解妙霞等[10]利用能量有限元方法研究圓柱殼體結構，并根據薄殼理論和能量平衡關系推導了圓柱殼彎曲振動時能量密度的控制方程。而目前關于FGM結構的能量有限元研究尚未見到報道。

本文基于能量有限元法對功能梯度梁以及耦合功能梯度梁的振動特性進行研究。通過FGM梁的基本動力學方程推導得到FGM梁的能量密度控制方程和有限元矩陣方程。在FGM梁的算例分析中，將本文的能量有限元法和傳統有限元法得到的能量密度解進行對比和分析，驗證能量有限元法的準確性。在此基礎上對耦合功能梯度梁結構的能量密度和能量流進行求解。

1 FGM梁的能量有限元方程

1.1 一維FGM梁結構的能量密度控制方程

FGM梁在定義的時候，通常假定材料的各項屬性沿著梁截面厚度方向以指數形式變化。坐標系如圖1所示，x為梁的長度方向，z為梁的高度方向。材料的屬性變化主要是體現在楊氏模量E(z)和材料密度ρ(z)這兩個方面

E(z)=E0eβz

(1)

ρ(z)=ρ0eβz

(2)

式中：E0和ρ0為z=0處的彈性模量和密度值；β=ln(E2/E1)/h為FGM梯度變化的常數。

圖1 功能梯度梁模型坐標系

考慮一個受橫向簡諧力激勵下的含阻尼的功能梯度梁，如圖2所示，采用簡單梁理論，其運動控制微分方程可以表示為[11]

(3)

圖2 受到簡諧激勵作用下的功能梯度梁模型

假設方程的通解為

w(x,t)=(A1e-ikfx+A2eikfx+A3e-kfx+A4ekfx)eiωt

(4)

式中：kf為復波數；A1、A2、A3和A4為方程待求解的系數，由邊界條件來確定。

方程的通解式(4)代入運動方程式(3)得到FGM梁的彎曲波波速的表達式為

(5)

由結構內部的能量平衡關系可知，橫向振動梁的能量密度和功率遠場解局部均值有如下的關系式[12]

(6)

通過加權余量法將上述偏微分方程轉化為代數方程組，將權函數作為形函數，則可通過Galerkin加權殘值法將梁的能量密度控制微分方程表示為[13]

(7)

同時，梁的能量流遠場解與能量密度的遠場解存在如下關系[14]

(8)

可以將式(7)用矩陣的形式表示出來

[Ke]{ee}={Pe}+{Qe}

(9)

式中：[Ke]為每個梁單元有關質量和剛度信息的系數矩陣；{ee}則是需要求解的各個節點的能量密度矩陣；{Pe}為在節點處輸入的激勵能量；{Qe}則為每個梁單元兩端進出的能量流。

因此，根據能量有限元法求解式(9)可得到節點的能量密度，即可進一步得到每個梁單元的能量流。

1.2 耦合梁結構的能量有限元方程

對于更加復雜的耦合結構，由于耦合結構會在耦合處存在材料或幾何的不連續，從而導致結構振動波在耦合處產生反射和透射。通常能量密度在耦合節點處是不連續的，而能量流是連續的。因此在建立能量有限元系數矩陣的同時，需要單獨考慮耦合節點處的能量傳遞和連接，可以通過連接矩陣[15]的形式來解決耦合結構的連接問題。

圖3 功能梯度梁模型中的能量傳遞過程

當僅有彎曲波由第一段梁向第二段梁入射時，兩段梁的橫向位移可以分別表示為

(10)

(11)

梁中彎矩M和剪力N可以表示為

(12)

(12)

式中：I為梁的截面慣性矩。

根據耦合節點處彎矩、剪力、位移以及轉角的連續，可以得到四個連續條件

M1=M2

(14)

N1=N2

(15)

W1=W2

(16)

(17)

通過求解上述四個方程可得

(18)

(19)

(20)

(21)

(22)

忽略近場效應，可以得出梁耦合節點處的能量透射系數τ12和反射系數γ11分別為

(23)

(24)

在耦合節點處，能量流有正負兩個方向，+代表正方向、-代表負方向，如圖4所示，并且有如下的關系式[16]。

圖4 耦合節點處的能量流

(25)

(26)

波在耦合節點處會同時產生反射和透射，主要有如下的關系式

(27)

(28)

式中：γ和τ分別為耦合梁處的反射系數和透射系數，其滿足關系式τ+γ=1。將式(27)和式(28)聯立，再利用耦合節點處功率連續性q1m+q2m=0的條件可以得到

(29)

同時，在耦合節點處，能量密度是由正負兩個方向的能量密度值之和組成，有如下關系式

(30)

(31)

將式(25)～式(28)代入式(30)、式(31)可得

(32)

聯立式(29)和式(32)可以得到

(33)

將關于能量流的項進行統一的整理，可以將系統的整體矩陣表示為

[Ke+J]{ee}={Pe}+{Qe}

(34)

由式(34)即可對耦合FGM梁的能量密度求解，從而得到振動能量。

2 FGM梁結構的能量有限元振動分析

2.1 FGM梁的能量密度計算驗證

本節按照上節推導的公式，對某FGM梁的能量密度進行計算，得到結構中能量密度值的分布。同時也通過傳統有限元模型對該梁的能量密度進行計算并對比，從而驗證EFEA方法的正確性。

考慮圖1所示的簡支梁模型，其長度為4 m，截面是邊長為0.02 m的正方形，梁的上表面材料為陶瓷，密度ρ為3 800 kg/m3，彈性模量E為3.8×1011Pa，μ為0.3，梯度變化常數β取85，因此可得到下表面的材料參數。結構的阻尼損耗因子η為0.02。在梁的中心處作用一輸入功率為1.12×10-3W的橫向激勵力，激勵頻率為3 000 Hz。

按照上一節的推導，對整個梁劃分能量有限元單元，單元數為10，形函數取n=3的Lagrange插值函數，計算得到能量密度分布如圖5所示，其中參考能量密度值為1×10-12J/m。

為了驗證本文能量有限元計算的準確性，同時也采用有限元法計算梁中的能量密度。有限元模型中的FGM梁采用8節點實體單元來模擬，梁的高度方向分成10份，每層梁單元按照梯度公式賦予不同的E和ρ來模擬功能梯度材料的變化，長度方向劃分100份，通過提取中性面節點的速度來計算能量密度，如圖5。從圖中可見由于能量有限元法計算得到的是時間和空間平均后的解，因此能量密度分布比較均勻，而有限元解反映出局部特征，其得到的能量密度隨著EFEM的解上下波動。但是兩種方法的局部誤差均不超過1 dB。

圖5 兩種方法沿梁長度方向能量密度分布的結果對比

由于能量有限元模型和傳統有限元模型的單元數量不同，且能量有限元方法求得的能量密度是對時間和空間平均后得到的相對平滑的解，因此也可將有限元法得到的能量密度值按照能量有限元法設定的單元長度選出對應包含的節點，將這些節點上的能量密度進行一次平均，然后來比較兩種方法的結果，如表1所示。從均值的量化結果來分析，兩種方法也吻合較好，從而說明本文針對FGM梁的能量有限元分析是準確的。

表1 EFEA和FEA(均值)結果量化對比分析

與傳統的有限元方法相比，能量有限元法的單元數量少，計算速度快，因此在計算效率上有了很大程度的提升，但是和有限元相比，其求解的能量反映的是時間和空間上的均值。與統計能量法相比，由于統計能量分析的基本變量是在一個子系統內對各組相似模態的能量密度在指定頻寬內頻率平均，而能量有限元則是同時對能量密度在一個波長的空間平均和在某個頻寬內的頻率平均，因此統計能量法只能求解子系統的能量，而能量有限元法能夠求解出每個單元的能量，且對結構模態數沒有要求，分析頻段可擴展到中高頻段。

2.2 耦合FGM梁的振動分析

如圖6所示，考慮兩段同軸的耦合FGM梁，每段梁長度為2 m，兩段梁截面面積不同。第一段梁的截面是邊長為0.01 m的正方形，第二段梁的截面是邊長為0.02 m的正方形。兩段梁的屬性和3.1節相同，在第二段梁上施加輸入功率為2.24×10-3W的橫向激勵力，激勵頻率為3 000 Hz。按照能量有限元法，將每段FGM梁劃分成十個單元，然后計算梁單元節點處的能量密度值和梁單元的能量流。將計算得到的能量密度和能量流沿梁的分布表示出來，如圖7和圖8所示。

圖6 耦合梁的EFEA分析模型

圖7 耦合梁的能量密度分布

從圖7和圖8可見，耦合梁在耦合位置處能量密度發生了突變，這是由于耦合梁的截面產生突變，造成能量密度在耦合節點處不連續，同時可見在能量輸入的位置附近能量密度較為集中。在圖8中由于能量平衡，因此能量流在耦合位置處是連續的，沒有發生突變，并且在激勵輸入處能量流達到峰值。同時，由于梁中阻尼的存在，遠離激勵位置的能量流逐漸衰減。

2.3 截面尺寸變化對耦合FGM梁的振動影響

考慮到耦合梁的截面參數變化會改變耦合面處的透射系數和反射系數，對連接矩陣產生影響，因此本節討論截面尺寸變化對耦合FGM梁的振動影響。本節算例中除了梁的截面邊長發生變化，其他參數均和上節相同。保持第一段梁的截面邊長不變，改變第二段梁的截面邊長，對比不同邊長比下的能量密度和能量流結果，如圖9和圖10所示。

圖8 耦合梁的能量流分布

圖9 不同邊長比下耦合梁的能量密度分布

圖10 不同邊長比下耦合梁的能量流分布

從圖9和圖10可見，兩段梁的邊長比變化對耦合梁的能量密度和能量流都會產生一定的影響。隨著邊長比的增大，兩段梁的能量密度和能量流的差別也逐漸增大。當邊長比達到3時，能量流在激勵點的位置達到最大，此時連接處截面突變最明顯，因此耦合處能量密度的突變也最明顯。這是由于邊長比變化會通過能量透射系數和反射系數引起連接矩陣的改變，隨著邊長比增大，能量在兩段梁的分布差別也會逐漸增加。

為了進一步研究不同邊長比下梁的振動規律，本節還計算了耦合梁單元的速度分布，如圖11所示。從圖中可以看出，由于耦合處截面的不連續，不同單元速度在耦合處依然是不連續的。相比第一段梁，由于第二段梁的截面發生變化，彎曲剛度也隨之變化，因此第二段梁單元的速度變化更加明顯。隨著邊長比的不斷增大，單元速度值也隨之減小。

圖11 不同邊長比下耦合梁的單元速度分布

3 結 論

本文基于能量有限元法對功能梯度梁和耦合功能梯度梁的振動特性進行了研究，推導了梁和耦合梁的能量有限元法基本方程，并推導了耦合FGM梁的連接矩陣和能量有限元矩陣方程。通過將EFEA計算得到的單元能量密度與傳統有限元FEM的結果進行對比，驗證了本文推導和求解的正確性。通過對耦合梁的能量密度和能量流計算可見，耦合FGM梁在耦合位置處能量密度發生了突變。能量流在耦合節點處保持連續，并且在能量的輸入位置能量流達到峰值。由于阻尼的存在，遠離激勵位置的能量流逐漸衰減。耦合梁的截面變化越明顯，耦合梁之間的能量密度變化也越明顯。本文的研究為能量有限元法應用到復雜功能梯度材料結構的振動研究提供了借鑒。

[1] 李信，劉海昌，滕元成，等. 功能梯度材料的研究現狀及展望[J]. 材料導報,2012(增刊1):370-373.

LI Xin, LIU Haichang, TENG Yuancheng，et al. Research status and future directions on functionally gradient materials[J]. Materials Review,2012(Sup1):370-373.

[2] 鄧昊，程偉. 一種功能梯度Timoshenko梁的損傷識別方法[J]. 北京航空航天大學學報,2016,42(10):2214-2221.

DENG Hao, CHENG Wei. A damage identification method for functionally graded Timoshenko beams[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016,42(10):2214-2221.

[3] YANG J, CHEN Y. Free vibration and buckling analyses of functionally graded beams with edge cracks[J]. Bulletin of the Polish Academy of Sciences Mathematics, 2008, 83(1):48-60.

[4] SIMSEK M. Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories[J]. Nuclear Engineering and Design, 2010,240(4): 697-705.

[5] ALSHORBAGY A E, ELTAHER M A, MAHMOUD F F. Free vibration characteristics of a functionally graded beam by finite element method[J]. Applied Mathematical Modeling, 2011,35(1):412-425.

[6] RAVIKIRAN K, KASHIF A, GANESAN N. Static analysis of functionally graded beams using higher order shear deformation theory[J]. Applied Mathematical Modeling, 2008,32(12): 2509-2525.

[7] 原凱，王建民，韓麗，等. 能量有限元在振動與噪聲預示中的研究進展[J]. 強度與環境,2015,42(3):10-19.

YUAN Kai, WANG Jianmin, HAN Li, et al. Review on prediction of vibration and noise using energy finite element analysis[J]. Structure & Environment Engineering,2015,42(3):10-19.

[8] 孫麗萍，聶武. 能量有限元法在船舶結構中的應用[J]. 哈爾濱工業大學學報,2008,40(9):1491-1494.

SUN Liping，NIE Wu. Application of energy finite element method in ship structures[J]. Journal of Harbin Institute of Technology,2008,40(9):1491-1494.

[9] 蔡忠云，唐文勇，張圣坤. 能量有限元方法在復合材料層合梁耦合結構振動分析中的應用[J]. 振動與沖擊,2010,29(10):23-27.

CAI Zhongyun, TANG Wenyong, ZHANG Shengkun. Application of energy finite element method in vibration analysis of coupled composite laminated beam structures[J]. Journal of Vibration and Shock,2010,29(10):23-27.

[10] 解妙霞，陳花玲，吳九匯. 圓柱殼高頻彎曲振動的能量有限元分析[J]. 西安交通大學學報,2008,42(9):1113-1116.

XIE Miaoxia, CHEN Hualing, WU Jiuhui. Energy finite elenment analysis to high-frequency bending vibration in cylindrical shell[J]. Journal of Xian Jiaotong University,2008,42(9):1113-1116.

[11] AYDIN K. Free vibration of functionally graded beams with arbitrary number of surface cracks[J]. European Journal of Mechanics - A/Solids, 2013, 42(42):112-124.

[12] 朱繼清. 基于能量流有限元的耦合結構振動傳遞特性分析[D].濟南：山東大學,2011.

[13] 楊旸. 能量有限元分析法的研究及其在復合材料層合板中的應用[D]. 寧波：寧波大學,2013.

[14] 孫麗萍. 能量有限元法研究及其應用[D]. 哈爾濱：哈爾濱工程大學,2004.

[15] 游進,李鴻光,孟光. 耦合板結構隨機能量有限元分析[J]. 振動與沖擊,2009,28(11):43-46.

YOU Jin，LI Hongguang，MENG Guang. Random energy finite element analysis of coupled plate structures[J]. Journal of Vibration and Shock, 2009,28(11):43-46.

[16] 徐福慧,鄭世杰.基于能量有限元法的耦合梁分析[J]. 江蘇航空,2012(增刊1):108-110.

XU Fuhui，ZHENG Shijie. Analysis of coupled beam based on energy finite element method[J]. Jiangsu Aviation, 2012(Sup1):108-110.

[17] 鄭國濤. 能量有限元方法若干問題的研究[D]. 上海：上海交通大學,2013.