變截面功能梯度旋轉軸的自由振動特性

李 震，關先磊，王青山，秦 斌

（1.中南大學 高性能復雜制造國家重點實驗室，長沙410083；2.中南大學 交通運輸工程學院 軌道交通安全教育重點實驗室，長沙410075）

旋轉軸作為動力傳輸的重要部件，在數控加工和航空航天等工程領域有著廣泛應用。隨著生產力發展，各向同性材料已經不能滿足旋轉軸高速旋轉的工作要求，迫切需要具有特殊力學性能的新型材料來取代各向同性材料。功能梯度材料(FGM)作為一種新型材料具有改善質量分布、減少熱應力、提高強度等一系列優點，引起國內外學者的關注。

葛仁余等[1]基于歐拉梁理論，采用插值矩陣法（IMM）研究了軸向力作用的軸向功能梯度歐拉梁的自由振動特性；汪亞運等[2]采用切比雪夫多項式得到軸向非均勻變截面梁的自由振動特性，并分析了梯度參數、邊界條件等對固有頻率的影響；韓偉等[3]基于歐拉-伯努利梁理論和動態剛度分析方法推導動態剛度矩陣，研究旋轉變截面梁的自由振動特性；王囡囡等[4]基于哈密爾頓原理，建立中心剛體-柔性梁的動力學模型，研究旋轉變截面梁的自由振動特性；關先磊[5]采用微分求積法來求解梁板殼耦合振動特性；張永旺等[6]基于Raleigh 梁理論建立柔性自轉軸的動力學方程，運用Galerkin 方法求解系統的頻率和模態；黃小林等[7]基于復合材料薄板理論推導出軸向運動功能梯度材料薄板的運動方程，采用Galerkin方法計算其自由振動和屈曲特性；陶彥鳴等[8]采用Chebyshev 譜方法研究變截面歐拉梁的動態特性。Dong等[9]通過假設模態法和拉格朗日方程推導旋轉功能梯度錐形空心圓截面梁的動力學方程，求解梁的自然頻率與模態；Shabanlou等[10]運用高階剪切變形梁理論進行功能梯度旋轉圓柱梁的自由振動特性研究；Huang 等[11]采用Spectral-Tchebychev方法推導旋轉功能梯度錐形圓截面梁的控制方程，研究其旋轉頻率和臨界轉速；Nikbakht 等[12]詳細介紹了功能梯度結構的性能，并對其進行系統的概述。

本文基于鐵木辛柯梁理論得到變截面功能梯度旋轉軸的動能和勢能，將能量表達式用微分求積權系數矩陣表示，得到單元質量矩陣和剛度矩陣,利用有限元法對上述所求的單元矩陣進行組裝得到整體質量矩陣和剛度矩陣以及陀螺矩陣。將基于微分求積有限元(DQFEM)得到的變截面功能梯度軸的無量綱頻率與參考文獻對比，驗證本方法的準確性、穩定性和通用性；緊接著在簡支邊界條件下探究變截面錐度系數和梯度系數對變截面功能梯度旋轉軸固有頻率的影響；最后在簡支邊界條件下探究變截面錐度系數和梯度系數對變截面功能梯度旋轉軸臨界轉速的影響。

1 理論推導

1.1 模型描述

圖1所示為變截面功能梯度旋轉軸的動力學模型。該模型基于鐵木辛柯梁理論，考慮4個自由度，用us和vs分別表示沿x軸和y軸方向移動；θx和θy分別表示繞x軸和y軸方向旋轉。其材料屬性沿z軸方向連續變化，即材料密度ρ、彈性模量E、剪切模量G是關于z的變量。

圖1 變截面功能梯度旋轉軸

1.2 動 能

如圖2所示，在靜止坐標系下，以轉軸中心線上的一點為研究對象，其中rO0OP表示轉軸上一點的位移矢量，i0、j0、k0、i1、j1、k1分別表示靜止坐標系和連體坐標系下的方向向量，通過式(1)在靜止坐標系和連體坐標系下分別求解位移矢量rO0O1和rO1P,并進行矢量疊加得到在靜止坐標系下轉軸上任意一點的位置矢量r0op。

