孫軍帥, 李旭鵬, 孫汝雷, 溫濟(jì)銘, 李小暢, 田瑞峰

(1.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院, 黑龍江 哈爾濱 150001; 2.中國船舶及海洋工程設(shè)計研究院, 上海 200011; 3.哈爾濱工程大學(xué) 航天與建筑工程學(xué)院, 黑龍江 哈爾濱 150001)

核電站蒸汽發(fā)生器作為核電站關(guān)鍵設(shè)備,其安全性關(guān)乎核電站的安全運行。在長期運行中,蒸汽發(fā)生器換熱管極容易受流體作用而發(fā)生振動,當(dāng)振動劇烈時,換熱管會與相鄰管束、支承結(jié)構(gòu)等發(fā)生碰撞或摩擦,導(dǎo)致管束失效,造成巨大經(jīng)濟(jì)損失[1],有必要對換熱管的流致振動現(xiàn)象進(jìn)行深入研究。

大渦模擬方法(large eddy simulation, LES)在湍流計算方面具有優(yōu)勢,隨著計算機(jī)技術(shù)的迅猛發(fā)展,具有較好的使用前景[2]。Lourenco等[3]較早地采用了粒子圖像測速法(particle image velocimetry, PIV)技術(shù)展開了Re=3 900的圓柱繞流研究。文獻(xiàn)[4-5]采用PIV技術(shù)對尾流區(qū)流場以及剪切層不穩(wěn)定性等進(jìn)行了研究。由于流場的湍流特性受實驗條件影響較大,且通過實驗研究能夠捕獲的流場信息極為有限,因此大量學(xué)者通過數(shù)值手段展開相關(guān)研究。Beaudan等[6]較早地采用LES方法對Re=3 900的圓柱繞流問題進(jìn)行了數(shù)值模擬。Kravchenko等[7]基于B樣條方法的高精度LES方法,研究了Re=3 900圓柱繞流問題,通過與實驗結(jié)果[3]對比發(fā)現(xiàn),實驗可能存在環(huán)境干擾從而導(dǎo)致了剪切層較早轉(zhuǎn)捩。Parnaudeau等[5]除了進(jìn)行了實驗研究外,同樣采用LES方法展開了數(shù)值研究。戰(zhàn)慶亮等[8]、端木玉等[9]采用LES方法研究了Re=3 900時圓柱繞流流場特性,通過對比研究驗證了計算程序在湍流場計算的準(zhǔn)確性。一些關(guān)鍵流場參數(shù)方面仍存在一定的爭議,這可能是由于實驗技術(shù)差異以及數(shù)值計算模型不同導(dǎo)致的。

單根換熱管流致振動從流致振動機(jī)理來說,本質(zhì)上屬于單圓柱渦激振動的范疇。文獻(xiàn)[10-13]對圓柱渦激振動進(jìn)行了系統(tǒng)性的實驗研究。文獻(xiàn)[14-15]在研究低質(zhì)量阻尼參數(shù)的雙自由度圓柱渦激振動時發(fā)現(xiàn):當(dāng)質(zhì)量比m*=2.6時,橫向振動響應(yīng)幅值會產(chǎn)生“超級上端分支”,并出現(xiàn)了“2 T”形態(tài)的尾跡模式。Kang等[16]、Li等[17]基于RANS法建立了雙自由度圓柱渦激振動問題的二維計算模型,但由于二維模型忽視了流場在高雷諾數(shù)時具有高度三維化渦體結(jié)構(gòu),導(dǎo)致得到的最大橫向振幅和“超級上端分支”的范圍與實驗結(jié)果存在一定的偏差。Gsell等[18]和Li等[17]分別采用直接數(shù)值模擬(direct numerical simulation,DNS)方法和LES方法建立了雙自由度圓柱渦激振動問題的三維計算模型以模擬“超級上端分支”,但結(jié)果還是均低于實驗值;同時,由于計算工況較少,無法對該問題進(jìn)行較好描述。

