基于Hermitian小波有限元的葉輪葉片載荷識別

薛曉峰，陳雪峰，耿 佳，張興武

(西安交通大學 機械制造系統工程國家重點實驗室，西安 710049)

離心壓縮機的關鍵旋轉部件是葉輪。葉輪所承受的氣動載荷對離心式壓縮機的整體性能起著關鍵性的決定作用。葉輪結構的復雜、載荷多樣性及工作環境的惡劣對離心式壓縮機的安全性提出了很大的挑戰。工作過程中，葉輪葉片表面受到離心力、氣動激振力、各種應力及復雜的介質作用沖擊力等等各種作用的影響，這些復雜的載荷屬于典型的隨機激勵，因而各種復雜的載荷沖擊是造成葉輪破壞的主要原因。葉輪工作中存在的輸運介質如固體顆粒等，由于葉輪內部的高速氣流帶動，以高速，變角度的方式沖擊葉片，導致葉片變薄，甚至斷裂的事故。1997年404廠的透平壓縮機組二級葉輪的端蓋與葉片在顆粒介質的沖擊力作用下產生嚴重的沖刷與沖蝕腐蝕。2007年某石化廠的高速離心式壓縮機，二級葉輪葉片在復雜載荷的沖擊下發生葉片斷裂事故。2008年某部門一臺壓縮機一級葉輪葉片葉根處在小顆粒的長期沖擊下掉塊，導致葉輪葉片前緣發生嚴重變形。2012年沈陽鼓風機廠離心式壓縮機運行兩年后檢修，發現大葉片壓力面后根部發生嚴重磨損[1]。

由于復雜的交變載荷對葉輪產生著重要的作用。國內外對離心式壓縮機葉輪葉片的載荷也開展了很多研究。如Shibata等[2-3]通過研究葉片載荷對離心式壓縮機性能的影響，對葉輪結構進行了改進和優化設計，提高了葉輪的效率。Krain等[4]對離心式壓縮機的葉輪進行優化，通過激光測試技術測量了葉輪的內部氣動載荷流動場，至今仍有很多研究人員用Krain葉輪數據來驗證自己的數值模擬。Ubaldi等[5]也對離心式壓縮機的葉輪進行了實驗研究，測量了載荷的徑向速度、周向速度等，得到了葉輪內部的載荷流動規律，耗散特點。國內的胡小文等對離心式壓縮機的葉輪載荷分別進行了數值研究，表面載荷集中在葉片前端的葉輪性能優于載荷集中于葉片中后端。陶麗樺等[6]指出葉片后端的載荷峰值可以得到較高的氣動性能。

研究表明，葉片上的流體激勵引起的葉片振動導致葉片的疲勞破壞，而引起葉片疲勞破壞的載荷為離心式壓縮機所承受的氣動載荷。而葉輪內部交變載荷引發的流體激振主要是由葉輪機械受到的強迫振動和自激振動引起的。對葉輪內部的流場分析，主要是由流體力學解決的。由于實驗的復雜，同時伴隨著計算機的快速發展，針對離心式壓縮機葉輪這種復雜的旋轉機械，利用有限元軟件進行結構分析適應于復雜的幾何結構并具有可靠的理論基礎，便于快速實現。Toshio提出對葉輪有限元分析時分為葉片、輪盤和輪蓋三部分進行，將離心力轉化為軸對稱的徑向作用力進行分析。針對各種問題的求解，葉輪葉片被簡化為各種模型進行模擬，均得到比較好的結果。

葉輪葉片作為離心式壓縮機的核心部件，工作過程中載荷的作用又如此重要，離心式壓縮機葉輪葉片的載荷識別以某廠提供的葉輪為研究對象，利用商用軟件ANSYS進行建模，而在關鍵區域，采用高精度的小波有限元進行分析，以提高載荷識別的精度。利用第五章提出的混合單元方法實現對葉輪葉片的載荷識別技術，此方法對改進葉輪葉片的設計起著一定的理論基礎和重要的實踐意義。

