何 曦 徐勝利 劉海濤 王曉放 陳旭東

（1.大連理工大學(xué)能源與動力學(xué)院；2.南洋理工大學(xué)計算機科學(xué)與工程學(xué)院）

0 引言

核主泵葉輪是為冷卻劑輸送提供動力并將熱量排出的設(shè)備，是核電站一回路壓力邊界組成部分。葉輪是核主泵內(nèi)高速轉(zhuǎn)動部件，同時作為核主泵的核心部件關(guān)系到核主泵運行的完整性，要求在高溫、高壓、強輻射環(huán)境下長期、安全、可靠地運行，直接影響冷卻劑的輸送，進(jìn)而影響整個反應(yīng)堆系統(tǒng)的安全。葉輪強度分析的關(guān)鍵因素包括葉輪模態(tài)分析、葉輪受到水壓力脈動載荷特性、葉輪共振校核和振動響應(yīng)分析。

核主泵葉輪在水中的模態(tài)參數(shù)與空氣中模態(tài)參數(shù)不同。由于周圍水附加質(zhì)量的影響，在水中的模態(tài)頻率比空氣中的頻率要低。孔銘等對葉輪干濕模態(tài)進(jìn)行了實驗研究，分析對比了葉輪干濕模態(tài)頻率，得到由于附粘水的影響，每階濕模態(tài)頻率低于相應(yīng)干模態(tài)頻率的結(jié)論[1]。張新等基于葉輪水中流固耦合對軸流泵葉輪進(jìn)行了模態(tài)分析，分析了葉輪在水中以及空氣中的各階頻率以及振型，研究發(fā)現(xiàn)各階固有頻率下降系數(shù)范圍在0.1～0.38，下降系數(shù)與各階振型有關(guān)[2]。Liang[3]等討論了附加質(zhì)量對于模態(tài)的影響。鄭小波[2]等采用Subspace以及Block Lancos方法求解了空氣中模態(tài)，并且利用流固耦合求解了水中模態(tài)，而且對頻率、振型等模態(tài)參數(shù)進(jìn)行了對比。張學(xué)榮[5]等通過對葉片干濕模態(tài)有限元分析，比較了GE公司的水中頻率下降率修正系數(shù)與實際每階模態(tài)頻率下降率。

葉輪在水環(huán)境下受到的阻尼比可以分為兩部分，材料阻尼比以及流體介質(zhì)阻尼比。Lazan[6]通過大量實驗，總結(jié)了材料阻尼耗能與應(yīng)力幅值的關(guān)系。Rao[7]等對葉片的阻尼比進(jìn)行了討論分析，而且提出了一種計算等效粘性阻尼比的方法。李鵬飛等提出了材料阻尼比的識別方法與計算公式，并對Kelvin模型進(jìn)行了相應(yīng)的修正[8]。

吳俊男等研究了某離心葉輪的材料阻尼比與氣動阻尼比，并且給出了離心式葉輪的總阻尼比的非線性模型[9]。

葉輪在實際工作時，由于流場的不均勻，會受到流場施加的壓力脈動[10-11]。由于擾動，從而導(dǎo)致振動等一系列問題，影響到整體結(jié)構(gòu)的安全性。Chiang[12]研究了葉片受壓力脈動作用的響應(yīng)。Yao[13]等研究了雙吸泵的壓力脈動時域特性。Tanaka[14]等研究了高水頭泵的振動特性以及動應(yīng)力，同時Rao[15]等研究了汽輪機末級葉片的動應(yīng)力的計算。Xu[16]等研究了基于擬靜力的核主泵葉輪動應(yīng)力分析方法，分析了縮尺模型葉輪在壓力脈動荷載作用下的動力學(xué)響應(yīng)，但是沒有考慮壓力荷載的相位分布，只考慮了壓力荷載的大小，并且是通過擬靜力分析來求解動響應(yīng)，因此動力學(xué)分析存在一定的誤差。

本文首先對真機葉輪進(jìn)行了模態(tài)分析，得到葉輪濕模態(tài)參數(shù)，通過扇區(qū)模型與整體模型的對比，驗證了整體模型計算可靠性。并且比較了干濕模態(tài)固有頻率。通過流固耦合分析，計算了流體介質(zhì)阻尼耗功，得到了流體介質(zhì)阻尼比以及材料阻尼比，并且分析了流體介質(zhì)阻尼比與材料阻尼比隨著振幅變化關(guān)系。然后對葉輪壓力脈動進(jìn)行分析，對比了縮尺模型試驗與數(shù)值模擬水力性能，研究了壓力荷載時域、頻域信息，并且分析了壓力荷載幅值、相位的分布。最后通過葉輪干涉圖分析，找出共振區(qū)域，通過水中諧響應(yīng)分析，校核了葉輪共振，較為準(zhǔn)確的分析了葉輪動應(yīng)力。

