旋轉復合材料薄壁軸的自由振動與穩定性＊

任勇生，代其義，張興琦

（山東科技大學機械電子工程學院，山東 青島266590）

引 言

在包括能源、航空和汽車在內的許多工業技術領域內，廣泛地存在著一類旋轉柔性細長結構，它的轉動角速度矢量與其結構的縱軸保持平行，在能量傳動過程中發揮著無可替代的重要作用，力學上通常將其稱之為旋轉軸。例如燃氣輪發電機傳動軸、航空發動機傳動軸以及汽車后輪驅動軸等等。

纖維復合材料由于比強度和比剛度高、抗疲勞和減振性能好，在直升機尾傳動軸［1］以及汽車傳動軸［2］的結構設計中已經顯示出廣闊的應用前景。采用輕質纖維復合材料取代傳統的金屬材料不僅可以減輕結構的重量，同時還能夠減少噪聲、提高結構的抗振性能［3］。然而由于纖維復合材料力學性能的各向異性，加之軸的設計一般采取薄壁結構形式，因此在軸向拉伸、橫向彎曲以及扭轉變形之間，存在著顯著的彈性耦合。精確地分析復合材料軸的振動特性是對其穩定性評估和參數優化設計的基礎。傳統的基于各向同性假定的金屬材料軸的動力學模型與分析方法顯然已經不再適合于各向異性復合材料傳動軸，因此，需要建立更為先進的模型來指導復合材料軸的動力學設計。

迄今為止，在實心梁（Euler－Bernolli梁和Timshenko梁）理論、圓柱殼理論和薄壁梁理論的框架下，人們已經提出了復合材料軸的一些動力學分析模型。文獻［4］將Timoshenko梁理論和Donnell薄殼理論結合在一起，建立旋轉復合材料軸的剛度矩陣并采用有限元方法導出系統的運動方程，由此計算了薄壁復合材料軸的臨界轉速。文獻［5］基于等效模量梁理論（Equivalent modulus beam theory，EMBT）對復合材料軸進行了動力學建模，并且將臨界轉速的理論結果與實驗結果進行了對比。文獻［6］分別基于EMBT和分層梁理論（Layerwise beam theory，LBT）建立復合材料軸的動力學模型，研究發現對于非對稱鋪層方式，采用兩種模型得到的臨界轉速存在偏差。文獻［7］依據殼的一階近似理論推出了薄壁復合材料傳動軸的運動微分方程。并且應用該模型分析計算了不同類型的復合材料傳動軸的臨界轉速。文獻［8］采用Timoshenko梁理論建立了復合材料傳動軸的動力學方程。該模型考慮了陀螺效應和彎扭耦合的影響。文獻［9］根據Timoshenko梁理論和Hamilton原理推導了旋轉復合材料錐形軸的動力學方程。研究發現錐度比對軸的固有頻率有顯著影響。文獻［10］基于一階剪切梁理論提出了一個復合材料軸系統的有限元動力學模型，該模型除了軸還包含了剛盤以及軸承。文獻［11］和［12］基于Rehfied的復合材料薄壁梁理論［13］，建立了復合材料軸的振動微分方程，研究了矩形截面軸和圓形截面軸的固有頻率和穩定性特性，其中包括橫向剪切的影響。

變分漸進法（Variational asymptotic method，VAM）［14］是從二維殼能量函數出發，采用漸進分析建立的一種精細的復合材料薄壁梁理論，與其他薄壁梁理論相比，它的主要特點在于其位移場表達式除了包含經典的扭轉翹曲之外，同時還包含由于軸向拉伸以及橫向彎曲變形引起翹曲。VAM復合材料薄壁梁理論不僅適用于單閉室復合材料薄壁梁，而且也適用于含有多個閉室的復合材料薄壁梁。自上世紀90年代VAM復合材料薄壁梁理論問世以來，國外一些學者已經將其成功地應用于先進直升機復合材料葉片的結構建模和動力學分析［15－16］。然而，迄今為止，作者還沒有見到將VAM復合材料薄壁梁理論用于復合材料軸動力學分析的研究報道。此外，大家也注意到，現有的VAM復合薄壁梁理論尚未考慮橫向剪切變形的影響。