綜上所述,在單圓柱渦激振動問題上,主要存在以下問題:1)在流場計算方面,CFD模型存在差異,導(dǎo)致關(guān)鍵流場參數(shù)存在著一定爭議;2)在振動響應(yīng)計算方面,不論建立的計算模型是二維還是三維的,計算結(jié)果均與Jauvtis &Williamson[15]的實驗值存在較大差異,體現(xiàn)在超級上端分支范圍、最大橫向振幅、運動軌跡等方面。

本文基于LES方法建立了精細(xì)化流場計算模型,以雷諾數(shù)為3 900的固定圓柱繞流為例,進(jìn)行了網(wǎng)格敏感性分析,而后通過不同求解方法的對比研究選取HHT-α方法求解結(jié)構(gòu)動力學(xué)方程,最終結(jié)合局部動網(wǎng)格變形技術(shù),建立了換熱管渦激振動流固耦合計算模型,并對Jauvtis &Williamson[15](簡稱為J-W)的經(jīng)典渦激振動實驗的振動響應(yīng)進(jìn)行了預(yù)測。

1 流固耦合計算模型

1.1 大渦模擬控制方程

LES方法最早是由氣象學(xué)家Smagorinsky[19]提出的,其主要思路是對N-S方程進(jìn)行濾波處理,將湍流中渦分為大尺度渦和小尺度渦2種,對大渦直接計算,對小渦進(jìn)行模擬。對于不可壓縮流動,濾波后的N-S方程為:

(1)

(2)

亞格子應(yīng)力模型有很多形式[19-21],本文采用了WALE模型[21],該模型通過考慮湍流壁面及動量傳遞效應(yīng)等,保證了對近壁面流域模擬的準(zhǔn)確性。同時,采用SIMPLEC方法求解壓力速度耦合方程,采用二階隱式方法對時間項進(jìn)行離散求解。

1.2 結(jié)構(gòu)動力學(xué)模型

將單根換熱管簡化為剛性圓柱,并采用雙自由度的質(zhì)量-阻尼-彈簧模型對換熱管渦激振動力學(xué)模型進(jìn)行描述,則流向和橫向的控制方程分別為:

(3)

(4)

對上述方程采用HHT-α方法[22]進(jìn)行求解,該方法比傳統(tǒng)的Runge-Kutta方法和Newmark-β方法具有更優(yōu)的算法穩(wěn)定性和數(shù)值耗散性能等。以式(4)為例,進(jìn)行方程推導(dǎo):

(1-η)kyt+Δt+ηkyt=(1-η)Fdt+Δt+ηFdt

(5)

(6)

(7)

最終得到關(guān)于yt+Δt的方程:

(8)

其中:

(9)

(10)

由式(8)可得:

(11)

將式(6)～(11)編寫成C語言程序,以實現(xiàn)單根換熱管渦激振動計算中對結(jié)構(gòu)動力學(xué)部分的求解。

1.3 幾何模型

在固定圓柱繞流計算和單根換熱管渦激振動計算中,坐標(biāo)原點位于換熱管中心,流向為x坐標(biāo),橫向為y坐標(biāo),展向為z坐標(biāo)。二者計算域的相對大小相同,但直徑D不同。換熱管渦激振動計算的計算域如圖1所示,寬為20D,長為30D?；谖墨I(xiàn)[4,7,9,23]對不同展向長度的圓柱繞流問題進(jìn)行的研究,選取展向長度為πD。圓柱距離入口10D,距離上下邊界10D。入口為均勻速度入口條件,上下邊界以及展向的前后邊界均為對稱邊界,出口為壓力出口,圓柱表面為無滑移壁面。

圖1 渦激振動計算的計算域示意Fig.1 Schematic diagram of the computational domain for vortex-induced vibration calculation

2 網(wǎng)格敏感性分析

表1 敏感性分析的網(wǎng)格模型設(shè)置參數(shù)

圖2 G3網(wǎng)格模型xoy截面示意Fig.2 Schematic diagram of the xoy section of G3 model