1 CAP1400葉輪濕模態(tài)分析與阻尼比分析

葉輪的模態(tài)分為干模態(tài)與濕模態(tài)。干模態(tài)是指葉輪在空氣中的模態(tài)，濕模態(tài)是葉輪在水中模態(tài)。由于水的附加質(zhì)量的影響，濕模態(tài)特性與干模態(tài)特性相比，有顯著的差別，比如頻率等模態(tài)參數(shù)。而葉輪實際工作在水環(huán)境中，為了研究實際運行工況時的振動，需要考慮葉輪濕模態(tài)特性?？諝庵信c水環(huán)境中的固有頻率可以用式（1）、（2）表示。

式中,fa,fw分別表示在空氣中、水環(huán)境中的固有頻率；k表示模態(tài)剛度；M為模態(tài)質(zhì)量；Ma為附加質(zhì)量?？梢钥闯?，由于水的附加質(zhì)量的影響，在水中的模態(tài)頻率比在空氣中的低。濕模態(tài)的基本思想是以結(jié)構(gòu)動力學(xué)方程為求解對象，為了體現(xiàn)流體對結(jié)構(gòu)的影響，動力學(xué)計算時候，把流體的附加質(zhì)量加到結(jié)構(gòu)的質(zhì)量陣上。如式（3）所示。

式中,Ms為結(jié)構(gòu)的質(zhì)量矩陣；Ma為流體的附加質(zhì)量矩陣；C為阻尼矩陣；R?,R?以及R分別為結(jié)構(gòu)的加速度矢量、速度矢量以及位移矢量。

1.1 CAP1400葉輪濕模態(tài)分析

為了校核核主泵葉輪動力學(xué)響應(yīng)，首先需要對葉輪進(jìn)行濕模態(tài)分析。葉輪材料參數(shù)如表1所示。

表1 葉輪材料屬性Tab.1 Properties of the impeller material

葉輪濕模態(tài)計算模型如圖1所示。葉輪模型采用ANSYS中SOLID187單元、水域采用FLUID221單元。

圖1 數(shù)值模擬模型Fig.1 Numerical simulation model

對于循環(huán)對稱模型，為了保證結(jié)果的準(zhǔn)確性，常采用扇區(qū)模型進(jìn)行計算。葉輪濕模態(tài)需要考慮葉輪與流體之間的耦合作用，模態(tài)提取方法為非對稱提取，因此葉輪濕模態(tài)計算采用葉輪整體模型。為了驗證整體模型的準(zhǔn)確性，對比扇區(qū)與整體模型模態(tài)計算結(jié)果。圖2為扇區(qū)結(jié)構(gòu)與整體結(jié)構(gòu)干模態(tài)固有頻率對比結(jié)果。從圖2中可以看出，整體結(jié)構(gòu)與扇區(qū)結(jié)構(gòu)固有頻率結(jié)果差別不大。0節(jié)徑下，兩者誤差為0.37%，1節(jié)徑下，誤差為0.37%，2節(jié)徑下誤差為0.21%。因此，采用整體結(jié)構(gòu)計算葉輪模態(tài)是可行的。圖3為葉輪各階干濕模態(tài)頻率。從圖3中可以看出，每節(jié)徑下，葉輪濕模態(tài)頻率比干模態(tài)頻率要低，并且不同節(jié)徑，頻率的下降率不同。0節(jié)徑下，頻率下降率為11.1%，1節(jié)徑下，頻率下降率為3.8%，2節(jié)徑下，頻率下降率為16.1%。

圖2 扇區(qū)結(jié)構(gòu)與整體結(jié)構(gòu)固有頻率對比Fig.2 Comparison of natural frequency of sector structure and overall structure

圖3 葉輪干濕模態(tài)固有頻率Fig.3 Air and wet modal natural frequency of impeller

1.2 CAP1400葉輪阻尼比分析

