何 鵬，劉占生，劉鎮(zhèn)星

(哈爾濱工業(yè)大學 能源學院，哈爾濱 150001)

汽輪機和燃氣輪機在工作過程中，受到工作介質(zhì)及環(huán)境的影響，轉(zhuǎn)子上存在一定的溫度分布。在溫度作用下，轉(zhuǎn)子發(fā)生熱彎曲，同時轉(zhuǎn)子材料的力學性能有所變化，影響轉(zhuǎn)子的振動特性。

在溫度作用下，轉(zhuǎn)子材料發(fā)生熱膨脹，引起轉(zhuǎn)子熱變形。轉(zhuǎn)子周向溫度分布不均勻時，將產(chǎn)生熱彎曲，影響轉(zhuǎn)子振動狀態(tài)。轉(zhuǎn)子的熱彎曲振動問題在工程中時常出現(xiàn)［1－3］。轉(zhuǎn)子熱彎曲成因可以分為周向介質(zhì)及上下缸溫差［4］、動靜碰摩產(chǎn)生的溫差［5－8］和轉(zhuǎn)子在滑動軸承中渦動引起的Morton效應［9－10］。對于轉(zhuǎn)子熱彎曲所產(chǎn)生的振動問題已經(jīng)有深入研究［11－12］，并且已經(jīng)形成了有效的分析診斷方法和處理措施［13－15］。

溫度不僅能引發(fā)轉(zhuǎn)子熱變形，還會引起轉(zhuǎn)子材料楊氏模量的改變。不同材料的楊氏模量隨溫度變化規(guī)律不同，文獻［16］通過測量不同溫度下振動頻率的變化，得到了含鎳鋼材的楊氏模量隨溫度變化曲線。文獻［17］基于固體力學理論，推導了常用合金材料楊氏模量隨溫度的變化規(guī)律。由溫度引起材料楊氏模量的變化對轉(zhuǎn)子振動特性也有影響。晏水平等［18］引入剪切模量隨溫度的變化規(guī)律，分析了軸向溫度分布對汽輪機轉(zhuǎn)子軸系扭轉(zhuǎn)振動的影響。朱向哲等［19－20］通過引入熱應力效應，考慮溫度作用下轉(zhuǎn)子熱變形引起剛度矩陣變化，分析了溫度分布對轉(zhuǎn)子振動特性的影響。劉少權等［21］引入楊氏模量隨溫度的變化，采用有限元軟件分析了溫度場對燃氣輪機拉桿轉(zhuǎn)子臨界轉(zhuǎn)速的影響。雖然對于溫度導致材料楊氏模量變化引起的振動問題已經(jīng)進行了一些研究，目前尚未有考慮溫度引起楊氏模量變化的轉(zhuǎn)子有限元模型，以分析在不同軸向溫度分布下轉(zhuǎn)子振動特性。

傳統(tǒng)的轉(zhuǎn)子動力學有限元模型［22］無法考慮轉(zhuǎn)子軸向溫度分布的影響，本文主要研究考慮軸向溫度分布的有限元建模。首先擬合了楊氏模量隨溫度的變化關系，結(jié)合有限元方法，推導了考慮單元內(nèi)溫度分布的有限元剛度矩陣，建立了溫度分布對轉(zhuǎn)子振動特性的影響模型。通過數(shù)值仿真分析了在溫度均勻分布、軸向線性分布和軸向二次分布情況下轉(zhuǎn)子振動特性的變化。

1 楊氏模量隨溫度的變化規(guī)律

材料的楊氏模量是材料原子間作用力的體現(xiàn)，反映了材料應力與應變的比例關系，是材料重要的力學性能參數(shù)。對于某一材料，楊氏模量只與材料的溫度有關。