為了得到圓柱繞流問題的合理統(tǒng)計數(shù)據(jù),根據(jù)Franke等[24]的研究,統(tǒng)計時間應(yīng)大于40個漩渦脫落周期,因此本文的總計算時間約為100個漩渦脫落周期,選取穩(wěn)定后的70個漩渦脫落周期數(shù)據(jù)作為有效數(shù)據(jù)進(jìn)行統(tǒng)計。

表2 網(wǎng)格敏感性分析結(jié)果對比Table 2 Comparison of grid sensitivity analysis results

4套網(wǎng)格的圓柱表面時均壁面速度梯度系數(shù)與實驗數(shù)據(jù)[28]的對比如圖3所示。壁面速度梯度系數(shù)VG定義為VG=|ζ|/Re0.5,|ζ|為無量綱速度梯度的絕對值?？梢钥闯?4套網(wǎng)格的時均壁面速度梯度系數(shù)變化均呈現(xiàn)先增后減的趨勢,并在θ約為50°處取得極大值,其中分離點為曲線與橫軸交點,隨著網(wǎng)格數(shù)量的增加,分離點逐漸收斂,曲線也逐漸收斂于實驗數(shù)據(jù),參數(shù)θsep逐漸收斂。但需要注意的是,在θ為180°處,計算得到的時均壁面速度梯度系數(shù)約為0,明顯小于實驗數(shù)據(jù)[28]。由于圓柱繞流尾渦交替脫落,θ為180°處的瞬時壁面速度梯度系數(shù)應(yīng)呈現(xiàn)周期性變化,作時均處理后其值應(yīng)近似為0,驗證了計算的正確性以及數(shù)據(jù)統(tǒng)計處理的合理性。

圖3 不同網(wǎng)格模型時均壁面速度梯度系數(shù)分布Fig.3 Time-averaged wall velocity gradient coefficient distributions for different grid models

4套網(wǎng)格的圓柱尾流區(qū)不同截面時均流向速度分布與相關(guān)實驗數(shù)據(jù)[3,5,29]及數(shù)值結(jié)果[5]的對比如圖4所示。為了將不同截面時均流向速度分布進(jìn)行對比,對各個截面的數(shù)據(jù)進(jìn)行了偏移處理?？梢钥闯?時均流向速度分布由近尾流區(qū)(x/D<3)的“U”型向遠(yuǎn)尾流區(qū)(3

圖4 不同網(wǎng)格模型圓柱尾流區(qū)不同截面時均流向速度分布Fig.4 Time-averaged streamwise velocity distributions of the centerline beyond cylinder for different grid models

圖5給出了G3網(wǎng)格模型的圓柱繞流瞬時流場渦量,圖中的渦量等值面采用Q準(zhǔn)則[30]表示,用無量綱流向速度u/U0進(jìn)行染色。可以看出,Re=3 900的圓柱繞流問題具有豐富的流動結(jié)構(gòu),顯示出復(fù)雜的漩渦脫落現(xiàn)象。同時,附著在圓柱壁面的剪切層呈現(xiàn)二維結(jié)構(gòu),而分離區(qū)和尾流區(qū)呈現(xiàn)出高度三維化的渦體結(jié)構(gòu),具有明顯的湍流特征。

圖5 圓柱繞流瞬時三維渦量等值面圖(Q=10)Fig.5 Isosurface diagram of instantaneous 3D vorticity of flow around a cylinder (Q=10)

綜上所述,對圓柱繞流問題而言,近尾流區(qū)流動情況復(fù)雜,特別是回流區(qū)對網(wǎng)格的要求較高,粗網(wǎng)格會導(dǎo)致剪切層較早分離,回流區(qū)長度變短,速度型提前由“U”型轉(zhuǎn)變?yōu)椤癡”型;而遠(yuǎn)尾流區(qū)對網(wǎng)格要求較低,不同網(wǎng)格的速度型與實驗值吻合較好,僅在峰值處存在誤差。也就是說,近尾流區(qū)流場對網(wǎng)格更敏感,因此需要對近尾流區(qū)的網(wǎng)格進(jìn)行加密,才能準(zhǔn)確地模擬圓柱繞流現(xiàn)象。因而,認(rèn)為G3可以準(zhǔn)確地模擬圓柱繞流現(xiàn)象,并且滿足網(wǎng)格無關(guān)性要求,則后續(xù)工作將基于G3網(wǎng)格模型開展。