葉輪阻尼比對葉輪動力學(xué)特性有一定影響，為了精確研究葉輪動力學(xué)特性，需要計算出葉輪阻尼比。核主泵葉輪的阻尼比分為材料阻尼比與流體介質(zhì)阻尼比。通過模態(tài)分析，根據(jù)公式（4）、（5）材料阻尼比可以計算出來。

式中,D0為材料阻尼耗功；σa為應(yīng)力幅值；σf為疲勞極限；U0為模態(tài)應(yīng)變能。

圖4為材料阻尼比隨著振幅的變化曲線，可以看出，材料阻尼比隨著振幅變化而變化，而且為非線性變化關(guān)系。

圖4 材料阻尼比隨著振幅變化曲線Fig.4 Material damping ratio with different amplitude

葉輪振動時，流體耗功，消耗了能量。這部分由于流體耗功產(chǎn)生的阻尼比為流體介質(zhì)阻尼比。流體介質(zhì)阻尼比計算公式為（6）、（7）、（8）。計算流體介質(zhì)阻尼比最關(guān)鍵的部分是計算流體介質(zhì)阻尼比耗功。而流體介質(zhì)耗功可以通過流體與固體的耦合計算得出。

式中,W為材料阻尼耗功；p為壓力;?為速度；?為單位法向量；A為面積；t為時間；σa為應(yīng)力幅值；E為彈性模量。

模態(tài)振型作為振動邊界條件，通過流固耦合計算，求出流體介質(zhì)阻尼耗功，通過公式（6）、（7）、（8）得出葉輪流體介質(zhì)阻尼比。圖5為流固耦合計算模型。

圖5 流固耦合計算模型Fig.5 Fluid-solid interaction calculation model

從圖6可以看出，同一模態(tài)下，流體介質(zhì)阻尼比基本不隨著振幅變化而變化。比較圖4與圖6可以看出，同一模態(tài)下，流體介質(zhì)阻尼比比材料阻尼比大一個數(shù)量級，因此計算總阻尼比時可以忽略材料阻尼比影響。

圖6 流體介質(zhì)阻尼比隨著振幅變化曲線Fig.6 Hydrodynamic damping ratio with different amplitude

2 CAP1400葉輪壓力荷載分析

葉輪在實際運行中，由于動靜干涉以及流場不均勻性會產(chǎn)生一定的壓力脈動。為了分析葉輪水中動應(yīng)力，需要對葉輪所受到的壓力脈動進(jìn)行分析。圖7為葉輪非定常計算流體域模型。首先分析縮尺葉輪模型非定常計算，并且與試驗進(jìn)行對比，表2為縮尺模型試驗與數(shù)值模擬水力性能對比，其中流量誤差為1.8%，揚程誤差為1.8%，效率誤差為2.8%，誤差均在要求范圍內(nèi)，驗證了計算平臺的可靠性。然后分析真機葉輪流場特性，表3為真機葉輪流場計算邊界條件。

圖7 葉輪CFD計算模型Fig.7 Impeller channel model for CFD calculation

表2 縮尺模型試驗與數(shù)值模擬水力性能對比Tab.2 The boundary condition of C flow field

表3 CAP1400葉輪流場計算邊界條件Tab.3 The boundary condition of CAP1 400 impeller flow field

由于非定常性以及動靜干涉等原因，葉輪內(nèi)部壓力分布較為不均勻。圖8為葉片前緣、尾緣點的壓力脈動在兩個周期內(nèi)的時域圖，圖9為葉片前緣、尾緣點的壓力脈動頻域圖。從圖8可以看出，相比于葉片前緣，葉片尾緣的壓力脈動的峰值較多。這是因為在葉片尾緣處，由于葉輪與導(dǎo)葉之間的動靜干涉，造成了一個高頻的脈動，并且頻率大小與導(dǎo)葉個數(shù)有關(guān)。在本文中，導(dǎo)葉數(shù)為13，因此這個高頻脈動的頻率為13倍轉(zhuǎn)頻。從圖9中可以很清晰看出，由于尾緣靠近導(dǎo)葉，在尾緣處壓力脈動頻域圖中會出現(xiàn)一個較高的13倍轉(zhuǎn)頻的壓力峰值。而葉片前緣由于與導(dǎo)葉相隔較遠(yuǎn)，13倍轉(zhuǎn)頻的壓力幅值相應(yīng)較小。

圖8 葉片前緣與尾緣壓力脈動時域圖Fig.8 Pressure fluctuation time-domain diagram of blade leading edge and blade trailing edge