離心式壓縮機葉輪葉片所受的載荷復雜性、隨機性、不可測量特性，及其葉輪工作過程中的惡劣的環境如高速、高壓、腐蝕、疲勞等。研究離心式壓縮機葉輪葉片所受的載荷對提高離心式壓縮機的性能有重大的意義。而載荷的不可測量使得研究者往往無法深入研究其對離心式壓縮機造成的危害。當然，對離心式壓縮機的振動響應測量是比較容易的事情。有限元載荷識別的算法就是通過逆Newmark算法構造載荷激勵與振動響應直接的傳遞矩陣，然后利用傳感器測量出來的振動響應對葉輪葉片所受的載荷進行反求。有限元載荷識別技術為識別離心式壓縮機葉輪葉片所受的載荷提供了一個新的思路。通過采用有載荷識別技術，識別出離心式壓縮機葉輪葉片所受的載荷，對研究離心式壓縮機具有非常重要的意義。當然，針對傳統商業軟件的效率低、收斂慢等等因素，從提高分析精度的角度出發，利用高精度的Hermitian小波有限元[7-8]對葉輪葉片進行局部分析，而對整體結構使用商業軟件ANSYS進行分析。使用小波單元嵌入商業軟件中所組成的混合有限元對葉輪葉片進行載荷識別的分析。

1 Hermitian小波有限元及載荷識別理論

通過構造Hermitian小波殼單元，將小波單元的剛度矩陣和質量矩陣代入逆Newmark算法中，建立激勵載荷和響應的傳遞矩陣，然后根據響應和傳遞矩陣來重構載荷。首先介紹Hermitian小波殼單元的插值函數。

1.1 Hermitian小波插值函數

Hermitian小波尺度函數φ1,k(x)，如圖1(a)所示并且表達式為

(1)

由于邊界條件的端點效應，Hermitian小波尺度函數無法直接作為有限元插值函數進行插值。針對這種情況，對原始的Hermitian小波尺度函數進行改造，修正的Hermitian小波尺度函數滿足有限元插值的條件，并且這種修正保留了Hermitian小波作為插值函數的一切優良特性。

(2)

Hermitian小波函數ψj,k，如圖2所示并且表達式為

(3)

圖1 在空間V1的尺度函數Fig.1 Scaling functions in scaling space V1

Hermitian小波基函數的特性是

(4)

并且

(5)

在空間Vj的小波基可以寫為

Φj={φ1,ψ1,ψ2,…,ψj-1}

(6)

圖2 在空間W的小波函數Fig.2 Wavelets in wavelet space W

圖3 Hermitian插值函數的張量積Fig.3 Tensor product HCSWI interpolation function

1.2 Hermitian小波殼單元

Hermitian小波殼單元可以看做是Hermitian小波彈性板單元[9]與中厚板單元[10]的組合。如圖4所示

節點自由度

(7)

剛度矩陣

(8)

質量矩陣

(9)

圖4 殼單元求解域及節點自由度排列Fig.4 Rectangle solving domain inshell element.

1.3 載荷識別公式

Newmark算法是一種隱式時間步積分方法，廣泛用于結構的動態分析。響應作為已經量來求解載荷。通過求解激勵點和響應點的傳遞矩陣來建立載荷識別算法。

Newmark算法是基于時間 [t,t+Δt]的假設

(10)

(11)

重寫等式(11)，加速度公式為

(12)

線性結構的運動微分方程可以表達為

(13)

式中：M是整體剛度矩陣，C是整天阻尼矩陣，K是整體剛度矩陣，F是整天激勵向量。

結合等式(10)，式(12)和式(13)，我們可以得到

(14)

等式(14)可以寫為矩陣形式

(15)

重寫等式(15)，m點位置的激勵和n點位置的響應可以寫為

(16)

求解激勵的病態矩陣采用截斷奇異值方法來處理。

載荷識別采用均方根值誤差(RMS)來表示

(17)

2 離心式壓縮機葉輪葉片載荷識別實驗

如圖5所示的離心式壓縮機葉輪是研究對象，具體思路是使用我們提出的新的Hermitian小波殼單元對葉片的待識別載荷區域進行建模，然后計算小波單元的剛度矩陣、質量矩陣。針對葉輪葉片的整體模型，我們使用商業軟件ANSYS進行建模，提取葉輪葉片的剛度矩陣、質量矩陣。使用小波單元計算出的的剛度矩陣、質量矩陣替換葉輪葉片中待識別載荷區域的剛度矩陣、質量矩陣，組成混合剛度矩陣。通過將混合剛度矩陣、質量矩陣代入逆Newmark算法中求解葉輪葉片載荷點與葉輪端蓋中測量加速度點的傳遞矩陣。最后根據傳遞矩陣和測量的加速度信號求解葉片的沖擊載荷。載荷識別的技術路線如圖6所示。Hermitian小波殼單元嵌入商業軟件ANSYS中構造出的混合單元可以高精度的反求葉輪葉片的激勵載荷。