本文從VAM復合材料薄壁梁理論出發［14］，提出了一個旋轉復合材料軸的動力學分析模型，其中引入了橫向剪切變形的影響。基于Hamilton原理推導出復合材料傳動軸的振動偏微分方程組。采用Galerkin法求解得到耦合振動固有頻率和臨界轉速的近似解，并且給出了算例驗證。采用本文建立的模型計算方法，計算和分析了周向均勻剛度配置（CUS）［12］復合材料傳動軸的動力學特性。通過數值計算揭示了鋪層角、旋轉速度、長徑比、徑厚比以及橫向剪切對復合材料傳動軸自由振動固有頻率和穩定性的影響。

1 振動微分方程的建立

圖1表示長度為L的、繞其軸線以定常角速度Ω旋轉的封閉截面纖維復合材料薄壁梁。旋轉坐標系為（x，y，z），局部坐標系為（x，s，ξ），其中環向坐標s沿著薄壁梁中面切線逆時針方向，ξ沿著薄壁梁中面法線方向。

圖1 圓截面復合材料薄壁軸Fig.1 Composite thin－walled shaft of a circular cross section

從文獻［14］的位移場表達式出發，進一步引入橫向剪切變形的影響，于是薄壁軸橫截面上的任意一點沿著x，y，z方向的位移假設如下

式中U1（x，t），U2（x，t），U3（x，t）分別表示橫截面沿著x，y，z方向的剛體位移；φ（x，t），θy（x，t），θz（x，t）分別表示橫截面繞x軸的扭轉角以及繞y，z軸的扭轉角。y，z表示橫截面中心圍線上的點的坐標，是環向坐標s的函數。

假定薄壁梁的翹曲函數g（s，x，t）具有如下形式上述等式右端的四項依次是與扭轉、軸向拉伸、繞z軸彎曲和y軸彎曲有關的翹曲分量，其中G（s）的物理意義為扭轉率，g1（s）的物理意義為軸向應變，g2（s）和g3（s）的物理意義為沿y，z方向的彎曲曲率。

在方程（1）和（2）中，θy（x，t），θz（x，t）可以表示如下

根據文獻［14］幾何方程，由位移方程（1）可以導出橫截面正應變γxx和面內剪應變γxs的表達式，并且依照文獻［17］，對橫向剪應變γxξ的表達式也作出假設。因此，考慮剪切變形的薄壁軸的幾何方程可以寫成如下

定義薄壁軸內力

薄壁軸內力－應變本構關系

其中

為了導出復合材料軸的振動方程，利用Hamilton原理

式中U和T分別為應變能和動能，分別由下式確定

式中σxx，σxs，σxξ分別表示橫截面正應力、面內剪應力和橫向剪應力；εxx＝γxx，εxs＝2γxs，εxξ＝2γxξ是相應的工程應變。

式中ρ表示材料密度，V表示變形后的梁上任意一點的速度矢量，它與變形梁任意一點的位置矢量：r＝（y＋u2）i＋（z＋u3）j＋（x＋u1）k之間滿足關系：。

定義復合材料軸的軸力Fx，扭矩Mx，彎矩My和Mz，剪力Qy和Qz如下

將式（6）代入式（11），并且利用幾何方程（4），得

式中kij（i，j＝1，…，6）為復合材料軸橫截面的剛度系數，具體表達式由于篇幅所限，不再列出。

通過比較可以發現，在全部的36個剛度系數kij（i，j＝1，…，6）中，其中的16個剛度系數kij（i，j＝1，2，…，4）的表達式與文獻［14］中不計剪切變形的復合材料薄壁梁的剛度系數是一致的，而其余的20個剛度系數kij（i＝1，2，…，6；j＝5，6），kji（i＝1，2，…，4；j＝5，6）是由于計及剪切變形新增加的剛度系數。

假設復合材料軸具有周向均勻剛度配置（CUS）構型［12］，即滿足

θ（y）＝θ（－y），θ（z）＝θ（－z）

式中θ表示由正向s軸進行度量的纖維鋪層角。

將方程（9）和（10）代入（8），由 Hamilton原理可以導出運動方程如下

如果在方程（13）中的第二、三、五、六個方程中，令U″2－θ′y＝0，U″3－θ′z＝0，U′2－θy＝0，U′3－θz＝0，則可以導出不考慮剪切的復合材料軸的運動方程，此處不再寫出。

方程（13）中的第二、三、五、六個方程，構成彎－剪耦合振動方程

值得注意的是，在方程（14）中除了存在彈性變性產生的耦合，也存在由于剛性旋轉產生的耦合，包括與轉速一次方Ω相關的項，以及與轉速平方Ω2相關的項。此外，在方程（13）中的第一和第四兩個方程構成拉－扭耦合振動方程，此類耦合問題不在本文中進行討論。