3 換熱管渦激振動計算

3.1 結(jié)構(gòu)動力學(xué)方程求解方法對比研究

表3 不同求解方法的對比Table 3 Comparison of different methods

振幅響應(yīng)反映的是算法的耗散性能,即算法導(dǎo)致的能量衰減程度。Kane等[31]的研究表明,隨著計算時間的延長,Newmark-β類方法的能量衰減要明顯小于四階Runge-Kutta方法。本質(zhì)上HHT-α方法是由Newmark-β方法衍生而來,具有更優(yōu)的數(shù)值性能。這直接導(dǎo)致方法A3計算得到的振幅響應(yīng)更接近實驗值,而方法A2次之,方法A1最次。

頻率響應(yīng)反映的是算法的彌散性能,即算法導(dǎo)致的周期延長程度。選取了12個振動周期的橫向位移y/D的時程變化曲線,起始時刻相位相同,3種計算方法的對比如圖6所示?？梢钥闯?隨著計算時間的延長,誤差逐漸積累導(dǎo)致方法A1和方法A2的振動周期延長,與方法A3相比,方法A1計算的振動周期延長了6.60%,方法A2計算的振動周期延長了23.2%。也就是說,在長時間進(jìn)行計算時,方法A1和方法A2將會帶來較大的計算誤差,最終導(dǎo)致計算結(jié)果會偏離物理事實。因此,最終選取方法A3(HHT-α方法)對換熱管渦激振動計算中結(jié)構(gòu)動力學(xué)部分進(jìn)行求解。

圖6 不同計算方法的橫向位移y/D的時程變化曲線Fig.6 Time-history curves of transverse displacement y/D for different methods

3.2 單根換熱管渦激振動計算

基于建立的換熱管流致振動流固耦合計算模型對單根換熱管渦激振動問題進(jìn)行計算,計算參數(shù)與J-W的實驗參數(shù)保持一致。圓柱直徑D=0.038 1 m,質(zhì)量比m*=m/md=2.6,質(zhì)量阻尼參數(shù)(m*+ca)ζ=0.013,水中固有頻率fn=0.4 Hz,工質(zhì)為水,根據(jù)實驗數(shù)據(jù)推算出流體密度ρ=995.4 kg/m3,動力粘度μ=0.000 781 Pa·s。時間步長Δt取為0.005 s,保證每個結(jié)構(gòu)固有周期包括大約500個時間步。計算工況為2

考慮到計算中圓柱周圍網(wǎng)格需要吸收圓柱邊界的運動,采用局部動網(wǎng)格變形技術(shù)處理網(wǎng)格變形。對圓柱周圍直徑為8D范圍內(nèi)采用O形網(wǎng)格進(jìn)行網(wǎng)格加密,將其內(nèi)側(cè)直徑為3D范圍的網(wǎng)格設(shè)置為剛體網(wǎng)格區(qū)域(A區(qū)域),剩余部分的網(wǎng)格設(shè)置為變形網(wǎng)格區(qū)域(B區(qū)域);加密區(qū)以外的網(wǎng)格設(shè)置為靜止網(wǎng)格區(qū)域(C區(qū)域)。變形前后的網(wǎng)格如圖7所示,可以看出,采用局部動網(wǎng)格變形技術(shù)變形后的網(wǎng)格具有較好的網(wǎng)格拓?fù)浣Y(jié)構(gòu),網(wǎng)格質(zhì)量較高,不會出現(xiàn)“負(fù)體積”現(xiàn)象,可以保證計算的有效進(jìn)行。同時,由于只有變形網(wǎng)格區(qū)域的網(wǎng)格參與變形計算,可以大大節(jié)省計算資源。

圖7 xoy截面變形前后的網(wǎng)格Fig.7 Grid of the xoy section before and after deformation

3.2.1 振動響應(yīng)