不同材料的楊氏模量隨溫度的變化規(guī)律不同。文獻［18］中，轉(zhuǎn)子材料剪切模量隨溫度變化規(guī)律為G=G0(1－c(t－t0)/100)，其中 c為常數(shù)。在文獻［19－20］中，楊氏模量隨溫度的變化規(guī)律為E(t)=207.8－0.069t－3.8 ×10－5t2。在文獻［21］中，給出了所計算壓氣機和渦輪轉(zhuǎn)子材料在不同溫度下的楊氏模量值。通過曲線擬合，可以得到文獻［21］中壓氣機轉(zhuǎn)子材料楊氏模量隨溫度的變化規(guī)律為E(t)=E0(－3.862×10－7t2－1.199 ×10－4t+1.002)，其中擬合的均方根誤差為0.001923;渦輪轉(zhuǎn)子材料楊氏模量隨溫度變化規(guī)律為E(t)=E0(－5.213×10－7t2－8.347×10－5t+1.002)，其中擬合的均方根誤差為0.003989。在文獻［16］中，作者通過測量不同溫度下叉狀物的振動頻率，得到了不同溫度下含鎳鋼材楊氏模量隨溫度的變化曲線。通過曲線擬合，可得該材料楊氏模量隨溫度變化規(guī)律為:

可以看出，不同材料楊氏模量隨溫度的變化規(guī)律略有不同，但大多轉(zhuǎn)子材料的楊氏模量隨溫度的變化規(guī)律可以用二次多項式形式表示:

其中，對于確定的材料，γ，β和α為常數(shù)。

2 考慮楊氏模量變化的轉(zhuǎn)子模型

根據(jù)楊氏模量隨溫度的變化規(guī)律，結(jié)合轉(zhuǎn)子軸向溫度分布，可以建立轉(zhuǎn)子有限元單元內(nèi)楊氏模量的變化趨勢，得到考慮軸向溫度分布的單元剛度矩陣及轉(zhuǎn)子有限元模型。

2.1 單元內(nèi)楊氏模量的變化規(guī)律

本文主要考慮轉(zhuǎn)子軸向溫度分布引起楊氏模量的變化對轉(zhuǎn)子振動特性的影響，故假設轉(zhuǎn)子各個截面溫度分布均勻，并且不考慮溫度分布引起的轉(zhuǎn)子的熱變形。

轉(zhuǎn)子有限元單元如圖1所示。假設單元內(nèi)溫度軸向線性分布，兩節(jié)點溫度分別為Ti和Ti+1，則單元內(nèi)任意位置處溫度可以表示為:

需要說明這里可以根據(jù)實際轉(zhuǎn)子軸向溫度分布情況，假設單元內(nèi)溫度具有高階多項式分布形式。

結(jié)合材料楊氏模量隨溫度的變化規(guī)律，將式(3)代入式(1)中，可以得到單元內(nèi)楊氏模量的軸向變化規(guī)律:

其中f(x)與單元溫度分布及軸向坐標有關:

2.2 考慮溫度分布單元剛度矩陣

對于如圖1所示轉(zhuǎn)子有限元單元，假設軸向方向為z軸，則節(jié)點i和節(jié)點i+1各具有四個自由度［22］，即{x，y，－ θx，θy}T。在 XOZ 平面內(nèi)，單元的自由度為{xi，θyi，xi+1，θyi+1}T。單元剛度矩陣為［22］:

其中，E為材料楊氏模量，I為單元截面慣性矩，N為單元型函數(shù)。

圖1 考慮溫度分布的單元模型Fig.1 Element with temperature

結(jié)合單元內(nèi)楊氏模量隨溫度分布的變化規(guī)律，將式(4)代入式(5)，可以得到考慮楊氏模量隨溫度變化的單元剛度矩陣:

根據(jù)f(x)的表達式，可以統(tǒng)一寫為如下形式:

其中，ci與單元節(jié)點溫度和單元內(nèi)的溫度分布相關;n與單元內(nèi)溫度分布形式及楊氏模量隨溫度變化規(guī)律相關。根據(jù)楊氏模量隨溫度二次變化規(guī)律和單元內(nèi)溫度線性分布的假設，這里n最大取2。

將式(7)代入式(6)中，可得考慮溫度分布的轉(zhuǎn)子單元的剛度矩陣:

各元素的通用表達式為:

一般轉(zhuǎn)子橫向剛度相等，在YOZ平面內(nèi)，單元的剛度矩陣與XOZ平面一致。綜合兩個平面內(nèi)的振動，轉(zhuǎn)子單元的橫向振動自由度為{xi，yi，－ θxi，θyi，xi+1，yi+1，－ θxi+1，θyi+1}T，單元的剛度矩陣為:

其中:

由于本文只考慮楊氏模量隨溫度變化對剛度矩陣的影響，不考慮單元的質(zhì)量矩陣、陀螺力矩矩陣受溫度的影響，可采用常規(guī)方法得到［22］。

