蔣愛華,章 藝,靳思宇,章振華,黃修長,華宏星
(1.上海交通大學(xué) 機械系統(tǒng)振動國家重點實驗室,上海 200240;2.中國船舶重工集團公司第704研究所,上海 200031)
離心泵的振動是廣泛關(guān)注的問題,流體激勵力是離心泵振動的主要原因之一。
從離心泵振動的研究歷史來看,對流體激勵力的研究還相對較少?,F(xiàn)有對離心泵流體激勵力的研究主要集中于葉輪蝸殼模型的簡化以及葉輪在流體激勵力作用下的渦動情況[1-2]、作用于葉輪上流體激勵力的簡化以及此時離心泵轉(zhuǎn)軸的臨界轉(zhuǎn)速[3]、流體激勵作用下離心泵蝸殼的結(jié)構(gòu)輻射噪聲與內(nèi)部流體噪聲[4-6]等。
其中,對葉輪蝸殼模型的簡化還主要集中于二維模型[7]。這些模型對促進認識葉輪在流體中的運動過程有非常重要的作用,但也有其局限。一方面該模型未考慮帶動葉輪轉(zhuǎn)動的轉(zhuǎn)軸,無法研究轉(zhuǎn)軸與葉輪之間的相互作用力;另一方面二維模型也無法研究葉輪的陀螺力矩。
對作用于葉輪上流體激勵力的簡化以及此時離心泵轉(zhuǎn)軸的臨界轉(zhuǎn)速的研究,還主要是基于傳統(tǒng)轉(zhuǎn)子動力學(xué)方法對流體激勵力進行簡化,即將流體對葉輪的激勵力簡化為葉輪的附加質(zhì)量,該附加質(zhì)量為葉輪內(nèi)所包含流體20%-40%的質(zhì)量,而后計算此轉(zhuǎn)子模型的臨界轉(zhuǎn)速。這種方法可以不用計算真實的流體激勵力,而只需從結(jié)構(gòu)力學(xué)考慮系統(tǒng)的機械特性,因此工程中得到了廣泛的應(yīng)用[8]。
但這種簡化方法有明顯的缺陷:由于葉輪上增加了附加質(zhì)量,所計算的轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速與實際值將有較大偏差;另外同樣流量下,隨著葉輪轉(zhuǎn)速的變化,實際葉輪所受的流體激勵也將發(fā)生變化,而此簡化模型中明顯沒有考慮轉(zhuǎn)速不同附加質(zhì)量的變化;且由于葉輪為軸對稱結(jié)構(gòu)而蝸殼流域為非對稱流場,在葉輪轉(zhuǎn)動過程中所受流體激勵力也將為非軸對稱形式。因此有必要對這種模型進一步改進。
對流體激勵作用下離心泵噪聲的研究主要集中于葉片與流體流固耦合誘發(fā)噪聲簡化模型的建立、汽蝕時離心泵內(nèi)流場的噪聲、蝸舌與流體流固耦合誘發(fā)噪聲的實驗研究[9-11]。采用CFD計算出離心泵蝸殼內(nèi)表面的流體力,而后計算蝸殼向空氣中聲輻射的研究也得到了一定關(guān)注[6]。
但需要指出的是,由于工程中離心泵蝸殼常與帶動葉輪轉(zhuǎn)動的動力機械通過基座連接在一起,而流體激勵葉輪,而后通過轉(zhuǎn)軸傳遞至動力機械基座也將引起蝸殼的振動,因此將無法分辨引起蝸殼輻射噪聲的主要激勵源。并且離心泵與動力機械的共同基座常與其他板殼結(jié)構(gòu)相連接,而離心泵向空氣中輻射的噪聲可能遠小于這些板殼結(jié)構(gòu)所輻射的噪聲,此時離心泵向基座傳遞的振動加速則成為板殼結(jié)構(gòu)向外輻射主要原因。
現(xiàn)階段國內(nèi)外還未查到研究離心泵內(nèi)部流體對基座激勵的相關(guān)文獻。雖然運用CFD計算離心泵內(nèi)部流場的研究已經(jīng)比較廣泛和深入[12-13],但由于流體力學(xué)的研究者主要關(guān)注離心泵內(nèi)流場的變化。同時,機械振動的研究者也主要關(guān)注激勵作用下結(jié)構(gòu)的響應(yīng)與變形等[14]。流固耦合的研究者偏向于基于連續(xù)介質(zhì)力學(xué)建立流體與固體區(qū)域的控制方程[15],而后通過流固邊界位移相等條件聯(lián)系兩組方程,并運用數(shù)值方法求解流場與固體結(jié)構(gòu)各節(jié)點的位移,這種方法比單獨先求解流場而后將流固邊界上作用力用于固體結(jié)構(gòu)的單向流固耦合方法具有更加準確的結(jié)果,但由于這種方法有一定難度,且選擇流體控制方程時還主要應(yīng)用勢流中的伯努利方程,還無法滿足離心泵這種流場中含有流動分離湍流流動的旋轉(zhuǎn)機械的需要。
鑒于現(xiàn)有研究的不足以及通過試驗研究流體激勵力的困難,本研究將基于CFD技術(shù),應(yīng)用CFX對離心泵內(nèi)流場進行非定常分析,而后通過積分最終得出運轉(zhuǎn)過程中流體對離心泵結(jié)構(gòu)的激勵力。本研究共分為兩部分,分別闡述流體對蝸殼與葉輪的激勵作用,本文為前者。
對于不可壓常粘性流體湍流流動的離心泵內(nèi)流場,其控制方程為雷諾方程。張量形式的雷諾方程如式(1)所示。
SSTk-w模式是一種雙方程的湍流模式,在邊界層以內(nèi)采用湍動能方程k與湍動頻率方程w,在邊界層以外采用湍動能方程k與耗散率方程e,因此實際數(shù)值計算過程中總共需要求解的方程為六個。為了得到較準確的計算結(jié)果,采用SSTk-w模式時邊界層以內(nèi)需要有多層網(wǎng)格。
離散求解以上方程時,常采用有限體積法(FVM)[17-18]。本研究中采用的工具為CFX,是基于有限單元的有限體積法,更加容易得到收斂解。
所研究對象為一單級單吸式離心泵,其額定工況下比轉(zhuǎn)速為139。運用Por-E建立了離心泵內(nèi)各流體域模型,共建立了七個流體域,包括離心泵進口流域、離心泵葉輪內(nèi)流域、蝸殼內(nèi)流域、平衡室流域、葉輪前蓋板與蝸殼間隙流域(包括泄漏環(huán)處間隙)、葉輪后蓋板與蝸殼間隙流域與蝸殼出口延長段流域。其中葉輪流域為旋轉(zhuǎn)流域,而其余流域為靜止流域。
延長蝸殼出口流域是為了使流場能夠充分發(fā)展并達到穩(wěn)定狀態(tài),以更加準確計算出離心泵內(nèi)流場。另外,將平衡室流域也作為研究的對象是為了分析其對離心泵運轉(zhuǎn)過程中軸向力的平衡作用。
各流域網(wǎng)格的劃分均采用ICEM完成,所劃分網(wǎng)格均為六面體,以加快計算過程的迭代收斂速度與穩(wěn)定性。劃分網(wǎng)格時,為了能夠確保壁面底層網(wǎng)格在邊界層以內(nèi),通過預(yù)估計算得出所劃分網(wǎng)格的壁面當(dāng)量厚度Yplus值,并調(diào)整底層網(wǎng)格使得Yplus值處于11.63與500之間,以確保底層網(wǎng)格處于對數(shù)率層以內(nèi)[19-20],并對底層網(wǎng)格進行了加密。
各流域網(wǎng)格劃分完成后,導(dǎo)入CFX中通過Interface界面連接為一個整體,如圖1所示即為所劃分的離心泵內(nèi)流域網(wǎng)格,其中坐標(biāo)原點位于蝸殼對稱平面與葉輪軸心交點處,空間方向如圖中箭頭所示,表1為各流域網(wǎng)格、節(jié)點數(shù)量與體積。當(dāng)增加網(wǎng)格密度使得網(wǎng)格數(shù)量為表中網(wǎng)格數(shù)量的約1.5倍時,所得預(yù)估穩(wěn)態(tài)流場與表中網(wǎng)格數(shù)量下所獲得預(yù)估穩(wěn)態(tài)流場的蝸殼對稱平面壓力分布與速度矢量圖基本相同,則認為此時所得流場已為網(wǎng)格無關(guān)解。
為了能夠準確計算出葉輪轉(zhuǎn)動過程中離心泵內(nèi)流場的變化情況,需要給定準確的初值,因此首先計算穩(wěn)態(tài)流場以作為瞬態(tài)流場的初值,而穩(wěn)態(tài)流場計算過程中各參數(shù)初始值均為零。
圖1 離心泵流域網(wǎng)格Fig.1 Mesh of fluid volume in centrifugal pump
表1 各流域節(jié)點數(shù)、單元數(shù)與體積Tab.1 Nodes number,elements number and volume of different fluid parts
計算離心泵工況如圖2所示,葉輪轉(zhuǎn)速為額定轉(zhuǎn)速2 930 n/min,流量為額定流量100 t/h,其吸入口軸線與主軸垂直于水平面,吸入口與水面垂直距離為395 mm,水面有一個標(biāo)準大氣壓。因此在CFX中設(shè)置參考壓力為1 atm,進口邊界條件分別為靜壓3 871 Pa(亦即395 mm水柱),出口邊界條件設(shè)置為流量,其余邊界設(shè)為固壁無滑移條件,湍流模式選擇SST模式,計算過程中對式(2)積分所得公式進行二階中心差分,收斂標(biāo)準采用各參數(shù)平均殘差值小于1E-5。
圖2 離心泵工況示意圖Fig.2 Sketch map of centrifugal pump working condition
如表2中所示分五個階段計算了離心泵瞬態(tài)流場,第一個階段流場初始值為穩(wěn)態(tài)流場計算結(jié)果,其余階段初始值均以上一階段最后一個時間步所得結(jié)果。采用了時間步長與收斂標(biāo)準不同的五個階段來計算流場是為了盡量消除數(shù)值計算過程的不穩(wěn)定性與從穩(wěn)態(tài)流場計算中葉輪位置不變到葉輪轉(zhuǎn)動過程中流場的不穩(wěn)定性,并提高結(jié)果的準確性。
從表2中可以看出,第一階段的5個時間步中,葉輪每轉(zhuǎn)動0.4°為一個時間步,整個過程中葉輪流域相對靜止流域轉(zhuǎn)動了2°;第二階段的30個時間步中,葉輪每轉(zhuǎn)2°為一個時間步,葉輪流域相對靜止流域轉(zhuǎn)動了60°;第三階段的55個時間步中,葉輪每轉(zhuǎn)12°為一個時間步,葉輪流域轉(zhuǎn)動了660°;第四階段采用了與第二階段相同的時間步長計算了60步,葉輪流域共轉(zhuǎn)動了120°;最后階段的計算中采用與第四階段相同的時間步長,共計算了180步亦即葉輪流域轉(zhuǎn)動360°。其中,前四個階段迭代計算中以殘差小于1E-5為收斂標(biāo)準,而為了用更加準確的流場得出流體激勵力,第五個階段收斂標(biāo)準設(shè)置為殘差小于1E-6。因此,所有計算過程完成以后,葉輪流域共轉(zhuǎn)動了1 202°,共計330個時間步,此時葉輪以2 930 n/min的轉(zhuǎn)速轉(zhuǎn)動了0.068 5 s,葉輪每轉(zhuǎn)一度需要的時間約為 5.688 3e-5s。
如圖3示即為第四階段瞬態(tài)流場計算中,葉輪轉(zhuǎn)動60°所得XY截面的絕對壓力分布變化。所顯示的四個時刻中,每相鄰兩個時刻葉輪轉(zhuǎn)角相差20°,初始時刻XY平面中軸心葉尖連線與通過葉輪軸心的蝸舌圓弧切線夾角(以下簡稱蝸舌切線)約為14°,而蝸舌切線與X軸夾角為30°。從圖中可以看出第1步與第30步所得到的絕對壓力分布有較高的相似性。
表2 瞬態(tài)計算時間步設(shè)置Tab.2 Setup on time step of transient simulation
圖3 葉輪轉(zhuǎn)動60°過程中流場絕對壓力變化Fig.3 Changes on absolute pressure within 60 degrees of rotating impeller
表3 各方向上合力的統(tǒng)計參數(shù)Tab.3 Statistical parameters of forces composition on different directions
基于以上所得離心泵內(nèi)流場,運用CFX后處理模塊對蝸殼流域固體界面表面力進行積分,該表面力包括了流體對固體界面的總壓以及由固體界面處流體剪切運動所帶來的粘滯力,從而得出離心泵內(nèi)部固體界面所受三個方向的流體合力。如圖4示即為330步瞬態(tài)計算過程中蝸殼所受空間三個方向上流體合力的變化情況與所建立的蝸殼質(zhì)量-彈簧振動模型。
從圖4出,蝸殼各方向上流體合力在葉輪轉(zhuǎn)動過程中均圍繞一定值上下波動。特別是時間步大于90步以后的波形其波動周期更加明顯,且每30個時間步為一個波動周期,對應(yīng)葉輪轉(zhuǎn)角為60°,亦即葉輪每個葉片經(jīng)過蝸舌切線時各向合力會出現(xiàn)一次波動。另外,X向流體合力相對小于其他兩個方向上的合力,波形圍繞零值上下波動。Y向流體合力保持為負值,表明蝸殼所受Y向合力始終從蝸殼指向基座方向,亦即垂直于蝸殼出口與流體流出方向相反。Z向流體合力保持為正值,表明蝸殼所受Z向流體力始終從離心泵進口指向平衡室方向,亦即垂直于蝸殼進口方向與流體流入方向相同。
從圖4中還可以看出,各向合力從第36至第90個時間步的波形相對于90步以后的波形有更小的波動周期。這是因為第36至第90個時間步以葉輪流域每旋轉(zhuǎn)12°作為一個時間步,而第90個時間步以后均以葉輪流域每旋轉(zhuǎn)2°作為一個時間步,當(dāng)以時間步數(shù)為橫軸時,同樣時間步數(shù)前者葉輪轉(zhuǎn)角相當(dāng)于后者的六倍。
另外,圖中前35個時間步波形周期不明顯。這一方面是由于所對應(yīng)的葉輪轉(zhuǎn)角較小且時間步不一致,看不出波動周期,另一方面也是由于相鄰的時間步采用不同時間步長進行迭代計算時,將會出現(xiàn)數(shù)值的不穩(wěn)定,第35、第36個時間步與第90、第91個時間步處,Z向合力的兩個尖峰更明確的說明了這一點。因此,在進行蝸殼的流體激勵力響應(yīng)時,運用遠離時間步長發(fā)生變化處的流體激勵力將更加準確。
圖4 蝸殼三個方向上流體激勵合力Fig.4 Three direction inciting force on volute by fluid
以第五階段180個時間步所得結(jié)果進行頻譜分析,其功率譜如圖5所示,統(tǒng)計參數(shù)如表3所示。從圖5中可以看出,蝸殼各向流體合力均以293 Hz為基頻,與葉片通過頻率相同;與其它兩向流體合力相比,Y向合力具有更大的波動幅值,其二階頻率的譜值就已經(jīng)高于其余各向合力基頻譜值。
圖5 X、Y與Z向合力功率譜Fig.5 Power spectrum of forces composition on X,Y and Z directions
由2.1可知,葉輪以一定轉(zhuǎn)速轉(zhuǎn)動過程中,蝸殼所受各向流體合力均有明顯的周期變化,若能夠根據(jù)該規(guī)律建立不同轉(zhuǎn)角下蝸殼受流體力的數(shù)學(xué)模型,對進一步研究離心泵中的流體激振將有重要意義。由于第五階段所得結(jié)果更加準確,因此以該階段的180個時間步所得結(jié)果進行數(shù)學(xué)模型的建立。
為了使所建立模型將流體激勵力與葉輪的轉(zhuǎn)角相對應(yīng),本研究分別采用最小二乘多項式與傅里葉級數(shù)對原始信號進行逼近,所有過程均運用Matlab中完成現(xiàn)有函數(shù)完成。
首先將三個方向流體合力按照其周期進行平均,每30個時間步為一個周期,對應(yīng)葉輪轉(zhuǎn)角為60°,得出用于逼近的原始信號,而后采用九次多項式對所得信號進行了最小二乘的擬合。
運用傅里葉級數(shù)對流體三個方向上的激振力逼近時,采用未加窗FFT變換后信號的激振力基頻幅值作為正弦波的幅值,根據(jù)正弦波波峰與原始波形的波峰修正正弦波的相位。
三向流體激勵力的原始信號、九次多項式擬合波形與傅里葉級數(shù)波形如圖6(a)中所示。擬合信號與原始信號的標(biāo)準偏差如表4所示。從圖與表中可以看出,與基頻正弦波相比較,九次多項式擬合與原始信號有更小的偏差,這是因為一方面基頻正弦波只應(yīng)用了基頻幅值作為逼近波形的幅值,具有較大的能量損失,另一方面原始波形與標(biāo)準正、余弦波也相差較大。但同時從圖中也可以看出,九次項式擬合方法所得到的波形與原始信號仍然有較大的差距,并且多項式所需要擬合的階次也比較高,對下一步運用質(zhì)量彈簧模型計算離心泵蝸殼的振動帶來不便。因此,還需要尋找其他建模的方法。
表4 不同擬合信號與原始信號的標(biāo)準偏差Tab.4 Standard error between original signal and different fitting signals
從圖6(a)中還可以看出,三個波形均呈現(xiàn)出先起始點單調(diào)降低至波谷,而后由波谷單調(diào)上升至波峰,最后由波峰單調(diào)下降至末點的趨勢,因此如果采用分段曲線擬合則可以用較低的階次得到較好的結(jié)果。
圖6(b)所示即分別為對X向、Y向與Z向靜止流域固體界面合力的分段擬合。擬合時以起始點到波谷為第一段,波谷至波峰為第二段,剩余的為第三段,各段擬合是根據(jù)具體情況選擇所需要擬合的階數(shù),并且各橫軸也由時間轉(zhuǎn)化為對應(yīng)的葉輪轉(zhuǎn)角。從圖中可以看出X向流體合力波谷與波峰所對應(yīng)的轉(zhuǎn)角為18°與42°處,由于初始時刻XY平面的軸心葉尖連線與蝸舌切線之間的夾角為14°,亦即X向流體合力最小與最大值分別出現(xiàn)于葉尖繞軸心轉(zhuǎn)動通過蝸舌4°與28°處;類似地Y向流體合力幅值最大出現(xiàn)于葉尖處于軸心蝸舌連線時,最小幅值出現(xiàn)于軸心葉尖連線繞葉輪軸心轉(zhuǎn)動通過蝸舌切線后14°處;Z向流體合力最小值對應(yīng)的葉輪轉(zhuǎn)角與Y向流體合力最大幅值處葉輪轉(zhuǎn)角相同,但最大值對應(yīng)的葉輪轉(zhuǎn)角與X向最大值處葉輪轉(zhuǎn)角相同。
圖6 蝸殼三向流體激力擬合Fig.6 Polynomial fitting of three-direction composite fluid forces on volute
從圖6(b)中可以看出,各向合力在分三段擬合以后的結(jié)果明顯好于圖6(a)中所示擬合結(jié)果,擬合信號與原始信號的標(biāo)準偏差分別為 1.09、3.51 與 1.31,也遠小于表4中所示兩種方法所獲得的標(biāo)準偏差,并且分三段擬合時各段多項式最高僅為三階,三個方向上流體激勵合力波形分三段擬合所獲得多項式分別如式(3)、式(4)與式(5)所示。
(1)運用CFD計算了離心泵葉輪轉(zhuǎn)動過程中的瞬態(tài)流場,并基于瞬態(tài)流場積分得出了蝸殼所受流體激勵合力,結(jié)果表明,離心泵運轉(zhuǎn)過程中蝸殼受到空間三個方向葉片通過頻率的周期性流體激勵力。
(2)對比分析整個計算過程中蝸殼流體合力并進行功率譜分析表明,三個方向的流體合力中,垂直于蝸殼出口平面方向的合力與流體流出方向相反并具有最大的波動幅值,垂直于蝸殼進口平面方向的合力與流體流入方向相同,正交于前兩個方向的合力方向隨著波動周期而改變并具有最小的波動幅值。
(3)對比分析單個周期內(nèi)各向流體合力變化可知,垂直于蝸殼出口表面方向的流體合力與垂直蝸殼進口表面方向的流體合力的波谷值,均出現(xiàn)于葉尖轉(zhuǎn)過蝸殼對稱平面內(nèi)穿過葉輪軸心的蝸舌圓弧切線處。
(4)分別運用九次多項式擬合、波動基頻正弦波與三段多項式擬合近似各向周期性流體激勵力,結(jié)果表明采用三段多項式擬合所建的數(shù)學(xué)模型與原始波形有最小的偏差,并且具有較低階次。
[1]Adkins D R.Analyses of hydrodynamic forces on centrifugal pump impellers.PHD dissertation[D].California Institute of Technology,1986.
[2]Adkins D R,Brennen C E.Analyses of hydrodynamic radial forces on centrifugal pump impeller[J].Transaction of ASME,1988(3):20-28.
[3]袁振偉,楮福磊,林言麗,等.考慮流體作用的轉(zhuǎn)子動力學(xué)有限元模型[J].動力工程,2005(8):457-461.
[4]Langthjem M A,Olhoff N.A numerical Study of Flowinduced noise in a two-dimensional centrifugal pump[J].Part I:Hydrodynamics.Journal of Fluids and Structures,2004,(19):349-368.
[5]Langthjem M A,Olhoff N.A numerical study of flow-induced noise in atwo-dimensionalcentrifugalpump.PartII:Hydroacoustics[J].Journal of Fluids and Structures,2004,(19):369-386.
[6] Jiang Y Y,Yoshimura S,Imai R,et al.Quantitative evaluation of flow-induced structural vibration and noise in turbomachinery by full-scale weakly coupled simulation[J].Journal of Fluids and Structures,2007(23):531-544.
[7]Brennen C E.Hydrodynamics of pumps[M].Oxford University Press,1994.
[8]聞邦春,顧家柳,夏松波,等.高等轉(zhuǎn)子動力學(xué)[J].機械工業(yè)出版社,2000.
[9]Howe M S.On the estimation of sound produced by complex fluid-structure interactions,with application to avortex interacting with a shrouded rotor[J].Mathematical and Physical Sciences,1991(443):573-598.
[10] Schmitz S.Reducing Pump noise in cooling tower applications[J].World Pumps,2004(9):24-29.
[11] Dong R,Chu S,Katz J.Effect of modification to tongue and impeller geometry on unsteady flow,pressure fluctuations,and noise in a centrifugal pump[J]. Journal of Turbomachinery,1997(119):506-515.
[12]郭鵬程.水力機械內(nèi)部復(fù)雜流動的數(shù)值研究與性能預(yù)測[D].西安:西安理工大學(xué),2006,9.
[13]徐朝暉.高速離心泵內(nèi)全流道三維流動及其流體誘發(fā)壓力脈動研究[D].北京:清華大學(xué),2004,4.
[14]李德葆,蘇銘德,諸葛洪程,等.流體激發(fā)結(jié)構(gòu)振動響應(yīng)預(yù)測模型及激振力的識別[J].振動與沖擊,1991,10(1):6-13.
[15]邢景棠,周 盛,崔爾杰.流固耦合力學(xué)概述[J].力學(xué)進展,1997(2):19-38.
[16] Menter F R.Two-equation eddy-viscosity turbulence models forengineering applications[J]. American Institute of Aeronautics and Astronautics.1994(8):1598-1605.
[17] Eymard R,Gallouet T,Herbin R.Finite volume method[M].Handbook of Numerical Analysis,1997.
[18] Ferziger J H,Peric M.Computational methods for fluid dynamics[M].Springer.-Verlag.2001.
[19] Fluent Inc.,User's Guide[M].Fluent Inc,2003.
[20] Versteeg H K, Malalasekera W. An introduction to computational fluid dynamics:the finite volume method[M].Wiley,New York,1995.