圖9 葉片前緣(a)尾緣(b)壓力脈動頻域圖Fig.9 Pressure fluctuation frequency-domain diagram of blade leading edge(a)and trailing edge(b)

葉片壓力荷載在葉片不同位置的幅值與相位角均不相同，圖10為葉片壓力荷載在1節(jié)徑頻率下葉片壓力面的幅值與相位分布圖。圖11為葉片壓力荷載在1節(jié)徑頻率下葉片吸力面的幅值與相位分布圖。可以看出，葉片不同位置的壓力脈動幅值不同，葉片壓力面幅值大于吸力面幅值。同樣的，在葉片不同位置，壓力脈動相位也是不一樣的，存在著相位差。因此，為了準(zhǔn)確計算壓力荷載對葉輪動力學(xué)特性影響，需要考慮葉片不同位置壓力荷載的幅值與相位。

圖10 1節(jié)徑1模態(tài)固有頻率下葉片壓力面壓力幅值、相位分布Fig.10 Blade pressure amplitude(a)and phase distribution(b)of pressure side in 1-ND first mode frequency

圖11 1節(jié)徑1模態(tài)固有頻率下葉片吸力面壓力幅值(a)相位(b)分布Fig.11 Blade pressure amplitude(a)and phase(b)distribution of suction side in 1-ND first mode frequency

3 CAP1400葉輪諧響應(yīng)分析與共振校核

通過分析葉輪濕模態(tài)、阻尼比以及壓力荷載，得到求解葉輪動應(yīng)力時需要的輸入量。圖12為葉輪諧響應(yīng)分析流程圖。

圖12 諧響應(yīng)分析流程圖Fig.12 Harmonic response analysis process

核主泵葉輪工作在液體環(huán)境中，對葉輪進(jìn)行動力學(xué)分析（諧響應(yīng)分析）需要考慮周圍流體介質(zhì)的影響。為了校核葉輪強度，需要對葉輪進(jìn)行共振校核。圖13為葉輪干涉圖。從圖中可以看出，葉輪在1節(jié)徑1模態(tài)下固有頻率與激振力頻率接近，因此，需要進(jìn)行1節(jié)徑1模態(tài)共振分析。

圖13 葉輪干涉圖Fig.13 Interference diagram of impeller

葉輪所受到的荷載為流場中的壓力荷載，通過流體域節(jié)點與固體域節(jié)點的信息傳遞，將此前分析的壓力荷載幅值以及相位信息施加到固體相應(yīng)位置。圖14為葉片流體域節(jié)點與固體域節(jié)點分布圖。

對葉輪進(jìn)行諧響應(yīng)分析時，需要考慮周圍液體環(huán)境的影響。設(shè)置水域中葉輪表面為流固耦合交界面，同時設(shè)置自由表面，剛性壁面。圖15為流體域節(jié)點以及流固耦合交界面示意圖，紅色部分的流體域節(jié)點為流固耦合交界面，此處的單元是將流體與固體聯(lián)系在一起的單元，是結(jié)構(gòu)至流體的過渡。通過諧響應(yīng)分析，得到了1節(jié)徑1模態(tài)葉輪共振動力學(xué)響應(yīng)。

圖14 葉片流體域節(jié)點與固體域節(jié)點分布Fig.14 Blade fluid domain nodes and solid domain nodes distribution

圖15 流體域節(jié)點與流固耦合交界面Fig.15 Fluid domain nodes and fluid-solid interaction interfaces

圖16為葉輪位移分布圖與動應(yīng)力分布圖，可以看出，在受到流場中壓力脈動荷載作用下，葉輪1節(jié)徑1模態(tài)共振時，最大動應(yīng)力為2.02MPa，并且最大應(yīng)力點出現(xiàn)在葉根部分。

圖16 葉輪1節(jié)徑1模態(tài)共振等效應(yīng)力云圖Fig.16 Equivalent stress of impeller in 1-ND first mode resonance

4 結(jié)論

本文通過比較扇區(qū)模型與整體模型的模態(tài)頻率，驗證了整體模型的準(zhǔn)確性。通過模態(tài)分析，比較了干濕模態(tài)頻率。發(fā)現(xiàn)葉輪濕模態(tài)頻率比干模態(tài)頻率低，并且每階模態(tài)下，頻率下降率不同。通過模態(tài)分析，根據(jù)材料阻尼比計算公式計算出材料阻尼比。同時，通過流固耦合分析，計算出葉輪流體介質(zhì)阻尼比。對比材料阻尼比與流體介質(zhì)阻尼比發(fā)現(xiàn)，材料阻尼比比流體介質(zhì)阻尼比小一個數(shù)量級，可以忽略材料阻尼比的影響。

