長葉片相鄰圍帶或凸肩接觸轉(zhuǎn)速計(jì)算方法

陳德祥， 徐自力， 曹守洪， 范小平， 吳其林

（1.西安交通大學(xué) 航天學(xué)院，機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，西安710049；2.東方汽輪機(jī)有限公司，德陽618000）

成圈葉片只有在“三重點(diǎn)”條件下才會(huì)產(chǎn)生共 振，在固有頻率和氣流激振力頻率相等而不滿足“三重點(diǎn)”的振型條件下，葉片并不會(huì)產(chǎn)生共振［1－2］，因此現(xiàn)代汽輪機(jī)動(dòng)葉中的長葉片多采用圍帶或圍帶加凸肩的結(jié)構(gòu)（圖1），使葉片在工作時(shí)處于自鎖的成圈狀態(tài).葉片圍帶有2種成圈方式：一是采用預(yù)扭設(shè)計(jì)［3］，這種方式采用的圍帶幾何形狀能使葉片裝配完成后相鄰的葉片圍帶相互接觸并擠壓；二是采用間隙設(shè)計(jì)［4－6］，這種方式初始狀態(tài)下相鄰圍帶之間存在一定的間隙，隨著轉(zhuǎn)速的加快，葉片發(fā)生扭轉(zhuǎn)恢復(fù)，最終使圍帶相互接觸并擠壓（圖2）.第二種方式下，葉片在裝配以及更換時(shí)比預(yù)扭葉片簡單.

圖1 低壓末級(jí)葉片F(xiàn)ig.1 Low－pressure last－stage blade

圖2 葉片扭轉(zhuǎn)恢復(fù)引起的圍帶接觸Fig.2 Shroud contact caused by blade untwisting

理論上，圍帶或凸肩之間的接觸法向力或接觸應(yīng)力定量地反映了葉片的壓緊程度和成圈性，是葉片設(shè)計(jì)中的重要參數(shù)之一，它的取值一方面與其所遵從的設(shè)計(jì)理念有關(guān)，是側(cè)重于成圈整體性還是側(cè)重于摩擦阻尼減振性，另一方面也與葉片的長度、工作條件有關(guān).但是接觸法向力或接觸應(yīng)力的大小難以在試驗(yàn)中直接測(cè)量.一般可通過測(cè)量安裝完成后的圍帶預(yù)扭角度來確定預(yù)扭設(shè)計(jì)葉片圍帶的壓緊程度，但對(duì)間隙設(shè)計(jì)的葉片只有在轉(zhuǎn)動(dòng)狀態(tài)下才形成成圈結(jié)構(gòu)，葉片之間的扭轉(zhuǎn)角度無法通過試驗(yàn)直接測(cè)量.圍帶開始接觸時(shí)的轉(zhuǎn)速宏觀上反映了間隙設(shè)計(jì)葉片相鄰圍帶之間的壓緊程度，進(jìn)而反映了葉片的自鎖成圈性，葉型一定時(shí)接觸轉(zhuǎn)速越小，則工作狀態(tài)下接觸面上相互擠壓越緊.重要的是，接觸轉(zhuǎn)速可以通過試驗(yàn)進(jìn)行測(cè)量，表現(xiàn)在坎貝爾圖上接觸轉(zhuǎn)速點(diǎn)的頻率發(fā)生跳躍.接觸轉(zhuǎn)速的測(cè)量值與設(shè)計(jì)值的符合程度反映了該級(jí)葉片的整體制造質(zhì)量，因此計(jì)算圍帶開始接觸時(shí)的轉(zhuǎn)速對(duì)設(shè)計(jì)成圈葉片有著重要意義.