圖2 矢量示意圖

對于轉軸上一點轉動的描述，可將旋轉軸類似剛體進行處理，即將轉軸上一點轉動轉化為剛體定點轉動。本文采用卡爾丹角描述剛體定點轉動，結果如圖3所示：

(1)坐標系xyz繞x軸轉動α得到坐標系x1y1z1；

(2)坐標系x1y1z1繞y1軸轉動β得到坐標系x2y2z2；

(3) 坐標系x2y2z2繞z2軸轉動γ得到坐標系x3y3z3(ξηζ)；

通過上述3次轉動實現坐標系xyz到坐標系ξηζ之間的轉換，式（2）是實現xyz到ξηζ的坐標變換矩陣。

圖3 卡爾丹角

采用卡爾丹角描述剛體的定點運動，角速度的合成也是采用矢量疊加，在連體坐標系下任意一點角速度分別用ωb表示。

由于α、β、γ值較小，參考《多剛體動力學基礎》[13]對變換矩陣進行簡化，得到矩陣A1。

通過簡化后坐標變換矩陣，可以求得在連體坐標系下轉軸上任意一點相對于靜止坐標系的角速度ωrb。

對于剛體上一點轉速的計算，參考《理論力學》[14]中剛體運動學中運動坐標與牽連速度的描述，計算轉軸上一點轉速。

將根據式(6)求得轉速代入到式(7)計算轉軸動能。

從而通過計算可以得到靜坐標系下轉軸動能的表達式(8)，其中S表示軸的截面積，I表示截面慣性矩。

1.3 勢能

基于鐵木辛柯梁理論，考慮剪切變形和彎曲變形，轉軸勢能表達式見式(9)，具體推導過程參考文獻[15]，其中E表示彈性模量；I表示截面慣性矩；κ表示剪切因子；G表示剪切模量；S表示截面積；其轉軸勢能表達式可以表示為

1.4 微分求積的原理

微分求積是基于數值分析方法將描述邊值問題的微分方程轉化成一組代數方程的形式，實現快速求解。設自變量x的一維函數為f(x)，在區間[-1，1]上連續可導。在區間[-1，1]上取N個不同節點，則函數在xi的1階導數可以表示為

式中：A(ij1)，i，j=1，2，…，N，是1階權系數，上角標“(1)”表示函數取1階導數。由權系數A(ij1)組成的矩陣A(1)為1階權系數矩陣，f為單元位移列向量。

微分權系數矩陣A(1)需要選取相應的插值基函數來獲得，采用拉格朗日作為插值基函數，因其具有較好的遞推關系。

類似于微分處理的節點劃分，用數值積分方法在[-1，1]的積分區間對函數進行積分，則有：

積分權系數矩陣C中的元素可用常規的數值積分方法求解，f是單元位移列向量。本文采用高斯洛巴托求積方法求解，其中權系數為

Pn(ξj)是n階勒讓德多項式，ξj(j=1，2，…，n-1，)是Pn-1(ξj)1階導的第(j-1)個零點，Cj表示積分權系數矩陣中的項。為了使權系數矩陣適用于不同長度軸單元的求解，需要對權系數矩陣進行修正：

其中le表示轉軸單元的長度。

1.5 微分求積的引入

將能量表達式當中的項用微分和積分權系數矩陣以及單元位移列向量進行表示。

式中：δe表示單元位移列向量；ue、ve、φex和φey分別表示各個自由度上的單元位移列向量。本文討論變截面功能梯度旋轉軸，其彈性模量E、剪切模量G、密度ρ、截面面積S、慣性矩I沿軸向變化，都是N維對角矩陣。例如：E=diag(E（z1），E（z2），…，E(zN))，其中diag表示對角矩陣。

將上述所有參數代入勢能表達式，得到下列微分求積形式的表達式，求解單元剛度矩陣：

式中：K1表示關于轉軸單元的剛度矩陣。

式中：M1和M2是關于轉軸單元的質量矩陣。