2 自由振動方程求解

假定彎曲位移U2（x，t），U3（x，t）和轉角θy（x，t），θz（x，t）具有下列形式

式中U2j（t），U3j（t），Θyj（t）和Θzj（t）表示廣義坐標；αj（x）和ψj（x）表示軸的振型函數。

將方程（15）代入（14），采用Galerkin近似求解方法，得

式中M為質量矩陣，C為旋轉陀螺效應產生的阻尼矩陣，K為由彈性變形剛度與旋轉引起的剛度疊加而成剛度矩陣

其中

廣義坐標矢量表示為

運動方程（16）可以進一步化為廣義特征值問題

式中

其中，0和I分別為零矩陣和單位矩陣。

3 數值結果與討論

3.1 模型驗證

算例1 為了檢驗本文建立的復合材料軸模型及其近似計算方法的正確性，下面首先針對一個石墨／環氧復合材料軸的固有頻率特性，進行了數值方法的收斂性檢驗。軸的幾何尺寸和材料參數分別取［12］：長度L＝2.023m，截面半徑r＝0.127m，單層厚度h＝0.063 5mm，截面鋪層方式為［θ］6；復合材料的性能參數為E1＝206.8GPa，E2＝E3＝5.17 GPa，G12＝3.1GPa，G23＝G13＝2.55GPa，ν21＝ν31＝0.006 25，ν32＝0.25，密度ρ＝1 528.15kg／m3。引入了標準化因子ω0＝138.85rad／s（表示非旋轉軸在纖維鋪層角θ＝0°的第一階固有頻率），則無量綱的固有頻率和轉速分別為ω＊＝ω×2π／ω0，Ω＊＝Ω×2π／（60ω0）。其中固有頻率ω和轉速Ω的單位分別為Hz和r／min。

表1表示兩端簡支軸的前六階固有頻率隨振型函數個數N的變化情況，結果表明，本文提出的近似計算方法具有很好的收斂性，例如，為了獲得前三階固有頻率，振型函數的個數只需取N＝6就可得到較高精度的結果（注：其中“－”表示大于保留最大模態個數N的高階頻率，由于已經超出了方程組（16）的階數，所以沒有結果顯示）。

表1 模態個數對于固有頻率的影響（Ω＊＝0，θ＝30°）Tab.1 Effect of model number Non natural frequencies（Ω＊＝0andθ＝30°）

表2表示不計剪切變形的懸臂復合材料軸的固有頻率計算結果的對比，其中計算參數和無量綱化方法同文獻［18］。由表2可以看出，本文結果與文獻［18］結果符合得很好。

表2 不計剪切的懸臂復合材料軸的固有頻率結果對比Tab.2 Comparison of the natural frequencies of a cantilevercomposite shaft without shear deformation

圖2表示兩端簡支復合材料軸的固有頻率隨轉速的變化曲線，其中考慮了剪切變形的影響。從圖2可以看到，由于旋轉的復合材料軸存在陀螺效應，所以固有頻率在轉速Ω＊≠0展示了分叉現象。其次，從圖2還可以看到，針對四種不同的鋪層角，本文結果與文獻［12］的結果非常相近。

圖2 不同鋪層角的復合材料軸固有頻率隨轉速的變化曲線Fig.2 The natural frequency of a simply supported composite shaft versus rotating speed for different ply angles

3.2 自由振動與穩定性分析

算例2 在下面的數值計算中，選取復合材料軸的幾何尺寸為：長度L＝1.67m，截面半徑r＝0.063 5m，單層厚度h＝0.132 1mm，截面鋪層方式為［θ］8，軸的兩端具有簡支邊界條件。彎曲振型和扭轉振型函數的具體表達式見文獻［8］和［9］。復合材料的性能參數如表3所示。

表3 復合材料性能參數Tab.3 Mechanical properties of composite material