筆者用接觸法向力來表征相鄰圍帶或凸肩之間的接觸關(guān)系，并通過迭代計(jì)算成圈葉片的接觸轉(zhuǎn)速.在給定轉(zhuǎn)速下，影響相鄰圍帶或凸肩上接觸法向力的因素復(fù)雜，難以給出其表達(dá)式，而帶接觸的三維有限元分析能考慮各種復(fù)雜因素［7－9］，故通過三維接觸有限元分析獲得相鄰圍帶或凸肩上的接觸法向力.根據(jù)整圈葉片幾何上具有循環(huán)對(duì)稱的特點(diǎn)，利用循環(huán)對(duì)稱算法來降階［10－11］，即通過計(jì)算1個(gè)扇區(qū)來獲得整圈的結(jié)果.由于通過三維接觸有限元分析只得到接觸法向力的函數(shù)值，沒有導(dǎo)數(shù)值，并且接觸法向力在未接觸時(shí)一直為0，所以牛頓法或割線法迭代不適用，故采用二分法進(jìn)行接觸轉(zhuǎn)速的計(jì)算.每迭代一次需要進(jìn)行一次有接觸的三維有限元分析，取上一次有限元分析的結(jié)果作為本次有限元分析初始值來提高收斂速度.

1 圍帶、凸肩接觸轉(zhuǎn)速的計(jì)算方法

1.1 相鄰圍帶、凸肩之間關(guān)系的量化表征

相鄰圍帶之間關(guān)系包括相互分離和相互接觸2種狀態(tài)，分離狀態(tài)下可用相鄰圍帶之間的間隙d來表征，接觸狀態(tài)下可用圍帶之間的法向相互作用力F來表征（圖2）.葉片的扭轉(zhuǎn)恢復(fù)主要是由離心力引起的，并且隨著轉(zhuǎn)速的加快而增加，因此d和F都是轉(zhuǎn)速的函數(shù)，即d＝d（ω），F(xiàn)＝F（ω），他們隨轉(zhuǎn)速變化的趨勢(shì)如圖3所示.當(dāng)ω＜ωc時(shí)，相鄰圍帶間隙隨著轉(zhuǎn)速加快單調(diào)減小，法向作用力保持為0；當(dāng)ω＞ωc時(shí)，法向作用力隨轉(zhuǎn)速加快單調(diào)增大，而間隙值保持為0，所以圍帶接觸轉(zhuǎn)速為：

圖3 圍帶之間接觸關(guān)系表征量隨轉(zhuǎn)速的變化Fig.3 Variation of characteristic parameters of shroud with rotational speed

根據(jù)圍帶之間接觸狀態(tài)的量化表征，確定圍帶接觸轉(zhuǎn)速最終歸結(jié)為尋找d（ω）＝0的最小根，或者尋找F（ω）＝0的最大根.相鄰?fù)辜缰g關(guān)系與圍帶的情況完全相同.

1.2 計(jì)算方法

影響相鄰圍帶、凸肩間隙和接觸法向力大小的因素復(fù)雜.首先，離心力與轉(zhuǎn)速為平方關(guān)系，間隙和接觸法向力隨轉(zhuǎn)速非線性變化；其次，圍帶、凸肩以及葉根處存在接觸，接觸問題本身具有非線性，并且接觸狀態(tài)受離心力的影響；再次，末級(jí)葉片一般柔度大，這導(dǎo)致幾何非線性問題；最后，如果葉片上同時(shí)有圍帶和凸肩結(jié)構(gòu)，它們相互之間也有影響.因此，圖3為間隙和接觸法向力隨轉(zhuǎn)速的變化趨勢(shì)，但無法表達(dá)成轉(zhuǎn)速的函數(shù).

筆者采用考慮接觸的三維有限元方法確定某一給定轉(zhuǎn)速下間隙或接觸法向力的值.但通過有限元計(jì)算只能得到間隙或接觸法向力的函數(shù)值，一般沒有導(dǎo)數(shù)信息，因此迭代中需要用到導(dǎo)數(shù)的梯度類算法，牛頓法不能用于該問題的求解.另外，從圖3可以看出，不論采用間隙還是接觸法向力進(jìn)行迭代，在全轉(zhuǎn)速范圍內(nèi)方程d（ω）＝0或F（ω）＝0的根不是唯一的，應(yīng)用非梯度類算法中的割線法能計(jì)算出其中的一個(gè)根，但不能得到d（ω）＝0的最小根或F（ω）＝0的最大根，因此割線法也不適用于該問題的求解.二分法不要求導(dǎo)數(shù)信息，也沒有割線法的問題，適用于接觸轉(zhuǎn)速的迭代計(jì)算.