由此得到了轉軸動能和勢能矩陣形式的表達式，為了求解轉軸單元的質量矩陣、陀螺矩陣和剛度矩陣，本文將上述動能和勢能的矩陣形式代入到拉格朗日方程求解系統總體的矩陣方程：

根據式（21）可以得到轉軸單元的矩陣方程，進一步確定轉軸單元的質量矩陣Me、剛度矩陣Ke以及陀螺矩陣Ge的組裝矩陣，具體的表達式如下：

根據式（22）可以確定轉軸單元的質量矩陣Me、陀螺矩陣Ge和剛度矩陣Ke。為了得到轉軸系統整體的質量矩陣M，陀螺矩陣G和剛度矩陣K，本文采用“對號疊加法”將轉軸的單元矩陣集成到總體矩陣中。對號疊加法的關鍵是根據節點號和自由度號得到單元矩陣在整體矩陣當中行號和列號，由于每個節點的自由度是4，所以節點號為n的第p個自由度在整體矩陣中的行號（或列號）為Q=（n-1）*4+p。本文采用微分求積的積分節點取值為n1=19，且n=n1+1。緊接著對已獲得的M、G和K采用組裝矩陣方式進行處理，其組裝矩陣表達式如式(23)所示：

緊接著借用MATLAB平臺采用eig函數求解相應的模態頻率和模態向量。相應的表達式如下：

式中：V和F分別表示模態向量和模態頻率。

1.6 變截面理論

變截面軸分為兩種，分別是面積線性變化和半徑線性變化。本文討論半徑沿軸向線性變化的變截面旋轉軸，變截面主要引起截面面積以及截面慣性矩的變化。設左側截面半徑為r0，右側截面半徑為r1，變截面軸的長度為L，變截面錐度系數為c（一般取值為0～0.8），z表示距離左側的軸向距離。

1.7 功能梯度理論

功能梯度材料是一種非各向同性材料，彈性模量和密度沿軸向呈冪等規律變化，其變化規律取決于梯度系數m（一般取0至10）。Es表示左側截面彈性模量，ρe表示左側截面密度，Es表示右側截面彈性模量，ρe表示右側截面密度，其余參數同上。

2 算例分析

為了便于編程求解變截面功能梯度旋轉軸的自由振動特性，本文給出了相應的算法流程圖，如圖4所示。

2.1 算法驗證

為了驗證本方法的收斂性和準確性，采用文獻中的經典案例進行驗證，為了便于與文獻進行對比令轉速v=0，同時假定功能梯度材料是由ZrO2和Al構成，其材料屬性沿著軸向呈冪等規律變化，具體的變化規律如式(24)所示。ZrO2材料的彈性模量和密度分別是Es=200 GPa和ρs=5 700 kg/m3；Al材料的彈性模量和密度分別是Ee=70 GPa 和ρe=2 702 kg/m3；變截面軸的長度L=1 m；軸的大徑r1=0.1 m；軸的細長比r=0.01；剪切因子為5/6；泊松比v=0.3；功能梯度系數m=2。本文均采用上述參數進行研究，為了便于與文獻進行對比，需要將求得的頻率無量綱化，引入固有頻率無量綱化的公式，關于細長比（回轉半徑）定義如下：

圖4 算法流程圖

計算結果如表1所示。

從表1看出。DQFEM算法和有限元算法[16]的計算結果基本吻合。可以驗證該算法的收斂性和正確性。

關于變截面功能梯度旋轉軸轉速的研究目前沒有能找到相關的文獻進行對比研究，本文在簡支邊界條件下，取變截面功能梯度旋轉軸的極端情況(m=0；m=∞)與ANSYS的仿真結果進行對比驗證(ANSYS 中選用BEAM188 單元，劃分100個單元)，變截面功能梯度旋轉軸的轉速為500 rad/s。

從表2中看出,采用DQFEM算法與ANSYS算法在兩種極端情況下得到的固有頻率誤差在1%以內，驗證施加轉速時的正確性。在此基礎上進行變截面功能梯度旋轉軸轉速的研究。

表1 不同邊界條件不同錐度系數的無量綱頻率

表1 不同邊界條件不同錐度系數的無量綱頻率

?