圖3表示復合材料軸的第一階固有頻率隨轉速的變化曲線。結果表明，當轉速為0時，每個鋪層角只對應于單獨的一個固有頻率，這是由于對于圓形截面軸而言，它在豎直（Z軸向）和水平方向（Y軸向）彎曲模態頻率是相同的。隨著轉速的增加，由于旋轉陀螺效應，固有頻率曲線分叉為上下兩支，上支隨固有頻率增加而增加，下支隨固有頻率增加而減少。分別稱為上固有頻率（Upper natural frequency，UNF）和下固有頻率（Lower natural frequency，LNF），它們分別對應于慣性坐標系下的正進動頻率和反進動頻率。當LNF隨轉速的增加而減小至0時，所對應的轉速稱為臨界轉速Ω＊c；從圖中還可以看出，無論是對于一階模態還是二階模態，臨界轉速都隨著鋪層角的增加而增加，最大的臨界轉速對應的鋪層角θ為90°。這主要是由于當纖維鋪設角沿軸的軸向鋪層時，軸的彎曲剛度變為最大。

表4給出不計剪切和計入剪切時的復合材料軸的第一階臨界轉速和第二階臨界轉速的計算結果，由此可以看出，不計剪切的臨界轉速高于計及剪切的臨界轉速，并且兩者的差別隨著鋪層角的增加而增加；在相同鋪層角下，計及剪切和不計剪切的第二階臨界轉速的相對誤差要高于相應的第一階臨界轉速的相對誤差。例如在鋪層角θ為90°時，這一誤差由2.40%增加到9.43%。因此，在預測鋪層角較大時的高階臨界轉速時，剪切變形的影響尤其就顯得十分重要。

圖3 不同鋪層角復合材料軸的第一階固有頻率隨轉速變化曲線Fig.3 The first natural frequency of a composite shaft versus rotating speed for different ply angles

表4 計及剪切與不計剪切復合材料軸的臨界轉速Tab.4 The critical speed of a composite shaft with and without shear deformation

圖4和5分別表示靜止和旋轉復合材料軸的前兩階固有頻率隨鋪層角的變化曲線。結果表明，在靜止狀態下，軸的第一階、第二階模態所對應的UNF曲線和LNF曲線是重合的，它們都隨著鋪層角的增加而增加，這與前面的結論是一致的。另一方面，相對于第一階固有頻率，第二階固有頻率隨鋪層角的增加而增加的趨勢更為明顯；在旋轉狀態下的前兩階固有頻率隨鋪層角的變化趨勢和在靜止狀態的變化趨勢是相同的，只不過此的UNF曲線和LNF曲線是分開的。

圖4 復合材料軸前兩階固有頻率隨鋪層角變化曲線（Ω＝0）Fig.4 The first two natural frequency of a simply supported composite shaft versus ply angle（Ω＝0）

圖6 復合材料軸的前兩階固有頻率隨長徑比的變化曲線（Ω＝2 000r·min－1，θ＝90°）Fig.6 The first two natural frequencies of a composite shaft versus ratio of length over radius（Ω＝2 000 r·min－1，θ＝90°）

圖6和7分別表示復合材料軸的前兩階旋轉固有頻率隨長徑比和徑厚比的變化曲線。由圖6可以看出，固有頻率隨著長徑比的增加而減小；由圖7可以看出，固有頻率隨著徑厚比的增加而增加。

圖7 復合材料軸的前兩階固有頻率隨徑厚比的變化曲線（Ω＝2 000r·min－1，θ＝90°）Fig.7 The first two natural frequencies of a composite shaft versus ratio of radius over thickness（Ω＝2 000r·min－1，θ＝90°）

圖8和9分別表示具有不同鋪層角的復合材料軸的前兩階臨界轉速隨長徑比和徑厚比的變化曲線。由圖8和9可以看出，臨界轉速隨著長徑比的增加表現為減小的趨勢，而隨著徑厚比的增加表現為增加的趨勢。

圖8 復合材料軸前兩階臨界轉速隨長徑比的變化曲線Fig.8 The first two critical speeds of a composite shaft versus ratio of length over radius

圖9 復合材料軸兩階臨界轉速隨徑厚比的變化曲線Fig.9 The first two critical speeds of a composite shaft versus ratio of radius over thickness

4 結 論