雖然理論上間隙值和接觸法向力都可以用于迭代計(jì)算，最終的收斂結(jié)果是相同的.但是進(jìn)入接觸之后，圍帶、凸肩的間隙都為0，間隙值不能反映成圈性的變化，而接觸法向力則可以直接反映壓緊程度，比間隙值更有工程價(jià)值，因此對(duì)接觸法向力進(jìn)行迭代計(jì)算，其算法流程圖見圖4.

假定ωc所在的轉(zhuǎn)速區(qū)間是（ωa，ωb），采用二分法計(jì)算時(shí)，每迭代一次，ωc所在轉(zhuǎn)速區(qū)間長度減小一半，直到轉(zhuǎn)速區(qū)間的長度滿足給定的收斂條件.收斂條件為：

式中：ωn為工作轉(zhuǎn)速；ωa、ωb分別為ωc所在轉(zhuǎn)速區(qū)間的上、下界；ε為給定的正常數(shù).

2 算例

2.1 計(jì)算對(duì)象及有限元網(wǎng)格

圖4 采用二分法計(jì)算接觸轉(zhuǎn)速流程圖Fig.4 Bisection iteration of contact speed calculation

采用以上方法計(jì)算某200MW機(jī)組末級(jí)動(dòng)葉片圍帶和凸肩的接觸轉(zhuǎn)速.葉片長度為816mm，同時(shí)具有圍帶和凸肩結(jié)構(gòu)（圖1），圍帶之間間隙為0.5 mm，凸肩之間間隙為0.2mm.利用循環(huán)對(duì)稱特點(diǎn)，取整級(jí)葉片和輪盤的1個(gè)扇區(qū)（圖5（a））進(jìn)行計(jì)算，計(jì)算中考慮相鄰?fù)辜缰g的接觸（圖5（b））、相鄰圍帶之間的接觸（圖5（c））以及葉根與輪盤之間的接觸（圖5（d））.因?yàn)樾枰紤]圍帶和凸肩的接觸，所以葉片上的循環(huán)對(duì)稱面選在葉片內(nèi)部，先對(duì)整個(gè)葉片進(jìn)行網(wǎng)格剖分，然后將圍帶和凸肩上的一部分網(wǎng)格轉(zhuǎn)動(dòng)1個(gè)扇區(qū)，轉(zhuǎn)動(dòng)之后自然形成逐點(diǎn)匹配的循環(huán)對(duì)稱節(jié)點(diǎn).采用8節(jié)點(diǎn)6面體單元進(jìn)行空間離散，單元總數(shù)為127 929個(gè)，其中葉片劃分了92 929個(gè)單元，輪盤劃分了35 000個(gè)單元；節(jié)點(diǎn)總數(shù)為148 360個(gè)，其中葉片節(jié)點(diǎn)數(shù)為106 150個(gè)，輪盤節(jié)點(diǎn)數(shù)為42 210個(gè).

2.2 計(jì)算結(jié)果

取收斂條件為ε＝0.01，圍帶、凸肩接觸轉(zhuǎn)速的迭代收斂過程見表1和表2.從表中可以看出，在迭代過程中，區(qū)間（ωa，ωb）的下限和上限分別從小于接觸轉(zhuǎn)速和大于接觸轉(zhuǎn)速的兩側(cè)向接觸轉(zhuǎn)速收斂.當(dāng)接觸力法向?yàn)?時(shí)，接觸轉(zhuǎn)速所在區(qū)間的下限逐漸增大；當(dāng)接觸力法向大于0時(shí)，接觸轉(zhuǎn)速所在區(qū)間的下限、上限均逐漸減小，對(duì)應(yīng)的接觸法向力也逐漸減小.將收斂時(shí)轉(zhuǎn)速區(qū)間的平均值作為接觸轉(zhuǎn)速，從表中可得到圍帶的接觸轉(zhuǎn)速為738r／min，凸肩的接觸轉(zhuǎn)速為1 184r／min.雖然該葉片圍帶間隙大于凸肩間隙，但圍帶比凸肩先進(jìn)入接觸狀態(tài).該葉片在工作轉(zhuǎn)速下圍帶和凸肩接觸法向力最大，對(duì)應(yīng)在圍帶上的平均接觸應(yīng)力為14.5MPa，凸肩上的平均接觸應(yīng)力為18.5MPa.

