張明根，郝小龍，趙經(jīng)明，郭軍剛，解俊良

(北京精密機(jī)電控制設(shè)備研究所，航天伺服驅(qū)動與傳動技術(shù)實驗室，北京 100076)

0 引 言

渦輪泵作為航天伺服系統(tǒng)核心元件，具有體積小、轉(zhuǎn)速高、比功率大的特點。其基本工作原理為燃?xì)獍l(fā)生器裝置點火后，噴射出高溫高壓燃?xì)怛?qū)動渦輪促使軸系高速旋轉(zhuǎn)，泵葉輪在泵腔中產(chǎn)生高壓油液推動活塞做功。目前渦輪泵主要應(yīng)用于液體火箭發(fā)動機(jī)、航天伺服系統(tǒng)、魚雷發(fā)射系統(tǒng)、航天飛機(jī)主發(fā)動機(jī)等領(lǐng)域。相關(guān)文獻(xiàn)主要集中于液體火箭發(fā)動機(jī)用渦輪泵[1-6]，包括支承剛度變化對轉(zhuǎn)子穩(wěn)定性影響、不平衡質(zhì)量影響下的瞬態(tài)響應(yīng)、軸承位置對臨界轉(zhuǎn)速影響等。魚雷發(fā)射系統(tǒng)及航天飛機(jī)用渦輪泵研究較少[7-8]。在伺服系統(tǒng)用渦輪泵方面，李振將等[9]分別使用等效剛度法、經(jīng)典理論法及傳遞矩陣法計算了渦輪泵轉(zhuǎn)子的一二階臨界轉(zhuǎn)速，并研究了滾動軸承剛度等因素對臨界轉(zhuǎn)速的影響；郭軍剛等[10]使用傳遞矩陣法對箔片軸承支承的超高速渦輪轉(zhuǎn)子系統(tǒng)進(jìn)行了分析，并針對試驗數(shù)據(jù)進(jìn)行了研究；筆者也對該類型渦輪泵進(jìn)行了熱力學(xué)研究，并分析了受溫度影響的轉(zhuǎn)子動力學(xué)響應(yīng)[11]?？偟膩碚f，針對該類型渦輪泵動力學(xué)的研究較少，現(xiàn)有的一些研究所使用方法單一，目前僅局限于傳遞矩陣法，難以進(jìn)行更深一步的研究。目前液發(fā)渦輪泵的動力學(xué)已經(jīng)有了較全面的研究成果，但由于伺服用渦輪泵的結(jié)構(gòu)形式及工作特點均與其存在較大的差異，因此相關(guān)研究僅具備一定的參考價值。為了全面掌握伺服系統(tǒng)用渦輪泵的動力學(xué)特性，提高產(chǎn)品在高轉(zhuǎn)速條件下的穩(wěn)定性，需要針對該類型渦輪泵開展大量專門研究。此外不同于以往的研究方法，筆者采用以有限元原理為基礎(chǔ)的ANSYS三維仿真軟件，通過APDL編制復(fù)雜軸系截面文件，并在主程序中調(diào)用的方法，實現(xiàn)了轉(zhuǎn)子動力學(xué)建模仿真過程參數(shù)化、分析通用化。以此為基礎(chǔ)進(jìn)行了渦輪泵的轉(zhuǎn)子臨界轉(zhuǎn)速計算，并針對軸系在不平衡質(zhì)量影響下的動力學(xué)特性進(jìn)行了研究。得到了渦輪泵軸系受質(zhì)量力影響下的振動響應(yīng)以及滾動軸承支承反力變化規(guī)律。

1 有限元法的基本原理

1.1 轉(zhuǎn)子系統(tǒng)的運動方程

根據(jù)有限元法的分析方法[12]，轉(zhuǎn)子系統(tǒng)的單元主要分為軸段、剛性圓盤和軸承座等，由于各自的受力、變形特點不同，因此運動方程存在差異。

通過建立軸段結(jié)點的廣義坐標(biāo)，假定形函數(shù)矩陣，推導(dǎo)出單元的移動慣性矩陣[MzT]、移動慣性矩陣[MzR]、回轉(zhuǎn)矩陣[Gz]=Ω[Jz]和剛度矩陣[Kz]，由Lagrange方程得到軸段單元的運動方程為：