然后分析葉輪流場，首先對比了縮尺葉輪模型試驗與數(shù)值模擬水力性能，兩者揚程誤差為1.8%，效率誤差為2.8%。并研究了葉輪所受到的壓力荷載。分析葉片前緣與尾緣的壓力脈動時域圖、頻域圖發(fā)現(xiàn)，由于動靜干涉的影響，尾緣處13倍轉(zhuǎn)頻的壓力幅值較大。而且在葉片不同位置，壓力荷載的幅值與相位均不相同。在研究葉輪動應(yīng)力時，需考慮壓力荷載的幅值與相位。

最后通過分析葉輪干涉圖發(fā)現(xiàn)，葉輪1節(jié)徑1模態(tài)處固有頻率與激振力頻率接近，因此校核了1節(jié)徑1模態(tài)共振。通過流體節(jié)點與固體節(jié)點的耦合，將壓力荷載信息傳遞到固體對應(yīng)位置。通過水中葉輪諧響應(yīng)分析，得到了在1節(jié)徑1模態(tài)共振下的葉輪動應(yīng)力分布，研究發(fā)現(xiàn)，最大動應(yīng)力為2.02MPa，并且出現(xiàn)在葉根部分。

[1]孔銘.核主泵葉輪干濕模態(tài)的實驗探究[D].大連理工大學(xué),2014.

[2]張新,鄭源,錢鈞,等.基于流固耦合的臥式軸流泵葉輪模態(tài)分析[J].水電能源科學(xué),2015(7):164-167.

[3]Liang Q W,Rodríguez C G,Egusquiza E,et al.Modal Response of Hydraulic Turbine Runners[C]//International Association of Hydra Engineering&Research,Symposium on Hydraulic Machinery and Systems.2006.

[4]鄭小波,羅興琦,鄔海軍,等.軸流式葉片的流固耦合振動特性分析[J].西安理工大學(xué)學(xué)報,2005,21(4):342-346.

[5]張學(xué)榮.葉片在水體下的模態(tài)分析[J].排灌機械工程學(xué)報,2002,20(5):10-12.

[6]B.J.Lazan,Damping of Materials and Members in Structural Mechanics,Pergamon Press,New York,1968.

[7]Rao J S,Saldanha A.Turbomachine blade damping[J].Journal of Sound and Vibration,2003,262(3):731-738.

[8]李鵬飛.應(yīng)力相關(guān)阻尼模型及其在梁式橋動力分析中的應(yīng)用[D].北京交通大學(xué),2014.

[9]吳俊男.離心壓縮機半開式葉輪阻尼及動應(yīng)力計算研究[D].大連理工大學(xué),2015.

[10]Fortes-Patella R,Longatte,Kueny L.Numerical analysis of unsteady flow in a centrifugal pump[J].ASME Fluid Machinery,1995(222):41-46.

[11]Gonzalez J,Fernandez J.Blanco,E.Numerical simulation of the dynamic effects due to impeller-volute interaction in a centrifugal pump[J]-Journal of Fluids Engineering,2002,124(2):348-355.

[12]Chiang H W D,Kielb R E.An analysis system for blade forced response[J].Journal of Turbomachinery,1993,115(4):762-770.

[13]Zhifeng Yao,Fujun Wang.Experimental Investigation of Time-Frequency Characteristics of Pressure Fluctuations in a Double-Suction Centrifugal Pump[J].ASME J.Fluids Eng.,2011.10,Vol.133/101303-1.

[14]H.Tanaka.Vibration behaviour and dynamic stress of runners of very high head reversible pump-turbines,in:Proceedings of the 15th IAHR Symposium,Belgrade,1990.

[15]Rao J S,Peraiah K C,Uday K S.Estimation of Dynamic Stresses in Last Stage Steam Turbine Blades under Reverse Flow Conditions[J].Advances in Vibration Engineering,Journal of Vibration Institute of India,2009,8(1):71.

[16]Shengli Xu,Xi He,Tao Sun,Xiaofang Wang.Research of Damping and Dynamic Stress for Impeller of Reactor Coolant Pump,in:17th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery,Maui,Hawaii,2017.