將不同折算流速下渦激振動響應(yīng)(橫向和流向振幅響應(yīng)、橫向振動頻率響應(yīng))與J-W[15]的實驗數(shù)據(jù)以及Li[17]和Gsell[18]的計算結(jié)果進(jìn)行對比,如圖8所示?？梢钥闯?計算得到的橫向振幅和流向振幅與實驗數(shù)據(jù)吻合較好,并且計算結(jié)果要優(yōu)于Li和Gsell的計算結(jié)果,計算可以成功捕捉到實驗中的初始分支、超級上端分支以及下端分支,各個分支的折算流速范圍與實驗吻合較好,其中超級上端分支最大橫向振幅為1.48D,略小于J-W的實驗值1.50D,對應(yīng)的折算流速Ur為8.35,與實驗值更為接近。橫向振動頻率響應(yīng)也存在著明顯的3個分支,分別與初始分支、超級上端分支和下端分支相對應(yīng),且與實驗數(shù)據(jù)吻合較好。3個分支間以“跳躍”的方式過渡,這種過渡方式與實驗變化趨勢相同。需要說明的是,在某些折算流速下,橫向振動可能存在多頻現(xiàn)象,這里只選取了橫向振動主頻。

圖8 不同折算流速下振動響應(yīng)對比Fig.8 Comparison of vibration response at different reduced velocities

當(dāng)折算流速較小時,流向振幅和橫向振幅相對較小,但在2

圖9 Ur=3時流向和橫向振幅頻譜圖Fig.9 Spectrogram of streamwise and transverse amplitude at Ur=3

隨著折算流速的增大,橫向振幅逐漸增大,但流向振幅逐漸減小。當(dāng)Ur>4時,振動逐漸進(jìn)入鎖定區(qū)間,流向振幅和橫向振幅逐漸增大,橫向振動頻率由漩渦脫落頻率逐漸趨近于圓柱水中固有頻率。在Ur=5時,橫向振幅明顯增大,振動進(jìn)入超上端分支。在超上端分支區(qū)間5

當(dāng)Ur>8.65時,橫向振幅和流向振幅迅速減小,振動進(jìn)入下端分支,橫向振動頻率跳躍至第3頻率分支,該頻率分支的頻率與圓柱在空氣中固有頻率相近。在整個下端分支區(qū)間8.65

3.2.2 運動軌跡

本文對不同折算流速下渦激振動的運動軌跡進(jìn)行研究,將計算得到的運動軌跡與實驗數(shù)據(jù)進(jìn)行對比,如圖10所示。運動軌跡左為實驗值[15],右為本文計算值。對計算得到的運動軌跡進(jìn)行縮放處理,但整體形狀不發(fā)生改變。可以看出,計算得到的運動軌跡形狀以及流向位移與橫向位移的相位差φxy與實驗數(shù)據(jù)吻合較好。

圖10 運動軌跡以及流向與橫向位移的相位差對比Fig.10 Comparison of motion trajectory and phase difference between streamwise and transverse displacement

當(dāng)振動處于初始分支和超級上端分支時,運動軌跡為順時針方向;而當(dāng)振動進(jìn)入下端分支時,運動軌跡轉(zhuǎn)變?yōu)槟鏁r針方向。在Ur=2時,運動軌跡類似于一條橫線,這是因為此時橫向位移較小而流向位移較大。隨著折算速度的增加,橫向振幅和流向振幅均逐漸增大,在Ur=3時,流向振動處于鎖定狀態(tài),二者振幅相當(dāng),并且流向振動頻率約為橫向振幅的2倍,運動軌跡表現(xiàn)為標(biāo)準(zhǔn)“8”字形。當(dāng)振動進(jìn)入超級上端分支區(qū)間,橫向振幅和流向振幅繼續(xù)增大,但橫向振幅的漲幅明顯大于流向振幅,振動軌跡形狀逐漸由瘦長“8”字形向傾斜“8”字形轉(zhuǎn)變。當(dāng)Ur=7.5時,流向振幅達(dá)到最大值,隨著折算流速的繼續(xù)增加,流向振幅開始逐漸減小,而橫向振幅卻繼續(xù)增大,振動軌跡形狀持續(xù)變化。當(dāng)Ur=8.35時,運動軌跡轉(zhuǎn)變?yōu)樾略滦?。?dāng)振動進(jìn)入下端分支時,橫向振幅和流向振幅均驟減,但橫向振幅要大于流向振幅,運動軌跡又逐漸變?yōu)椤?”字形,但此時運動軌跡為逆時針方向。