(1)

式中：u1z、u2z為廣義坐標(biāo)；Q1z、Q2z為廣義坐標(biāo)下的力。

已知剛性圓盤的質(zhì)量、赤道轉(zhuǎn)動慣量、極轉(zhuǎn)動慣量，推導(dǎo)出動能表達(dá)式，并得到動能的微量方程，變換得到圓盤單元的運動方程：

(2)

將滑動軸承簡化為質(zhì)量-彈簧-阻尼器模型，動力特性系數(shù)矩陣為：

(3)

(4)

軸承中心坐標(biāo)為xb、yb對應(yīng)軸頸中心結(jié)點的編號是z(j)，則軸頸中心的坐標(biāo)為xz(j)和yz(j)，則運動方程為：

(5)

具有n個結(jié)點，n-1個軸段連接而成的轉(zhuǎn)子系統(tǒng)的位移向量為：

(6)

綜合圓盤、軸段及軸承的運動方程，可以得到系統(tǒng)的運動方程：

(7)

式中：[M1]為整體質(zhì)量矩陣；Ω[J1]為回轉(zhuǎn)矩陣；[K1]為剛度矩陣。

1.2 臨界轉(zhuǎn)速及不平衡響應(yīng)分析

根據(jù)上面建立的系統(tǒng)運動方程，求解其齊次解，相應(yīng)的能夠得到自轉(zhuǎn)角速度Ω下的轉(zhuǎn)子渦動頻率ω，也就得到了當(dāng)Ω=ω時，轉(zhuǎn)子的臨界轉(zhuǎn)速及相應(yīng)振型。

對于滾動軸承支承的軸系來說，kxx=kyy,kxy=kyx=0，忽略阻尼及軸承座等效質(zhì)量，轉(zhuǎn)子系統(tǒng)運動微分方程的齊次式為：

(8)

給定自轉(zhuǎn)角速度后，得到頻率方程：

|-M1ω2+J1Ωω+K1|=0

(9)

式中：令Ω=ω，則可以計算出系統(tǒng)的各階臨界轉(zhuǎn)速。

轉(zhuǎn)子在不平衡力、力矩的激勵下所產(chǎn)生的振動為不平衡響應(yīng)。對于各向同性軸承如滾動軸承，不計阻尼和軸承座振動影響，系統(tǒng)的不平衡響應(yīng)方程為：

(10)

其中等式右邊為不平衡激勵的廣義力。令：

{z}={U1}+{U2}

(11)

不平衡響應(yīng)方程可以寫為：

(12)

求解上式得到不平衡響應(yīng)的特解：

{z}={A}eiΩt

(13)

其中，{A}=Ω2[-M1Ω2+J1Ω2+K1]-1{Q}

2 渦輪泵軸系動力學(xué)分析

2.1 APDL參數(shù)化仿真計算流程及方法

航天伺服系統(tǒng)用超高速渦輪泵軸系基本結(jié)構(gòu)為：兩端懸臂分別為泵葉輪、渦輪盤，中間為兩個尺寸為滾動軸承支承。由于輪盤結(jié)構(gòu)復(fù)雜，采用基于極轉(zhuǎn)動慣量的方法對軸系進(jìn)行模塊化建模，建模后的簡化軸系結(jié)構(gòu)如圖1所示。建模后的軸系結(jié)構(gòu)均由規(guī)則盤軸單元組成，有利于APDL參數(shù)化仿真。

圖1 渦輪泵軸系模型

按照以下方法進(jìn)行APDL轉(zhuǎn)子動力學(xué)建模仿真：

(1) 根據(jù)轉(zhuǎn)子的結(jié)構(gòu)特點，將轉(zhuǎn)子劃分為軸段，本文軸系結(jié)構(gòu)共劃分為12個軸段，見圖2。

圖2 轉(zhuǎn)子軸段劃分示意圖