圖5 離心式壓縮機葉輪Fig.5 The impeller of centrifugal compressor

測試系統為標準的錘擊測試，該葉輪是國內某注明壓縮機制造商加工而成，葉輪固定在轉子實驗平臺上，流線頭對葉輪進行軸向約束，如圖7所示。數采設備選擇德國HBM Genesis 2i 高速便攜式可視化數據記錄儀，使用力錘進行敲擊測試及PCB單軸加速度傳感器測試響應，力錘錘頭為鋁制錘頭，靈敏度為12.85 mV/N，加速度傳感器靈敏度為100 mV/g。系統采樣頻率是10.24 kS/s。葉輪中受載荷的部位主要是葉片，從實驗的可行性出發，通過對葉片的端部進行力錘敲擊獲取激勵的信號。對于錘擊測試，激勵方式是用錘擊法敲擊沖擊激勵，通過固定響應點，變換激勵點的方式來敲擊，反求載荷的方式是通過構造出傳感器與錘擊激勵的傳遞矩陣，然后利用傳感器獲得的加速度信號求解出激勵載荷。沖擊力錘擊法的敲擊點與傳感器響應點的布置如圖8所示。

圖6 離心式壓縮機葉輪葉片載荷識別技術路線Fig.6 Load identification technology roadmap for centrifugal compressor impeller blade

圖7 離心式壓縮機葉輪葉片載荷識別測試系統Fig.7 Load identification system for centrifugal compressor impeller blade

圖8 離心式壓縮機葉輪葉片載荷識別激勵與響應的布置方式Fig.8 The load identification arrangement about excitationand response for centrifugal compressor impeller blade

葉輪的材料參數是彈性模量214 600 MPa，剪切模量82 500 MPa，泊松比0.300 6，密度7 850 kg/m3。葉輪葉片所受的載荷一般都是在葉片端部比較嚴重，采取用力錘敲擊獲得沖擊力，然后用傳感器測試獲得加速度響應，通過代入加速度響應進逆Newmark算法中，反求激勵，與力錘敲擊的沖擊力相比較，用實體模型對Hermitian小波載荷識別算法進行驗證。傳感器貼在葉輪后端蓋的前端，選擇8個敲擊點用力錘進行敲擊，其中，敲擊點7是第一個大葉片的下端部，敲擊點2、8是第二個大葉片的上下端部，敲擊點1、3、5是三個小葉片的上端部，敲擊點4、6是第三、第四個大葉片的上端部，由于小葉片的下端部無法用力錘敲擊測試，因而測試的是大葉片的上、下端部和小葉片的上端部。

考慮到葉輪葉片均為薄壁結構，其中大葉片、小葉片厚度是9 mm，前端蓋厚度是13.7 mm，后端蓋厚度是16 mm，內端蓋厚度是22.08 mm。使用SolidWorks三維軟件分別抽去葉輪葉片厚度，得到的二維葉輪葉片如圖9(a)所示。將二維葉輪導入商業軟件ANSYS中，當然，導入的模型需要修整，刪除小面，短線等，ANASYS單元選擇SHELL63單元，SHELL63單元既具有彎曲能力又具有膜力，可以承受平面內載荷和法向載荷。本單元每個節點具有6個自由度：沿節點坐標系X，Y，Z方向的平動和沿節點坐標系X，Y，Z軸的轉動。網格模型中，共劃分16 931個四邊形單元，18 310個節點，如圖9(b)所示。這里的前端蓋、后端蓋及內端蓋因為不是分析的重點，所以劃分的網格數較少，但是也確保使用四邊形網格進行掃掠劃分。大葉片、小葉片是主要承受載荷的部位，網格數劃分的相對比較多。二維葉輪葉片在商業軟件中被拆分成25個零部件，分別是前端蓋1個、后端蓋1個、內端蓋1個、大葉片11個、小葉片11個，全部采用綁定接觸，這里需要說明的是，將二維葉輪拆分是為了劃分網格的方便，而綁定接觸相當于焊接的一種接觸方式，這種選擇也是比較符合實際的。約束方式采用將內端蓋與后端蓋的內徑表面全部固定約束，這種約束是為了與實驗測試的約束方式保持一致。