圖5 某200MW機(jī)組低壓末級(jí)葉片、輪盤模型Fig.5 Meshing of a low－pressure last－stage blade and disk for a 200MW unit

表1 計(jì)算圍帶接觸轉(zhuǎn)速的二分法迭代過程Tab.1 Bisection iteration process of shroud contact speed calculation

根據(jù)圍帶、凸肩每一次迭代的計(jì)算結(jié)果，給出圍帶、凸肩上的接觸法向力隨轉(zhuǎn)速（300～1 500 r／min）的變化情況，結(jié)果見圖6.由于圍帶接觸比凸肩接觸出現(xiàn)得早，從圍帶進(jìn)入接觸到凸肩進(jìn)入接觸的這一過程中，離心力產(chǎn)生的扭轉(zhuǎn)恢復(fù)力矩與圍帶之間法向力產(chǎn)生的力矩平衡，圍帶上的法向力單調(diào)增大.當(dāng)凸肩進(jìn)入接觸后，一部分扭轉(zhuǎn)恢復(fù)力矩與凸肩法向力產(chǎn)生的力矩平衡，因此圍帶上的法向力隨轉(zhuǎn)速加快而增大的速度在此時(shí)有所減小.

表2 計(jì)算凸肩接觸轉(zhuǎn)速的二分法迭代過程Tab.2 Bisection iteration process of lacing contact speed calculation

圖6 某低壓末級(jí)葉片圍帶、凸肩上接觸法向力隨轉(zhuǎn)速的變化Fig.6 Variation of normal contact force with rotational speed of shroud and lacing on a low－pressure last－stage blade

3 計(jì)算速度和收斂條件分析

以上方法的最大計(jì)算成本在于三維接觸有限元計(jì)算，對(duì)所研究的葉片共進(jìn)行了8次有限元計(jì)算，包括迭代開始之前工作轉(zhuǎn)速下的一次有限元計(jì)算.接觸問題具有非線性特點(diǎn)，其求解過程通常需要進(jìn)行迭代，非線性問題的迭代收斂速度與初始值選擇有關(guān)，初始值越接近真實(shí)解其收斂速度越快.由于上一次迭代已經(jīng)在某一轉(zhuǎn)速下達(dá)到平衡狀態(tài)，其計(jì)算結(jié)果比任意給定的一組值更接近本次迭代轉(zhuǎn)速下的真實(shí)解，因此筆者將上一次有限元計(jì)算結(jié)果作為本次接觸有限元迭代的初始值.采用8核CPU，16G內(nèi)存計(jì)算機(jī)進(jìn)行計(jì)算，圍帶或凸肩接觸轉(zhuǎn)速的迭代計(jì)算過程需7h，其中第一次接觸有限元分析需2h，之后每次接觸有限元分析約43min，計(jì)算速度顯著加快.

采用二分法迭代k次后，解所在區(qū)間的長度是初始區(qū)間長度的1／2k.對(duì)工作轉(zhuǎn)速為3 000r／min的機(jī)組，取收斂條件ε＝0.01，需要迭代7次，收斂時(shí)接觸轉(zhuǎn)速所在區(qū)間的長度為23.4r／min，取ωc＝（ωa＋ωb）／2，則接觸轉(zhuǎn)速的計(jì)算誤差小于11.7r／min.如果將收斂條件提高一個(gè)量級(jí)，取ε＝0.001，這時(shí)需要迭代10次，接觸轉(zhuǎn)速的計(jì)算誤差小于1.5 r／min.由于ε＝0.01給出的計(jì)算精度基本滿足描述成圈性對(duì)精度的要求，且當(dāng)法向力很小時(shí)，三維接觸有限元計(jì)算的誤差成為主要因素，增加迭代次數(shù)并不能提高計(jì)算精度，反而增加計(jì)算成本，因此，并非迭代次數(shù)越多越好，ε＜0.001是不必要的.筆者建議迭代收斂條件在0.01～0.001之間，即整個(gè)計(jì)算過程迭代7～10次.

4 結(jié) 論