(2) 由于采用極轉(zhuǎn)動慣量原則進(jìn)行建模時，在保證總質(zhì)量及抗彎剛度的前提下，材料的密度、彈性模量數(shù)值會發(fā)生改變，因此需要對每個截面分別進(jìn)行特性的屬性定義。編寫截面子程序，繪制截面，定義截面不同位置的屬性，并劃分網(wǎng)格，使用SECWRITE命令定義截面名，保存后用于主程序調(diào)用；

(3) 編制主程序，在命令/PREP7定義渦輪軸向尺寸，使用BEAM189單元定義軸系，各個軸段的截面尺寸及網(wǎng)格調(diào)用截面子程序；軸承采用COMBIN14單元(如果是氣浮軸承或滑動軸承可選用COMBI214)；

(4) 在求解軸系的臨界轉(zhuǎn)速及模態(tài)振型時，在/SOLU命令后，通過MODOPT和MXPAND定義求解模態(tài)的方法及模態(tài)階數(shù)，求解結(jié)束后在/POST1命令后輸出結(jié)果；

(5) 在求解軸系的不平衡響應(yīng)時，于/SOLU命令后定義不平衡質(zhì)量的分布及掃描頻率范圍，然后通過/POST26命令獲取關(guān)鍵點的振動幅值及節(jié)點反力，在/POST1命令下輸出結(jié)果。

2.2 渦輪轉(zhuǎn)子的臨界轉(zhuǎn)速及模態(tài)振型分析

使用APDL建模方法，得到轉(zhuǎn)子的動力學(xué)模型，如圖3所示。

圖3 渦輪泵軸系模型

根據(jù)軸承生產(chǎn)廠家的剛度試驗數(shù)據(jù)，支承剛度取值6.3×107N/m，計算得到渦輪轉(zhuǎn)子坎貝爾圖(見圖4)及模態(tài)振型圖(見圖5、6)。

圖4 坎貝爾圖

圖5 一階振型圖 圖6 二階振型圖

渦輪轉(zhuǎn)子前三階臨界轉(zhuǎn)速為62 315.229、172 740.763、625 382.764 r/min。而超高速渦輪泵額定工作轉(zhuǎn)速為130 000 r/min。根據(jù)要求[13]，工作轉(zhuǎn)速與一階臨界轉(zhuǎn)速比應(yīng)當(dāng)大于1.33，與二階臨界轉(zhuǎn)速比應(yīng)當(dāng)小于0.8。此次研究的渦輪泵轉(zhuǎn)速比值為2.09、0.75，可以看出渦輪泵的工作轉(zhuǎn)速滿足要求，但是與二階臨界轉(zhuǎn)速之間的余量并不大，在特殊工況下，轉(zhuǎn)速升高時，會出現(xiàn)接近二階臨界轉(zhuǎn)速的情況。單方面的降低工作轉(zhuǎn)速確實能夠降低轉(zhuǎn)子共振帶來的影響，但這樣會嚴(yán)重影響渦輪泵的輸出性能。由模態(tài)振型圖，一階轉(zhuǎn)速下主要的振動為渦輪盤；二階轉(zhuǎn)速下主要的振動為泵葉輪位置。即不同轉(zhuǎn)速下渦輪盤和泵葉輪受振動激勵下作用下的響應(yīng)存在差別。

2.3 渦輪泵轉(zhuǎn)子不平衡響應(yīng)研究

由于軸系始終存在殘余不平衡質(zhì)量，高速旋轉(zhuǎn)下，將成為一個振動激勵。為了降低由于渦輪轉(zhuǎn)子旋轉(zhuǎn)帶來的振動影響，除了控制轉(zhuǎn)速避開各階臨界轉(zhuǎn)速以外，控制殘余不平衡質(zhì)量的大小和分布，也是一個方法。由模態(tài)振型圖可以看出，渦輪盤和泵葉輪的振動受不平衡質(zhì)量影響最大。文中研究不平衡質(zhì)量大小、分布位置作用下的軸系響應(yīng)。根據(jù)模態(tài)振型圖可知，在轉(zhuǎn)速小于二階臨界轉(zhuǎn)速的范圍內(nèi)，軸系的主要振動為渦輪盤和泵葉輪位置。分別簡記泵葉輪為位置A、渦輪盤為位置B；兩個軸承分別記為Ⅰ、Ⅱ，如圖7所示。

