旋轉(zhuǎn)復(fù)合材料薄壁軸的自由振動與穩(wěn)定性＊

任勇生，代其義，張興琦

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

引 言

在包括能源、航空和汽車在內(nèi)的許多工業(yè)技術(shù)領(lǐng)域內(nèi)，廣泛地存在著一類旋轉(zhuǎn)柔性細長結(jié)構(gòu)，它的轉(zhuǎn)動角速度矢量與其結(jié)構(gòu)的縱軸保持平行，在能量傳動過程中發(fā)揮著無可替代的重要作用，力學(xué)上通常將其稱之為旋轉(zhuǎn)軸。例如燃氣輪發(fā)電機傳動軸、航空發(fā)動機傳動軸以及汽車后輪驅(qū)動軸等等。

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

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

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

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

1 振動微分方程的建立

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

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

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

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

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

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

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

定義薄壁軸內(nèi)力

薄壁軸內(nèi)力－應(yīng)變本構(gòu)關(guān)系

其中

為了導(dǎo)出復(fù)合材料軸的振動方程，利用Hamilton原理

式中U和T分別為應(yīng)變能和動能，分別由下式確定

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

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

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

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

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

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

假設(shè)復(fù)合材料軸具有周向均勻剛度配置（CUS）構(gòu)型［12］，即滿足

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

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

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

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

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

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

2 自由振動方程求解

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

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

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

式中M為質(zhì)量矩陣，C為旋轉(zhuǎn)陀螺效應(yīng)產(chǎn)生的阻尼矩陣，K為由彈性變形剛度與旋轉(zhuǎn)引起的剛度疊加而成剛度矩陣

其中

廣義坐標矢量表示為

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

式中

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

3 數(shù)值結(jié)果與討論

3.1 模型驗證

算例1 為了檢驗本文建立的復(fù)合材料軸模型及其近似計算方法的正確性，下面首先針對一個石墨／環(huán)氧復(fù)合材料軸的固有頻率特性，進行了數(shù)值方法的收斂性檢驗。軸的幾何尺寸和材料參數(shù)分別?。?2］：長度L＝2.023m，截面半徑r＝0.127m，單層厚度h＝0.063 5mm，截面鋪層方式為［θ］6；復(fù)合材料的性能參數(shù)為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（表示非旋轉(zhuǎn)軸在纖維鋪層角θ＝0°的第一階固有頻率），則無量綱的固有頻率和轉(zhuǎn)速分別為ω＊＝ω×2π／ω0，Ω＊＝Ω×2π／（60ω0）。其中固有頻率ω和轉(zhuǎn)速Ω的單位分別為Hz和r／min。

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

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

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

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

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

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

3.2 自由振動與穩(wěn)定性分析

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

表3 復(fù)合材料性能參數(shù)Tab.3 Mechanical properties of composite material

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

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

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

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

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

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

圖6 復(fù)合材料軸的前兩階固有頻率隨長徑比的變化曲線（Ω＝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分別表示復(fù)合材料軸的前兩階旋轉(zhuǎn)固有頻率隨長徑比和徑厚比的變化曲線。由圖6可以看出，固有頻率隨著長徑比的增加而減?。挥蓤D7可以看出，固有頻率隨著徑厚比的增加而增加。

圖7 復(fù)合材料軸的前兩階固有頻率隨徑厚比的變化曲線（Ω＝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分別表示具有不同鋪層角的復(fù)合材料軸的前兩階臨界轉(zhuǎn)速隨長徑比和徑厚比的變化曲線。由圖8和9可以看出，臨界轉(zhuǎn)速隨著長徑比的增加表現(xiàn)為減小的趨勢，而隨著徑厚比的增加表現(xiàn)為增加的趨勢。

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

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

4 結(jié) 論

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

［1］ 王震鳴，杜善義，張衡，等.復(fù)合材料及其結(jié)構(gòu)的力學(xué)、設(shè)計、應(yīng)用和評價（第二冊）［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，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.