周健斌，章俊杰，孟 光

(1.中國商用飛機有限責(zé)任公司 上海飛機設(shè)計研究院，上海 200232;2.上海交通大學(xué) 機械系統(tǒng)與振動國家重點實驗室，上海 200240)

現(xiàn)代大型客機設(shè)計廣泛采用機翼下吊發(fā)動機的布局形式。隨著飛機設(shè)計技術(shù)的發(fā)展，現(xiàn)代大型客機朝著加大機翼展弦比、采用輕質(zhì)材料、減小結(jié)構(gòu)剛度、采用高涵道比大推力發(fā)動機的方向發(fā)展，因此翼下吊發(fā)動機不僅影響飛機的氣動力布局，發(fā)動機轉(zhuǎn)子產(chǎn)生的陀螺效應(yīng)對飛機結(jié)構(gòu)動力學(xué)特性以及氣動彈性特性的影響也日益突出。王彬文［1］等用有限元方法研究了轉(zhuǎn)子陀螺效應(yīng)對帶有翼吊發(fā)動機系統(tǒng)振動特性的影響，結(jié)果表明陀螺效應(yīng)對其固有特性的影響是不可忽略的。

具有陀螺效應(yīng)的彈性體問題可歸結(jié)為陀螺彈性系統(tǒng)問題。Eleuterio［2-6］和 Hughes［7］首次提出了陀螺彈性系統(tǒng)問題并針對自旋航天器研究了一類具有分布陀螺力矩的彈性系統(tǒng)的動力學(xué)特性。Peck等［8，9］在此基礎(chǔ)上提出了運用嵌入在結(jié)構(gòu)中的角動量來實現(xiàn)結(jié)構(gòu)自適應(yīng)調(diào)諧，如頻移、模態(tài)耦合和解耦、相位調(diào)整等。針對一類帶有旋轉(zhuǎn)部件的機械臂，Yamanaka［10-11］分析了端部帶有轉(zhuǎn)子的均質(zhì)梁的動力學(xué)和穩(wěn)定性問題，其轉(zhuǎn)子轉(zhuǎn)軸沿梁的軸向。Li等［12］研究了端部帶有任意方向旋轉(zhuǎn)軸的轉(zhuǎn)子的柔性連接的穩(wěn)定性問題。

本文從陀螺彈性系統(tǒng)理論出發(fā)，以帶有橫向轉(zhuǎn)子的Bernoulli-Euler懸臂梁為研究對象，運用Hamilton變分原理建立了計及發(fā)動機陀螺效應(yīng)的機翼彎扭耦合動力學(xué)模型，分析了系統(tǒng)的動力學(xué)特性，研究了發(fā)動機陀螺效應(yīng)對系統(tǒng)固有特性的影響。

1 陀螺效應(yīng)的耦合作用

如圖1所示，轉(zhuǎn)子垂直于梁展向安裝，設(shè)轉(zhuǎn)子自轉(zhuǎn)角速度為Ω，極慣性矩為Jp，則轉(zhuǎn)子對其質(zhì)心的動量矩為:

當(dāng)梁發(fā)生y方向(面內(nèi))的彎曲振動以及扭轉(zhuǎn)振動時，轉(zhuǎn)子轉(zhuǎn)動方向?qū)㈦S梁振動而發(fā)生改變，此時轉(zhuǎn)子做繞y方向的橢圓進動。根據(jù)轉(zhuǎn)子動力學(xué)理論，該橢圓進動由正進動和反進動合成。設(shè)轉(zhuǎn)子繞y方向的進動角速度為ω，則轉(zhuǎn)子產(chǎn)生陀螺力矩:

其方向垂直于Ω與ω組成的平面。設(shè)Ω與ω夾角為φ(φ?1)，則陀螺力矩大小可表示為:

該力矩與φ成正比，相當(dāng)于彈性力矩。在正進動情況下(如圖1)，陀螺力矩作用使得梁的剛度提高;反之，在反進動情況下，梁的剛度減小。由此，陀螺力矩使得梁的面內(nèi)彎曲模態(tài)和扭轉(zhuǎn)模態(tài)發(fā)生耦合，并影響系統(tǒng)的剛度。

圖1 陀螺力矩耦合作用示意圖(正進動)Fig.1 Illustration of the mechanism of the gyroscopic coupling effect