圖9 二維離心式壓縮機葉輪及網格劃分Fig.9 The impeller and meshing for the two-dimensional

使用新的的Hermitian小波有限元取代商業軟件ANSYS的部分結構，進行載荷識別的分析。葉輪葉片ANSYS的節點剛度矩陣和質量矩陣通過Harwell-Boeing方法來提取。每個節點提取6個自由度的剛度矩陣和質量矩陣。葉片屬于承受載荷的關鍵部件，因此Hermitian小波單元會取代大葉片、小葉片的部分SHELL63單元，由于一個Hermitian小波單元有18×18個內節點，等價于商業軟件ANSYS的17×17個單元。為方便用一個小波單元替換商業軟件的17×17個單元，將小波單元的節點號排列按照商業軟件的單元節點排列方式重拍，這樣小波單元可以按照節點替換的方式直接替換商業軟件的相應節點。如圖10 所示為商業軟件ANSYS17×17個單元和1個Hermitian小波有限元的節點號排列方式。

圖10 有限元的節點號排列方式Fig.10 The node number arrangement based the finite element method

由于一個節點有6個自由度，因此小波單元相應的節點剛度矩陣和質量矩陣可以直接替換商業軟件ANSYS對應結構節點的剛度矩陣和質量矩陣。替換節點的方式是從商業軟件ANSYS提取的總剛度矩陣和質量矩陣中減掉大葉片前端、后端，小葉片前端的節點剛度矩陣和質量矩陣，然后添加相應的小波殼單元的剛度矩陣和質量矩陣。用1個Hermitian小波殼單元來分別替換大葉片前后端部和小葉片前端部的17×17部分的網格。

大葉片劃分969個四邊形單元，1 044個節點，如圖11所示。對大葉片劃分為17×47的網格。用兩個Hermitian小波單元分別替換大葉片網格的前端、后端SHELL63單元中劃分為17×17部分。通過將兩個Hermitian小波單元替換大葉片的前后端的SHELL63單元組成新的混合單元，混合單元代入逆Newmark算法構造敲擊點和傳感器測量的加速度響應點的傳遞矩陣，傳遞矩陣和加速度響應可以反求出大葉片前后端的激勵載荷。力錘敲擊的沖擊載荷作為精確解驗證算法的正確性。

圖11 離心式壓縮機二維大葉片及網格劃分Fig.11 The two-dimensional large blade and meshing for the centrifugal compressor

小葉片劃分510個四邊形單元，558個節點。如圖12所示。對小葉片劃分為17×30的網格。因此，用一個Hermitian小波單元替換小葉片的前端SHELL63單元中劃分為17×17部分。通過將一個Hermitian小波單元替換小葉片的前端的SHELL63單元組成新的混合單元，混合單元代入逆Newmark算法構造敲擊點和傳感器測量的加速度響應點的傳遞矩陣，傳遞矩陣和加速度響應可以反求出大葉片前后端的激勵載荷。力錘敲擊的沖擊載荷作為精確解驗證算法的正確性。

圖12 離心式壓縮機二維小葉片及網格劃分Fig.12 The two-dimensional small blade and meshing for the centrifugal compressor

結構阻尼矩陣采用線性Rayleigh阻尼

C=a1M+a2K

(18)

式中：a1=0.180 7和a2=6.8×10-4是阻尼系數[11]。

3 離心式壓縮機葉輪葉片載荷識別結果分析

沖擊實驗中，沖擊激勵是持續的記錄。使用力錘去錘擊并且測量載荷，加速度響應是通過加速度傳感器去測量。如圖13所示是在敲擊點7的沖擊激勵和相應的加速度響應。錘擊激勵是典型的沖擊載荷，相應的加速度信x號顯示出快速的衰減趨勢。每個敲擊點錘擊五次，分別截取一次錘擊激勵的加速度響應進行載荷識別。

(19)