2.3 考慮楊氏模量隨溫度變化的轉(zhuǎn)子動力學方程

根據(jù)單元的質(zhì)量矩陣、剛度矩陣和陀螺力矩矩陣，可以合并形成轉(zhuǎn)子系統(tǒng)的質(zhì)量矩陣、剛度矩陣和陀螺力矩矩陣，進而可得轉(zhuǎn)子系統(tǒng)的動力學方程:

其中［M］為轉(zhuǎn)子系統(tǒng)的質(zhì)量矩陣，［C］為包括陀螺力矩的轉(zhuǎn)子阻尼矩陣，［Ki］為考慮轉(zhuǎn)子系統(tǒng)溫度分布的剛度矩陣，{Q}為轉(zhuǎn)子節(jié)點載荷向量，{u}為轉(zhuǎn)子節(jié)點自由度向量。

3 結(jié)果分析

為了體現(xiàn)轉(zhuǎn)子軸向溫度分布對轉(zhuǎn)子振動特性的影響和便于驗證模型，下面采用簡單轉(zhuǎn)子模型，分析軸向溫度均勻分布、線性分布和二次分布情況下轉(zhuǎn)子的振動特性。轉(zhuǎn)子的結(jié)構如圖2所示，轉(zhuǎn)子長0.5 m，軸頸0.04 m，兩端剛性支承。轉(zhuǎn)子材料密度為7.85×103kg/m3，常溫下楊氏模量為205 GPa，楊氏模量隨溫度的變化規(guī)律采用式(1)函數(shù)表示。

圖2 模型轉(zhuǎn)子結(jié)構Fig.2 The structure of the rotor model

3.1 軸向溫度均布轉(zhuǎn)子振動特性變化

假設轉(zhuǎn)子溫度均勻分布，可以采用本模型分別計算不同溫度下轉(zhuǎn)子的臨界轉(zhuǎn)速。針對結(jié)構簡單，剛性支承轉(zhuǎn)子，臨界轉(zhuǎn)速可以采用如下解析式計算［23］:

其中:i為臨界轉(zhuǎn)速階數(shù)。

從表達式可以看出，解析式將轉(zhuǎn)子當作簡支梁處理，沒有考慮轉(zhuǎn)子陀螺力矩的影響。為了便于對比分析，在計算溫度均布轉(zhuǎn)子臨界轉(zhuǎn)速時，也不考慮陀螺力矩影響，結(jié)果如表1所示。

表1 模型轉(zhuǎn)子前兩階臨界轉(zhuǎn)速隨溫度變化Tab.1 Variation of model rotor critical speeds with temperature

圖3 模型轉(zhuǎn)子前兩階臨界轉(zhuǎn)速隨溫度的變化趨勢Fig.3 Variation of critical speeds with the rotor temperature

轉(zhuǎn)子臨界轉(zhuǎn)速隨溫度的變化趨勢及與解析解的對比也可以從圖3看出。

從表1和圖3看出，隨著溫度的增加，轉(zhuǎn)子材料楊氏模量減小，轉(zhuǎn)子臨界轉(zhuǎn)速降低，本文模型計算結(jié)果和解析解比較接近。轉(zhuǎn)子在500℃時，轉(zhuǎn)子材料楊氏模量變?yōu)槌刂档?3.39%，本文模型和解析解的前兩階臨界轉(zhuǎn)速均降低8.68%，這驗證了本文模型的正確性。

3.2 軸向溫度線性分布轉(zhuǎn)子振動特性變化

假設轉(zhuǎn)子溫度沿著軸向線性增加，其中低溫端固定為20℃，分別分析模型轉(zhuǎn)子在不同軸向溫差下的臨界轉(zhuǎn)速，結(jié)果見表2。

在軸向溫差作用下，轉(zhuǎn)子臨界轉(zhuǎn)速的變化規(guī)律如圖4所示。從圖中可以看出，隨著軸向溫差的增加，轉(zhuǎn)子的兩階臨界轉(zhuǎn)速是以拋物線規(guī)律降低。在轉(zhuǎn)子軸向溫差為500℃時，轉(zhuǎn)子一階臨界轉(zhuǎn)速降低了3.41%，二階臨界轉(zhuǎn)速降低了3.57%。