2 動力學(xué)模型

為研究發(fā)動機陀螺效應(yīng)對機翼動力學(xué)特性的影響，將機翼-發(fā)動機系統(tǒng)簡化為一端固支的彎扭耦合Bernoulli-Euler懸臂梁結(jié)構(gòu)，忽略發(fā)動機與機翼間的連接剛度，如圖2(a)所示。機翼翼展為l，半弦長為b，右手坐標(biāo)系o-xyz原點o位于機翼根部，ox軸與機翼剛軸重合，oy軸沿飛機航向。發(fā)動機重心位置為(xs，ys，zs)，機翼重心軸在-y方向的距離為e。圖2(b)所示為變形后x=xs處機翼弦向截面示意圖。o'為機翼剛軸與截面交點，截面位移為v、w，截面繞o'轉(zhuǎn)角為 θ。Cw為截面重心，發(fā)動機重心位置Cs，發(fā)動機重量為ms，轉(zhuǎn)子轉(zhuǎn)速為Ω，轉(zhuǎn)動矢量方向與剛軸垂直。

圖2 機翼-發(fā)動機系統(tǒng)示意圖Fig.2 Illustration of the wing-engine system

由Hamilton變分原理

可推導(dǎo)得系統(tǒng)動力學(xué)方程，其中 δU、δTw、δTs和 δW分別為機翼勢能、機翼動能、發(fā)動機動能和非保守外力做功的變分。由彎扭耦合懸臂梁動力學(xué)分析，機翼勢能的變分為:

其中，EI2為梁的面外彎曲剛度，EI3為梁的面內(nèi)彎曲剛度，GJ為梁的扭轉(zhuǎn)剛度。

機翼動能的變分可表示為:

其中，mw為機翼單位長度質(zhì)量，e為機翼截面重心與剛軸的距離(重心在剛軸前為正)，σ為機翼單位長度截面繞剛軸的慣性半徑。

對發(fā)動機動能分析可知，發(fā)動機的動能的變分由兩部分組成，δTs=δTm+δTr，其中 δTm為發(fā)動機集總質(zhì)量的動能的變分，δTr為發(fā)動機轉(zhuǎn)子的轉(zhuǎn)動動能的變分。

由轉(zhuǎn)子動力學(xué)理論，當(dāng)機翼變形時，向量Ω的偏轉(zhuǎn)角在oxy和oyz平面的投影分別為v'和θ，由此，發(fā)動機轉(zhuǎn)子的轉(zhuǎn)動的動能可表示為:

其中，Jp為轉(zhuǎn)子繞其轉(zhuǎn)軸的極轉(zhuǎn)動慣量，Jd為轉(zhuǎn)子的赤道轉(zhuǎn)動慣量。對Tr求變分可得:

將式(5)～(7)和(9)代入式(4)可得如下計及發(fā)動機陀螺效應(yīng)的系統(tǒng)動力學(xué)方程:

式中，δD()為 Dirac函數(shù)?！?”表示對X求導(dǎo)，“·”表示對τ求導(dǎo)。

3 動力學(xué)特性分析

采用extended Galerkin［13］方法對系統(tǒng)動力學(xué)方程進行離散化，令:

選取如下滿足邊界條件的函數(shù)作為基函數(shù):

代入方程(13)～方程(15)，積分可得離散后的系統(tǒng)方程:

為驗證系統(tǒng)模型，選取 Goland［14］機翼模型參數(shù)EI2=9.76 ×106N·m2，EI3=9.76 ×106N·m2，GJ=9.88 ×105N·m2，mw=35.75 kg/m，σ =0.49 m，e=0.18 m，l=6.10 m，計算發(fā)動機位于機翼根部且 Ω=0 rad/s時系統(tǒng)前五階固有頻率，結(jié)果如表1所示。根據(jù)梁的動力學(xué)特點，取N1=5、N2=5、N3=3。通過對比可知，本文計算結(jié)果與文獻(xiàn)相吻合，數(shù)學(xué)模型及計算方法具有較高的精度。

表1 Goland梁前五階固有頻率Tab.1 Modal frequency of the first 5 modes for the Goland wing

3.1 轉(zhuǎn)速對系統(tǒng)固有特性的影響