如圖7所示，將敲擊點7所在的大葉片稱為第一個大葉片，那么敲擊點1所在的小葉片就是第一個小葉片，敲擊點2、8所在的大葉片是第二個大葉片，敲擊點3所在的小葉片是第二個小葉片，敲擊點4所在的大葉片是第三個大葉片，敲擊點5所在的小葉片是第三個小葉片，敲擊點6所在的大葉片是第四個大葉片。

如圖14所示，敲擊點7是第一個大葉片的下端部，由于傳感器所在的位置距離敲擊點7比較近，所以Hermitian小波有限元和商業軟件ANSYS識別的精度都比較高，如圖14(a)中所示，Hermitian小波有限元精度高點，圖14(b)中誤差圖也說明了Hermitian小波有限元的誤差是商業軟件ANSYS的1/4左右。

圖13 力錘和加速度傳感器的測量數據Fig.13 The measurement data from the hammer and the acceleration sensor

圖14 敲擊點7的載荷識別Fig.14 Load identification at node 7

如圖15所示，敲擊點1是第一個小葉片的上端部，傳感器所在的位置距離敲擊點1略微遠點，但是Hermitian小波有限元和商業軟件ANSYS仍可以高精度的識別出沖擊載荷，圖15(b)中誤差圖顯示Hermitian小波有限元的誤差是商業軟件ANSYS的1/4左右。

如圖16所示，敲擊點2是第二個大葉片的上端部，傳感器所在的位置距離敲擊點2接近于傳感器距離敲擊點1的距離，Hermitian小波有限元和商業軟件ANSYS可以比較精確的識別出沖擊載荷，圖16(b)中誤差圖顯示Hermitian小波有限元與商業軟件ANSYS與圖16(b)相比，略微有點增加，且Hermitian小波有限元的誤差是商業軟件ANSYS的1/3左右。

圖15 敲擊點1的載荷識別Fig.15 Load identification at node 1

圖16 敲擊點2的載荷識別Fig.16 Load identification at node 2

如圖17所示，敲擊點8是第二個大葉片的下端部，傳感器所在的位置距離敲擊點8略微遠于傳感器距離敲擊點7的距離，Hermitian小波有限元和商業軟件ANSYS可以比較精確的識別出沖擊載荷，圖17(b)中誤差圖顯示Hermitian小波有限元與商業軟件ANSYS的誤差比較接近于敲擊點1的誤差，這是因為不同的敲擊點與傳感器所在位置所構成的傳遞矩陣不同，總體來說，激勵與傳感器的位置越近，識別出的載荷精度越高。Hermitian小波有限元的誤差是商業軟件ANSYS的1/4左右。

如圖18所示，敲擊點3是第二個小葉片的上端部，傳感器所在的位置距離敲擊點3略微遠于傳感器距離敲擊點2的距離，Hermitian小波有限元和商業軟件ANSYS可以識別出沖擊載荷，圖18(b)中誤差圖顯示Hermitian小波有限元與商業軟件ANSYS的誤差略微高于敲擊點2的誤差，這是由于敲擊點3與傳感器的距離遠與敲擊點2與傳感器的距離造成的。Hermitian小波有限元的誤差是商業軟件ANSYS的1/4左右。

如圖19所示，敲擊點4是第三個大葉片的上端部，傳感器所在的位置距離敲擊點4略微遠于傳感器距離敲擊點3的距離，Hermitian小波有限元和商業軟件ANSYS識別沖擊載荷的能力進一步減弱，圖19(b)中誤差圖顯示Hermitian小波有限元與商業軟件ANSYS的誤差略高于敲擊點3的誤差，這是由于敲擊點4與傳感器的距離遠與敲擊點3與傳感器的距離造成的。Hermitian小波有限元的誤差是商業軟件ANSYS的1/4～1/3之間。

圖17 敲擊點8的載荷識別Fig.17 Load identification at node 8

圖18 敲擊點3的載荷識別Fig.18 Load identification at node 3

圖19 敲擊點4的載荷識別Fig.19 Load identification at node 4

如圖20所示，敲擊點5、6分別是第三個小葉片的上端部、第四個大葉片的上端部，傳感器所在的位置距離敲擊點5、6相比較傳感器距離敲擊點4的距離，又遠了許多，Hermitian小波有限元和商業軟件ANSYS識別沖擊載荷幾乎被噪聲淹沒了。無論是Hermitian小波有限元還是商業軟件ANSYS，識別沖擊載荷的能力都是隨著傳感器與載荷激勵的位置不斷接近而增強的。當載荷激勵的位置距離傳感器的位置比較遠時，有限元識別出的載荷就會被噪聲淹沒，失真。