表2 模型轉(zhuǎn)子前兩階臨界轉(zhuǎn)速隨軸向溫差變化Tab.2 Variation of model rotor critical speeds with axial temperature difference

圖4 轉(zhuǎn)子前兩階臨界轉(zhuǎn)速隨線性分布軸向溫差變化趨勢Fig.4 Variation of the first two critical speeds with axial temperature difference of linear distribution

在軸向溫差作用下，轉(zhuǎn)子在高溫部位的剛度降低較多，低溫部位的剛度相對較大，這可能會使得轉(zhuǎn)子的振型發(fā)生變化。為此，對比分析了轉(zhuǎn)子在不考慮溫度分布和軸向溫差為500℃的線性分布兩種情況下轉(zhuǎn)子前兩階振型，結(jié)果如圖5所示。從圖中可以看出，在軸向溫差500℃的線性溫度分布下，轉(zhuǎn)子的振型與不考慮溫度情況之間的差別很小，故可以不考慮軸向溫度分布對轉(zhuǎn)子振型的影響。

圖5 轉(zhuǎn)子前兩階振型的變化Fig.5 Variation of the first two modal shape

3.3 軸向溫度二次多項式分布轉(zhuǎn)子振動特性變化

假設轉(zhuǎn)子的溫度在軸向以二次多項式分布，如圖6所示為不同的二次多項式溫度分布，采用有限元方法分別計算各種二次多項式溫度分布下轉(zhuǎn)子的臨界轉(zhuǎn)速。

圖6 不同溫差的轉(zhuǎn)子溫度軸向二次分布Fig.6 Several temperature difference quadratic distribution of rotor axial temperature

表3 不同溫差二次分布下轉(zhuǎn)子臨界轉(zhuǎn)速Tab.3 Rotor critical speeds under several temperature difference quadratic distribution of axial temperature

轉(zhuǎn)子軸向溫度在不同二次多項式分布下的前兩階臨界轉(zhuǎn)速如表3所示，臨界轉(zhuǎn)速的變化趨勢如圖7所示。轉(zhuǎn)子在軸向溫度二次分布下，臨界轉(zhuǎn)速變化比線性分布大。在最大溫差為500℃時，轉(zhuǎn)子的一階臨界轉(zhuǎn)速降低了7.11%，二階臨界轉(zhuǎn)速降低了5.46%?？梢娹D(zhuǎn)子軸向溫度分布不同，對轉(zhuǎn)子臨界轉(zhuǎn)速的影響差別較大。

圖7 不同二次多項式溫度分布下轉(zhuǎn)子臨界轉(zhuǎn)速變化趨勢Fig.7 Variation of first two critical speeds under different axial temperature quadratic distribution

從上面三種不同形式軸向溫度分布情況下轉(zhuǎn)子臨界轉(zhuǎn)速的變化可以看出，本文的模型能夠充分考慮轉(zhuǎn)子軸向溫度分布的不同，體現(xiàn)轉(zhuǎn)子軸向溫度分布對轉(zhuǎn)子材料楊氏模量的影響，為轉(zhuǎn)子熱彈耦合動力學分析奠定了基礎，可以提高轉(zhuǎn)子在實際工作中臨界轉(zhuǎn)速的計算精度。

需要說明的是，本文仿真中所采用的轉(zhuǎn)子軸向溫度分布都是理想化的，實際機組中轉(zhuǎn)子軸向溫度分布是十分復雜的。例如汽輪機在工作過程中，轉(zhuǎn)子的軸向溫度變化規(guī)律在各級之間是不同的，同時軸向溫度分布也受到汽輪機啟動過程的影響。為此，在應用本文模型計算實際轉(zhuǎn)子臨界轉(zhuǎn)速時，需要首先確定轉(zhuǎn)子的溫度分布，得到轉(zhuǎn)子有限元模型各節(jié)點溫度。此外，本模型基于轉(zhuǎn)子材料楊氏模量隨溫度的變化規(guī)律，在計算中需針對不同材料，采用不同的變化規(guī)律。

4 結(jié)論