表2 轉(zhuǎn)子結(jié)構(gòu)參數(shù)Tab.2 Structural parameters of the engine rotor

圖3 Goland機翼固有頻率和阻尼比隨發(fā)動機無量綱慣性矩變化曲線Fig.3 Modal frequency and damping ratio vs.nondimensional rotor moment of inertia of the Goland wing

由圖3可知，隨著發(fā)動機轉(zhuǎn)子慣性矩的增大，在發(fā)動機陀螺力矩作用下，當(dāng)系統(tǒng)二階扭轉(zhuǎn)模態(tài)與面內(nèi)二階彎曲模態(tài)相鄰時，該兩階模態(tài)相互耦合，系統(tǒng)二階扭轉(zhuǎn)模態(tài)頻率增大、面內(nèi)二階彎曲模態(tài)頻率減小，即在陀螺力矩作用下，系統(tǒng)二階扭轉(zhuǎn)模態(tài)剛度增大，而二階彎曲剛度減小;此時系統(tǒng)特征值為純虛數(shù)，各階振型有固定節(jié)線，即發(fā)動機陀螺力矩影響系統(tǒng)的剛度。當(dāng)轉(zhuǎn)子慣性矩增大至=1.56時，上述耦合模態(tài)對應(yīng)的特征值為共軛復(fù)數(shù)重根，模態(tài)頻率重合，模態(tài)阻尼分岔并出現(xiàn)負(fù)阻尼，系統(tǒng)處于動態(tài)失穩(wěn)狀態(tài);如圖4所示為=2.52時系統(tǒng)前六階振型，此時系統(tǒng)第四階模態(tài)與第五階模態(tài)相互耦合，具有相同的振型曲線，且振型無固定的節(jié)線。重合后該模態(tài)頻率隨著發(fā)動機轉(zhuǎn)速的增大而加速降低=5.55時，該耦合模態(tài)出現(xiàn)過阻尼，模態(tài)頻率為零，相應(yīng)模態(tài)阻尼進一步分岔。

圖4 H=2.52時系統(tǒng)前六階振型圖Fig.4 Mode shapes of the first 6 modes with=2.52

3.2 臨界轉(zhuǎn)速

以上如圖5所示為臨界無量綱慣性矩及=0時系統(tǒng)模態(tài)頻率ω/ω0隨機翼剛軸與重心軸距離e/l的變化曲線。由圖可知=0時系統(tǒng)面外彎曲模態(tài)及扭轉(zhuǎn)模態(tài)的頻率隨著e/l的增大而變化，當(dāng)其中某一階頻率曲線接近并與面內(nèi)彎曲模態(tài)頻率曲線相交時，將相交點臨近區(qū)域稱為頻率相交影響區(qū)域，如圖5中陰影區(qū)域Ⅰ和Ⅱ所示。在頻率相交影響區(qū)域內(nèi)，隨著發(fā)動機轉(zhuǎn)子轉(zhuǎn)速的提高，發(fā)動機陀螺效應(yīng)使得頻率相交模態(tài)參與耦合，系統(tǒng)動態(tài)失穩(wěn)，此時系統(tǒng)臨界轉(zhuǎn)速快速降低，且在頻率相交點臨界頻率達(dá)到最小值，如圖6所示為e/l=0.01時，系統(tǒng)模態(tài)頻率隨發(fā)動機轉(zhuǎn)子轉(zhuǎn)速變化Campbell圖，從圖中可以看出，由于系統(tǒng)結(jié)構(gòu)參數(shù)使得系統(tǒng)面內(nèi)一階彎曲頻率和面外一階頻率耦合，在陀螺效應(yīng)的作用下，系統(tǒng)第一、二階頻率在=0.02 時即出現(xiàn)共軛特征值，導(dǎo)致動態(tài)失穩(wěn)。在非頻率相交影響區(qū)域內(nèi)，系統(tǒng)臨界阻尼不隨e/l變化。

如圖7所示為臨界無量綱慣性矩及=0時系統(tǒng)模態(tài)頻率ω/ω0隨機翼慣性半徑σ/l的變化曲線。與前述情況類似，由于慣性半徑的變化系統(tǒng)在0≤≤0.15范圍內(nèi)出現(xiàn)頻率相交影響區(qū)域Ⅰ、Ⅱ、Ⅲ和Ⅳ，在這些區(qū)域內(nèi)系統(tǒng)臨界轉(zhuǎn)速快速下降，而在頻率相交影響區(qū)域內(nèi)，系統(tǒng)臨界轉(zhuǎn)速隨慣性半徑的增大而增大。