3.2.3 尾跡模式

選取xoy截面(z=0),對渦激振動尾跡模式進(jìn)行研究,圖11給出了不同折算流速下的尾跡模式,其中Ωz為z方向渦量,定義為Ωz=?v/?x-?u/?y。

圖11 不同折算速度下尾跡模式對比Fig.11 Comparison of wake mode at different reduced velocities

可以看出,渦激振動尾跡模式分為“SS”模式、“2S”模式、“2T”模式和“2P”模式。當(dāng)Ur=2時,由于橫向振幅較小,而流向振動占主導(dǎo)地位,振動處于“SS”分支,每個振動周期內(nèi)圓柱后側(cè)脫落一對對稱的反向漩渦,與Jauvtis等[14]進(jìn)行的可視化實驗結(jié)果相同。隨著折算流速的增加,尾跡模式發(fā)生轉(zhuǎn)變,當(dāng)Ur=3時,每個振動周期內(nèi)圓柱兩側(cè)交替脫落方向相反的2個漩渦,該泄渦模式被稱為“2S”模式。

當(dāng)折算流速繼續(xù)增大時,振動逐漸進(jìn)入超級上端分支,在這過渡區(qū)間上,尾跡模式由“2S”模式向一種新的泄渦模式“2T”模式轉(zhuǎn)變,即每半個振動周期內(nèi)圓柱一側(cè)脫落3個漩渦,其中一個漩渦的強(qiáng)度遠(yuǎn)小于另外2個漩渦的強(qiáng)度,且2個強(qiáng)度較大的漩渦方向相同,如圖11(c)所示。當(dāng)振動完全處于超級上端分支時,尾跡模式完全轉(zhuǎn)變?yōu)椤?T”模式,如圖11(d)和(e)所示。當(dāng)振動進(jìn)入下端分支時,尾跡模式發(fā)生轉(zhuǎn)變,在每半個振動周期內(nèi)圓柱一側(cè)脫落一對方向相反強(qiáng)度相近的漩渦,稱為“2P”模式。

4 結(jié)論

1)在進(jìn)行固定圓柱繞流問題計算時,圓柱近尾流區(qū)(x/D<3)流動參數(shù)對網(wǎng)格較為敏感,粗網(wǎng)格會使剪切層提前轉(zhuǎn)捩,導(dǎo)致回流區(qū)長度變短,而圓柱遠(yuǎn)尾流區(qū)(3

2)與傳統(tǒng)的Runge-Kutta方法和Newmark-β方法相比,本文采用的HHT-α方法在求解結(jié)構(gòu)動力學(xué)方程上更具優(yōu)勢,體現(xiàn)在較低的數(shù)值耗散、較小的數(shù)值彌散以及較高的數(shù)值精度,使得計算得到的振幅和振動周期更貼近實驗數(shù)據(jù)。

3)在渦激振動計算中,與其他學(xué)者的結(jié)果相比,本文建立的渦激振動計算模型可以準(zhǔn)確預(yù)測雙自由度圓柱渦激振動行為,計算可以成功捕捉到初始分支、超級上端分支和下端分支,以及相應(yīng)的“SS”、“2S”、“2T”以及“2P”尾跡模式,超級上端分支范圍和最大橫向振幅以及不同分支間的過渡方式與實驗數(shù)據(jù)一致。

由于換熱管主要以管束的形式存在于蒸汽發(fā)生器中,在今后的工作中,可以將本文建立的精細(xì)化換熱管渦激振動流固耦合模型應(yīng)用于管束結(jié)構(gòu)流致振動計算中,以支持蒸汽發(fā)生器管束結(jié)構(gòu)的設(shè)計與優(yōu)化工作等。