表2 變截面功能梯度旋轉軸極端情況下固有頻率對比

圖5 無量綱固有頻率變化

2.2 變截面錐度系數和梯度系數對固有頻率的影響

上述對比驗證本文方法的正確性，結合工程實際應用，本文僅在簡支邊界條件下探究變截面錐度系數以及梯度系數對功能梯度旋轉軸固有頻率的影響。同時指出此時的轉速為0，其他具體參數如2.1節所示，相應的結果如圖5所示。

如圖5所示，在簡支邊界條件下，變截面功能梯度軸的無量綱固有頻率隨著變截面錐度系數的增加呈現出降低的趨勢,隨著梯度系數增加，呈現出先逐漸增加，后保持不變的變化趨勢。從上述4階無量綱固有頻率的變化曲線中看出，當梯度系數大于2以后，梯度系數的增加對無量綱固有頻率幾乎不產生影響。

2.3 變截面錐度系數和梯度系數對臨界轉速的影響

由上所述，臨界轉速是工程應用中極為重要的參數，因此有必要研究變截面錐度系數和梯度系數對臨界轉速的影響。本文首先研究變截面功能梯度旋轉軸在簡支邊界條件下以及特定變截面錐度系數和梯度系數條件下轉軸前后渦動頻率隨轉速的變化規律，其中梯度系數m=2，變截面錐度系數c=0.8，其他具體參數如2.1節所示，得到如圖6所示的整體坎貝爾圖。

圖6 整體坎貝爾圖

從圖6中不難看出，變截面功能梯度旋轉軸渦動頻率（固有頻率）隨轉速增加出現分岔現象，分為前渦動和后渦動，前渦動頻率隨轉速的增大而增大，后渦動頻率隨轉速的增大而減少。圖6中斜線是描述渦動頻率與轉速相等的一條直線，該直線的表達式為ω=Ω。同時指出，該直線與前渦動頻率的交點所對應橫坐標值是臨界轉速。圖中ΩI、ΩII、ΩIII和ΩIV分別表示1階、2階、3階和4階臨界轉速。為了能夠直觀描述上述現象，本文將圖6中的坎貝爾圖進行局部放大，結果如圖7所示。

同時有必要給出變截面功能梯度旋轉軸臨界轉速的表達式，以便開展關于臨界轉速的研究。首先求解前渦動頻率直線的解析表達式，然后將表達式與解析式為ω=Ω的直線聯立求解得到臨界轉速。假設轉速為Ω=0時，前渦動頻率為ω1；轉速為Ω=1 000時，前渦動頻率為ω2，則渦動曲線的表達式為

圖7 局部放大坎貝爾圖

將式（28）與ω=Ω進行聯立可得臨界轉速的表達為

基于以上討論，確定了臨界轉速的計算方法，緊接著研究變截面功能梯度旋轉軸在簡支邊界條件下變截面錐度系數和梯度系數對臨界轉速的影響，具體參數見2.1節，相應的結果如圖8所示。

圖8 簡支邊界條件下臨界轉速的變化

從圖8可以看出，在簡支邊界條件下，研究變截面功能梯度旋轉軸臨界轉速隨變截面錐度系數以及梯度系數的變化，4階臨界轉速的變化趨勢相同，隨著變截面錐度系數的增加，臨界轉速逐漸降低，隨著梯度系數增加，臨界轉速先增加后保持不變。

3 結語

(1)本文計算結果分別與已有文獻結果和有限元ANSYS 計算結果的對比驗證證明本文方法具有準確性、穩定性和通用性，適用于求解變截面功能梯度旋轉軸的自由振動問題。

(2)在模型驗證的基礎上，開展了關于變截面功能梯度旋轉軸固有頻率的研究，主要探究變截面錐度系數和梯度系數對固有頻率的影響，結果表明上述參數對固有頻率有重要影響。

(3)根據已建立的模型，進一步研究了變截面功能梯度旋轉軸渦動頻率隨轉速的變化規律，給出了臨界轉速的定義，并進一步探究變截面錐度系數和梯度系數對臨界轉速的影響，對實際工程應用具有重要的意義。