圖20 敲擊點5、6的載荷識別Fig.20 Load identification at node 5,6

4 結 論

離心式壓縮機葉輪葉片的載荷識別，采用Hermitian小波有限元與逆Newmark算法構造力錘敲擊點與加速度傳感器所在位置的傳遞矩陣，根據加速度響應信號重構載荷，案例可以清晰的看見Hermitian小波有限元識別的載荷，誤差是商業軟件ANSYS的1/4～1/3。這種精度在工程中的意義很大。而不論是Hermitian小波有限元還是商業軟件ANSYS，當載荷激勵的位置距離傳感器的位置越接近時，利用加速度信號重構出來的載荷精度就越高。

所采用的小波有限元載荷識別技術，是通過力錘敲擊沖擊力，傳感器測試加速度響應來反求載荷，相對比較簡單，可以為離心式壓縮機葉輪葉片的載荷識別提供一定的思路。實際工況中的葉輪葉片所受的載荷非常復雜，且實際工況中葉輪葉片的載荷位置很難獲取，振動響應的測量也非常困難，這給載荷識別帶來很大的難度，目前采取的方式是將傳感器固定在輪盤外殼上測量振動響應，當然使用這種信號去識別載荷還有待于進一步研究。

參 考 文 獻

[1] 王光存.離心壓縮機葉輪沖蝕磨損機理和規律的研究[D].濟南：山東大學, 2015.

[2] SHIBATA T, YAGI M, NISHIDA H.Effect of impeller blade loading on compressor stage performance in a high specific speed range[J].Journal of Turbomachinery, 2010, 134(4): 1757-1768.

[3] SHIBATA T, YAGI M, NISHIDA H.Performance improvement of a centrifugal compressor stage by increasing degree of reaction and optimizing blade loading of a shrouded 3D impeller[J].Journal of Turbomachinery, 2009, 133(2): 1305-1315.

[4] KRAIN H.A CAD-method for centrifugal compressor impellers[J].Journal of Engineering for Gas Turbines & Power, 1984, 106(2): 482-488.

[5] UBALDI M, ZUNINO P, GHIGLIONE A.Detailed flow measurements within the impeller and the vaneless diffuser of a centrifugal turbomachine[J].Experimental Thermal & Fluid Science, 1998, 17(1/2): 147-155.

[6] 陶麗樺, 王鍵, 劉艷.葉片載荷分布對離心葉輪氣動性能影響的數值研究[J].風機技術, 2015, 57(6): 13-20.

TAO Lihua, WANG Jian, LIU Yan.Effect of the blade loading distribution on aerodynamic performance of centrifugal impeller[J].Compressor Blower & Fan Technology, 2015,57(6): 13-20.

[7] XUE Xiaofeng, ZHANG Xingwu, LI Bing.Modified hermitian cubic spline wavelet on interval finite element for wave propagation and load identification[J].Finite Elements in Analysis & Design, 2014, 91(15): 48-58.

[8] ZHANG X W, GAO R X, YAN R Q.Multivariable wavelet finite element-based vibration model for quantitative crack identification by using particle swarm optimization[J].Journal of Sound & Vibration, 2016, 375: 200-216.

[9] XUE Xiaofeng, CHEN Xuefeng, ZHANG Xingwu.Hermitian plane wavelet finite element method: Wave propagation and load identification[J].Computers & Mathematics with Applications, 2016, 72(12): 2920-2942.

[10] XUE Xiaofeng, CHEN Xuefeng, ZHANG Xingwu.Hermitian mindlin plate wavelet finite element method for load identification[J].Shock and Vibration, 2016(6): 1-24.

[11] JEHEL P, LéGER P, IBRAHIMBEGOVIC A.Initial versus tangent stiffness-based Rayleigh damping in inelastic time history seismic analyses[J].Earthquake Engineering & Structural Dynamics, 2014, 43(3): 467-484.

[12] YANG Zhibo , CHEN Xuefeng , XIE Yong.Hybrid two-step method of damage detection for platelike structures[J].Structural Control & Health Monitoring, 2015, 23(2): 267-285.