結(jié)合轉(zhuǎn)子材料楊氏模量隨溫度的變化規(guī)律，推導了考慮轉(zhuǎn)子溫度分布的單元剛度矩陣，建立了考慮轉(zhuǎn)子軸向溫度分布的有限元模型。仿真分析了模型轉(zhuǎn)子在軸向溫度均布、線性分布和二次多項式分布下，轉(zhuǎn)子臨界轉(zhuǎn)速的變化趨勢，結(jié)果表明模型能夠較好考慮轉(zhuǎn)子軸向溫度分布對轉(zhuǎn)子臨界轉(zhuǎn)速的影響，提高高溫工作轉(zhuǎn)子臨界轉(zhuǎn)速計算精度。

［1］陸 山，趙 明，任平珍，等.某型發(fā)動機轉(zhuǎn)子熱彎曲變形及其影響數(shù)值分析［J］.航空動力學報，1997，12(3):243－246.LU Shan，ZHAO Ming，REN Ping-zhen，et al.A numerical analysis of thermal bending deformation and its influence on rotor［J］.Journal of Aerospace Power，1997，12(3):243 －246.

［2］莊景菁.轉(zhuǎn)子熱彎曲引起汽輪機振動的分析和處理［J］.華東電力，2002，30(7):61－64.ZHUANG Jing-qing.Analysisand treatmentofturbine vibration caused by thermal bending of rotor［J］.East China Electric Power，2002，30(7):61 －64.

［3］吳子斌.汽輪發(fā)電機轉(zhuǎn)子突發(fā)性振動原因［J］.東北電力技術，2003，24(1):28－29.WU Zi-bin，The cause of burst vibration of the rotor of turbogenerator［J］.Northeastern Electric Power Technology，2003，24(1):28 －29.

［4］戚先萍，胡璧剛，許都純，等.轉(zhuǎn)子熱彎曲瞬態(tài)溫度場及熱彎曲變形試驗與分析［J］.燃氣渦輪試驗與研究，1997，10(2):32 －36.QI Xian-ping，HU Bi-gang，XU Du-chun，et al.Experiment and analysis of the transient temperature distribution and the thermal bending of rotor［J］.Gas Turbine Experiment and Research，1997，10(2):32 －36.

［5］Goldman P，Muszynska A，Bently D E.Thermal bending of the rotor due to rotor-to-stator rub［J］.International Journal，2000，6(2):91 －100.

［6］黃葆華，楊建剛.摩擦熱彎曲對轉(zhuǎn)子穩(wěn)定性的影響［J］.振動工程學報，2001，14(1):23－28.HUANG Bao-hua，YANG Jian-gang.Influence of thermal bow of the stability of rubbing bearing-rotor system［J］.Journal of Vibration Engineering，2001，14(1):23－28.

［7］田永偉，楊建剛.基于特征值分析的摩擦轉(zhuǎn)子熱穩(wěn)定性研究［J］.振動與沖擊，2010，29(7):66－68.TIAN Yong-wei，YANG Jian-gang.Thermal stability of a rubbing rotor based on eigenvalue analysis［J］.Journal of Vibration and Shock，2010，29(7):66 －68.

［8］賀 威，袁惠群，朱向哲.瞬態(tài)傳遞矩陣法分析熱彎轉(zhuǎn)子系統(tǒng)的熱碰摩故障［J］.東北大學學報(自然科學版)，2009(6):861－864.HE Wei， YUAN Hui-qun， ZHU Xiang-zhe. Transient transfer matrix method for thermal rub-impact faults in hot bending rotor system［J］.Journal of Northeastern University(Natural Science)，2009(6):861－864.

［9］Keogh P S，Morton P G.Journal bearing differential heating evaluation with influence on rotor dynamic behaviour［J］.Proceedings oftheRoyalSocietyofLondon.SeriesA:Mathematical and Physical Sciences，1993，441(1913):527－548.

［10］Keogh P S，Morton P G.The dynamic nature of rotor thermal bending due to unsteady lubricant shearing within a bearing［J］.Proceedings of the Royal Society of London.Series A:Mathematical and Physical Sciences，1994，445(1924):273.

［11］任平珍，柴衛(wèi)東，胡壁剛，等.航空發(fā)動機轉(zhuǎn)子熱彎曲穩(wěn)態(tài)響應計算方法研究［J］.燃氣渦輪試驗與研究，1996(3):27－32.REN Ping-zhen， CHAI Wei-dong， HU Bi-gang， et al.Research on the computation method of the steady state response of the thermal bending aeroengine rotor［J］.Gas Turbine Experiment and Research，1996(3):27－32.