本文基于復合材料薄壁梁理論建立了旋轉軸的動力學模型。結構的位移場采用VAM模型進行描述，并且在其中引入了橫向剪切變形的影響；采用Hamilton原理導出旋轉軸的控制方程并采用Galerkin法對數學模型進行數值求解，獲得旋轉軸耦合振動的固有頻率及其臨界轉速的近似分析解。本文模型的正確性通過收斂性分析以及與文獻結果對比，得到了驗證。采用本文建立的模型與計算方法，研究了旋轉復合材料軸的自由振動與動力穩定性。數值結果表明：（1）當復合材料軸的鋪層角或者徑厚比增加時，固有頻率隨之增加；當復合材料軸的長徑比增加時，固有頻率隨之減小，在軸的剛性旋轉作用下，固有頻率分叉為UNF和LNF曲線，其中UNF隨轉速的增加而增加，LNF隨轉速的增加而降低；（2）不計剪切變形的臨界轉速高于計入剪切變形的臨界轉速，并且兩者的差別隨著鋪層角的增加而增加，在相同鋪層角下，計及剪切和不計剪切的低階臨界轉速的之間相對誤差要高于相應的高階臨界轉速之間的相對誤差；（3）復合材料軸的臨界轉速隨著鋪層角或者徑厚比的增加而增加；隨著長徑比的增加而減小。本文提出的模型可以為復合材料軸的動力學研究與設計，提供一種可供選擇的分析理論與方法。

［1］ 王震鳴，杜善義，張衡，等.復合材料及其結構的力學、設計、應用和評價（第二冊）［M］.哈爾濱：哈爾濱工業大學出版社，1998.Wang Zhenming，Du Shanyi，Zhang Heng，et al.Mechanics，design，application and evaluation of composites and structures（partⅡ）［M］.Haibin：Harbin Institute of Technology Press，1998.

［2］ Hetherington E L，Kraus R E，Darlow M S.Demonstration of a super critical composite helicopter power transmission shaft［J］.Journal of American Helicopter Society，1990，35（1）：23—28.

［3］ Faust H，Mack J，Spencer B.A composite rotor shaft for the Chinook［J］.Journal of the American Helicopter Society，1984，29：54—58.

［4］ dos Reis H L M，Goldman R B，Verstrate P H.Thinwalled laminated composite cylindrical tubes：Part IIICritical speed analysis［J］.Journal of Composites Technology and Research，1987，9：58—62.

［5］ Zinberg H，Symonds M F.The development of an advanced composite tail rotor drive shaft［A］.Presented at the 26th Annual National Forum of the American Helicopter Society［C］.June，Washington，DC，1970.

［6］ Singh S P，Gupta K.Composite shaft rotordynamic analysis using a layer wise theory［J］.Journal of Sound and Vibration，1996，191（5）：739—756.

［7］ Kim C D，Bert C W.Critical speed analysis of laminated composite hollow drive shafts［J］.Composites Engineering，1993，3：633—643.

［8］ Bert C W，Kim C D.Whirling of composite－material drivshafts including bending－twisting coupling and transverse shear deformation［J］.Journal of Vibration and Acoustics，1995，117：17—21.

［9］ Kim W，Argento A，Scott R A.Free vibration of a rotating tapered composite Timoshenko shaft［J］.Journal of Sound and Vibration，1999，226（1）：125—147.

［10］Chang M Y，Chen J K，Chang C Y.A simple spinning laminated composite shaft model［J］.Int.Journal of Solids and Structures，2004，41：637—662.

［11］Song O，Librescu L.Anisotropy and structural coupling on vibration and instability of spinning thinwalled beams［J］.Journal of Sound and Vibration，1997，204：477—494.

［12］Song O，Jeong N－H，Librescu L.Implication of conservative and gyroscopic forces on vibration and stability of an elastically tailored rotating shaft modeled as a composite thin－walled beam ［J］.J.Acoust.Soc.Am.，2001，109（31）：972—981.

［13］Rehfield L W.Design analysis methodology for composite rotor blades［A］.Proc.Seventh DoD／NASA Conf.on Fibrous Composites in Structural Design［C］.AFWAL－TR－85－3094，1985：1—15.

［14］Berdichevsky V，Armanios E，Badir A.Theory of anistropic thin－walled closed－cross－section beams［J］.Composites Engineering，1992，2：411—432.

［15］Cesnik C E S，Shin S J.On the modeling of integrally actuated helicopter blades［J］.International Journal of Solids and Structures，2001，38：1 765—1 789.

［16］Park J S，Kim J H.Design and aeroelastic analysis of active twist rotor blades incorporating single crystal macro fiber composite actuators［J］.Composites：Part B，2008，39：1 011—1 025.

［17］Song O，Librescu L.Structural modeling and free vibration analysis of rotating composite thin－walled beams［J］.Journal of the American Helicopter Society，1997：358—369.

［18］Banerjee J R，Su H.Development of a dynamic stiffness matrix for free vibration analysis of spinning beams［J］.Comput.Struct.，2004，82：2 189—2 197.