汽輪機(jī)末級(jí)長葉片圍帶、凸肩的接觸轉(zhuǎn)速反映了整級(jí)葉片的成圈性，對(duì)葉片成圈性的定量評(píng)價(jià)具有重要意義.用接觸法向力來表征相鄰圍帶、凸肩的接觸關(guān)系，采用二分法計(jì)算圍帶和凸肩的接觸轉(zhuǎn)速.迭代中取上一次迭代的有限元分析結(jié)果作為本次迭代中接觸有限元分析的初始值，提高了迭代的收斂速度.計(jì)算接觸法向力時(shí)采用的三維有限元法考慮了離心力、接觸非線性、幾何非線性、圍帶和凸肩相互作用等多種因素的影響，具有較高的計(jì)算精度.應(yīng)用上述方法對(duì)某200MW汽輪機(jī)組末級(jí)葉片圍帶、凸肩的接觸轉(zhuǎn)速進(jìn)行了計(jì)算，經(jīng)過7次迭代，接觸轉(zhuǎn)速的計(jì)算誤差小于11.7r／min.本文方法可用于長葉片設(shè)計(jì)中相鄰圍帶或凸肩接觸轉(zhuǎn)速的計(jì)算.

［1］HUANG W H.Free and forced vibration of closely coupled turbomachinery blades［J］.AIAA Journal，1981，19（7）：918－924.

［2］徐芬，程凱，彭澤瑛.圍帶結(jié)構(gòu)對(duì)整圈自鎖葉片振動(dòng)特性的影響［J］.動(dòng)力工程學(xué)報(bào)，2010，30（5）：347－351.XU Fen，CHENG Kai，PENG Zeying.Influence of shroud configuration on vibration characteristics of the full circle self－lock blade［J］.Journal of Chinese Society of Power Engineering，2010，30（5）：347－351.

［3］楊鑫，馬艷紅，洪杰.基于接觸狀態(tài)的葉冠預(yù)扭設(shè)計(jì)［J］.航空發(fā)動(dòng)機(jī)，2008，34（4）：11－15.YANG Xin，MA Yanhong，HONG Jie.Pretwisted design of shroud at contact state［J］.Aeroengine，2008，34（4）：11－15.

［4］KANEKO Y，MORI K，OHYAMA H.Development and verification of 3 000rpm 48inch integral shroud blade for steam turbine［J］.JSME International Journal Series B，2006，49（2）：205－211.

［5］KANEKO Y，OHYAMA H.Analysis and measurement of damping characteristics of integral shroud blade for steam turbine［J］.Journal of System Design and Dynamics，2008，2（1）：311－322.

［6］FUKUDA H，OHYAMA H，MORI TMK，et al.Development of 3，600－rpm 50－inch／3，000－rpm 60－inch ultra－long exhaust end blades［J］.Mitsubishi Heavy Industries Technical Review，2009，46（2）：18－25.

［7］PAPANIKOS P，MEGUID S，STJEPANOVIC Z.Three－dimensional nonlinear finite element analysis of dovetail joints in aeroengine discs［J］.Finite Elements in Analysis and Design，1998，29（3／4）：173－186.

［8］MEGUID S，KANTH P，CZEKANSKI A.Finite element analysis of fir－tree region in turbine discs［J］.Finite Elements in Analysis and Design，2000，35（4）：305－317.

［9］BEISHEIM J，SINCLAIR G.On the three－dimensional finite element analysis of dovetail attachments［J］.Journal of Turbomachinery，2003，125（2）：372－378.

［10］張錦，劉曉平.葉輪機(jī)振動(dòng)模態(tài)分析理論及數(shù)值方法［M］.北京：國防工業(yè)出版社，2001.

［11］董廣明，陳進(jìn)，孟慶集.循環(huán)對(duì)稱算法及其在汽輪機(jī)葉盤系統(tǒng)振動(dòng)分析中的應(yīng)用［J］.振動(dòng)工程學(xué)報(bào)，2004，17（z1）：126－129.DONG Guangming，CHEN Jin，MENG Qingji.Cyclic symmetry algorithm with its application in vibration analysis of turbine bladed discs［J］.Journal of Vibration Engineering，2004，17（z1）：126－129.