圖7 Goland機翼臨界無量綱慣性矩隨變化曲線Fig.7 Critical nondimensional rotor moment of inertia vs. σ of the Goland wing

4 結(jié)論

本文運用Hamilton原理建立了計及發(fā)動機陀螺效應(yīng)的機翼動力學(xué)數(shù)學(xué)模型，通過系統(tǒng)特性分析、研究得到如下結(jié)論:

(1)受發(fā)動機陀螺力矩的耦合作用，隨著發(fā)動機陀螺慣性矩的增大，機翼面內(nèi)彎曲模態(tài)頻率增大，扭轉(zhuǎn)模態(tài)頻率降低，且模態(tài)振型相耦合。

(2)當(dāng)面內(nèi)彎曲模態(tài)與扭轉(zhuǎn)模態(tài)相鄰時，且發(fā)動機陀螺無量綱慣性矩大于臨界穩(wěn)定值時，系統(tǒng)特征值出現(xiàn)共軛復(fù)重根，參與耦合的模態(tài)頻率重合，系統(tǒng)動態(tài)失穩(wěn)。

(3)在非頻率相交影響區(qū)域內(nèi)，系統(tǒng)動態(tài)臨界穩(wěn)定無量綱慣性矩隨機翼慣性半徑σ的增大而增大，而剛軸與重心軸的距離e對無影響。

［1］王彬文，孫俠生，齊丕騫.轉(zhuǎn)子陀螺效應(yīng)對系統(tǒng)振動特性的影響［J］.振動工程學(xué)報，2010，23(增刊):63-67.

［2］ D'Eleuterio G M T.On the theory of gyroelasticity［J］.Journal of Applied Mechanics，Transactions ASME，1988，55(2):488-489.

［3］D'Eleuterio G M T，Hughes P C.Dynamics of gyroelastic spacecraft［J］.Journal of Guidance，Control，and Dynamics，1987，10(4):401-405.

［4］D'Eleuterio G M T，Hughes P C. General motion of gyroelastic vehicles in terms of constrained modes［C］.1985，384-390.

［5］ D'Eleuterio G M T，Hughes P C.Dynamics of gyroelastic continua［J］.Journal of Applied Mechanics，Transactions ASME，1984，51(2):415-422.

［6］ D'Eleuterio G M T.Dynamics of gyroelastic vehicles［R］.UTIAS Report(University of Toronto，Institute for Aerospace Studies)，1986.

［7］Hughes P C，D'Eleuterio G M T.Modal parameter analysis of gyroelastic continua［J］.Journal of Applied Mechanics，1986，53(4):918-924.

［8］Peck M A，Cavender A R.Structural tuning through embedded angular momentum［C］.44th AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Conference，2003，2:1462-1469.

［9］Peck M A，Cavender A R.Practicable gyroelastic technology［C］.27th Annual AAS Rocky Mountain Guidance and Control Conference，2004，118:239-253.

［10］Yamanaka K，Heppler G R，Huseyin K.On the dynamics and stability of a beam with a tip rotor［C］.the 35th AIAA/ASME/ASCE/AHS/ASC Structures， Structural Dynamics，and Materials Conference，1994，2:1031-1038.

［11］ Yamanaka K， Heppler G R， Huseyin K. Stability of gyroelastic beams［J］.AIAA Journal，1996，34(6):1270-1278.

［12］ Li L，Heppler G R，Huseyin K.Stability of a flexible link with an arbitrarily oriented tip rotor and a conservative tip load［C］.2000，2:1472-1477.

［13］ Fletcher C.Computational Galerkin methods［M］.Springer-Verlag New York，1984.

［14］ Goland M.The flutter of a uniform cantilever wing［J］.Journal of Applied Mechanics，1945，12(4):197-208.

［15］ Banerjee J，Su H.Dynamic stiffness formulation for beams undergoing free transverse and lateral vibration with torsional coupling［C］.44th AIAA/ASME/ASCE/AHS Structures，Structural Dynamics，and Materials Conference，2003.