［12］任平珍，陸 山，趙 明.轉(zhuǎn)子熱彎曲變形及其影響的數(shù)值分析方法［J］.機械科學與技術，1997，16(2):279－282.REN Ping-zhen，LU Shan，ZHAO Ming.Numerical analysis method of thermal deflection and its affection on vibration response of rotor［J］.Mechanical Science and Technology，1997，16(2):279 －282.

［13］陸 山，杜生廣，趙 明，等.降低航空發(fā)動機轉(zhuǎn)子熱彎曲及其影響方法研究［J］.航空發(fā)動機，1997(3):23－27.LU Shan，DU Sheng-guang，ZHAO Ming，et al.Research on reducing the thermal bending and its influence of aeroengine rotor［J］.Aeroengine，1997(3):23 －27.

［14］張連祥.航空發(fā)動機轉(zhuǎn)子熱彎曲引發(fā)的典型故障分析［J］.振動與沖擊(增刊)，2008，27(S):7－9.ZHANG Lian-xiang.Analysis of typical aeroengine vibration fault related to thermal bowed rotor［J］.Journal of Vibration and Shock，2008，27(S):7－9.

［15］張麗娟，楊新華，馬呈霞，等.轉(zhuǎn)子熱彎曲引起發(fā)電機組振動的診斷與處理［J］.大電機技術，2009(4):13－16.ZHANG Li-juan，YANG Xin-hua，MA Cheng-xia，et al.Remedy on vibration of turbogenerator due to thermal-caused bending of rotor［J］.Large Electric Machine and Hydraulic Turbine，2009(4):13－16.

［16］Kimball Jr A L，Lovell D E.Variation of young's modulus with temperature from vibration measurements［J］.Physical Review，1925，26(1):121 －124.

［17］潘留仙，焦善慶，杜小勇.高溫下常用合金材料線脹系數(shù)，楊氏模量與溫度的關系［J］.湖南師范大學自然科學學報，2000，23(2):47－51.PAN Liu-xian，JIAO Shan-qing，DU Xiao-yong.The relation between linear expansion coefficient and yang's modulus of alloy materials ＆ temperature at high temperature［J］.Journal of Natural Science of Hunan Normal University，2000，23(2):47 －51.

［18］晏水平，黃樹紅.汽輪發(fā)電機組轉(zhuǎn)子溫度分布對其扭轉(zhuǎn)振動的影響［J］.中國電機工程學報，2000，20(11):10－13.YAN Shui-ping，HUANG Shu-hong.Impact of temperature distribution on the torsional vibration of turbogenerator rotors［J］.Chinese Society for Electrical Engineering，2000，20(11):10－13.

［19］朱向哲，袁惠群，賀 威.穩(wěn)態(tài)熱度場對轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速的影響［J］.振動與沖擊，2007，26(12):113－116.ZHU Xiang-zhe，YUAN Hui-qun，HE Wei.Effect of steady thermal field on critical speeds of a rotor system［J］.Journal of Vibration and Shock，2007，26(12):113－116.

［20］朱向哲，賀 威，袁惠群.穩(wěn)態(tài)溫度場對轉(zhuǎn)子系統(tǒng)振動特性的影響［J］.東北大學學報:自然科學版，2008(1):113－116.ZHU Xiang-zhe，HE Wei，YUAN Hui-qun.Effects of steady temperature field on vibrational characteristics of a rotor system［J］.Journal of Northeastern University(Natrual Science)，2008(1):113－116.

［21］劉少權，張艷春，杜兆剛，等.溫度場對燃氣輪機拉桿轉(zhuǎn)子臨界轉(zhuǎn)速的影響［J］.燃氣輪機技術，2011，24(2):20－23.LIU Shao-quan，ZHANG Yan-chun，DU Zhao-gang，et al.Prediction of the influence of temperature field on the critical speeds of a rod-fastened rotor［J］.Gas Turbine Technology，2011，24(2):20 －23.

［22］鐘一諤，何衍宗，王 正，等.轉(zhuǎn)子動力學［M］.北京:清華大學出版社，1987:176－185.

［23］GentaG. Consistent matricesin rotordynamic［J］.Meccanica，1985，20(3):235 －248.