圖7 渦輪轉(zhuǎn)子示意圖

依據(jù)設(shè)計規(guī)范[14]，該渦輪泵轉(zhuǎn)子動平衡精度取G2.5，不平衡量計算：

式中：轉(zhuǎn)子總質(zhì)量M=0.35 kg，轉(zhuǎn)速取n=130 000 r/min；則最大允許不平衡量U=0.064 3 g·mm，通過表1方式對不平衡量進(jìn)行分配。

表1 不平衡質(zhì)量分布表

通過上述APDL計算流程，得到A、B點在不同頻率下的振動曲線。

由于渦輪泵額定工作轉(zhuǎn)速為130 000 r/min，通過圖8、9可以看出，在該轉(zhuǎn)速范圍內(nèi)，A點振動幅值不超過5 μm，B點不超過4 μm。受一階臨界轉(zhuǎn)速影響，轉(zhuǎn)速到達(dá)62 000 r/min時，A、B點振動曲線均出現(xiàn)峰值，該峰值隨著A點的不平衡質(zhì)量力減小出現(xiàn)一定程度的降低，但是當(dāng)不平恒質(zhì)量力繼續(xù)降低時，峰值并不呈降低趨勢。對于A點位置來說，其振動不完全取決于該位置的不平衡質(zhì)量力，還受到整個軸系振動的影響，即軸系其他位置(如B點)的振動同樣會傳遞至A點，圖8曲線為綜合影響下的振動。通過圖8、9各曲線，可以發(fā)現(xiàn)“分布4”即A點0.4U、B點0.6U條件下，一階臨界轉(zhuǎn)速下的峰值最小。

圖8 A點振動曲線圖

圖9 B點振動曲線圖

圖10 Ⅰ軸承反力曲線圖

圖11 Ⅱ軸承反力曲線圖

當(dāng)渦輪泵轉(zhuǎn)速超過160 000 r/min后，A、B點的振動急劇增大，一直到172 000 r/min，即二階臨界轉(zhuǎn)速，達(dá)到峰值并隨后逐漸降低。振動受不平衡質(zhì)量力影響與一階臨界轉(zhuǎn)速時一致。

對于軸承反力曲線來說，由于Ⅰ、Ⅱ軸承分別靠近A、B點位置，圖8、10，圖9、11完全一致。較大的振動幅值，將造成臨近軸承出現(xiàn)較大應(yīng)力(即支承反力)。

3 結(jié) 論

通過采用APDL參數(shù)化轉(zhuǎn)子動力學(xué)仿真技術(shù)，針對渦輪泵動力學(xué)進(jìn)行了研究，可以得出如下結(jié)論。

(1) 渦輪泵軸系在130 000 r/min轉(zhuǎn)速下工作，能夠有效避免一、二階臨界轉(zhuǎn)速影響，第三階臨界轉(zhuǎn)速由于遠(yuǎn)遠(yuǎn)高于工作轉(zhuǎn)速，可以忽略其對動力學(xué)的影響。

(2) 渦輪泵軸系轉(zhuǎn)子的一、二階模態(tài)振型中的主要振動位置分別為泵葉輪、渦輪盤，在使用過程中，可以通過采集兩個位置的振動以獲取軸系振動數(shù)據(jù)。

(3) 在G2.5平衡精度下，渦輪泵轉(zhuǎn)子在額定工況下工作，剩余不平衡質(zhì)量力的分布形式不會對軸系振動產(chǎn)生較大影響，當(dāng)不平衡質(zhì)量力分布接近A點0.4U、B點0.6U時，整個軸系振動受質(zhì)量力影響最小。

(4) 在渦輪泵結(jié)構(gòu)不變的情況下，如果需要提高工作轉(zhuǎn)速，需嚴(yán)格控制泵葉輪端的不平衡質(zhì)量力。

航天伺服燃?xì)鉁u輪泵的可靠性是伺服系統(tǒng)正常工作的保障，通過文中對渦輪泵軸系動力學(xué)的研究分析，一方面可以為后續(xù)渦輪泵產(chǎn)品設(shè)計提供參考，另一方面也可以為現(xiàn)有產(chǎn)品可靠性提升改進(jìn)提供方向。