馮志鵬, 臧峰剛, 劉 帥, 黃 旋, 齊歡歡

(中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實(shí)驗室,成都 610213)

在核工程領(lǐng)域,流致振動仍然是影響蒸汽發(fā)生器性能的主要限制因素,流彈失穩(wěn)(流體彈性不穩(wěn)定性)是最具破壞性的流致振動機(jī)理[1]。尤其是2012年美國圣奧諾弗雷核電站(SONGS)因流致振動而導(dǎo)致蒸汽發(fā)生器失效后,人們又開始重點(diǎn)關(guān)注管束的流彈失穩(wěn)。

流彈失穩(wěn)(fluid-elastic instability, FEI)是由兩種流體-結(jié)構(gòu)相互作用機(jī)制引起的：①與相鄰管振動的流體耦合有關(guān),稱為流體剛度控制機(jī)制;②與管速度相關(guān)的流體力分量相關(guān),稱為流體阻尼控制機(jī)制。已有研究表明[2],阻尼機(jī)制在較低質(zhì)量阻尼參數(shù)下占主導(dǎo)地位,剛度機(jī)制在較高質(zhì)量阻尼參數(shù)下占主導(dǎo)地位。壓水堆中的蒸汽發(fā)生器管束屬于低質(zhì)量阻尼參數(shù)范圍,容易發(fā)生流體阻尼控制型流彈失穩(wěn)。

長期以來,管束流彈失穩(wěn)研究主要依賴于試驗,提出的許多半經(jīng)驗?zāi)Ｐ秃雎粤送牧鲝?qiáng)度、雷諾數(shù)等因素的影響,但在實(shí)際工程中,管束中的流動通常是黏性和湍流的,在沒有經(jīng)驗數(shù)據(jù)的情況下對流動進(jìn)行詳細(xì)描述的最佳途徑是采用計算流體力學(xué)(computational fluid dynamics, CFD),如De Pedro等[3]研究了阻尼控制型流彈失穩(wěn)的數(shù)值模型,Khalifa等[4]利用ANSYS-CFX建立了管束的CFD模型,預(yù)測了流彈失穩(wěn)的穩(wěn)定性,馮志鵬等[5]通過理論建模和CFD相結(jié)合的方式,發(fā)展了流彈失穩(wěn)行為的數(shù)值預(yù)測方法。

盡管自SONGS事故以來,對流彈失穩(wěn)的分析方法與建模工具等方面都進(jìn)行了改進(jìn),如考慮到了之前忽略的一些因素、采用更復(fù)雜的公式來解釋兩相流效應(yīng)和流向不穩(wěn)定的影響,但導(dǎo)致這種不穩(wěn)定現(xiàn)象的根本因素尚未確定,仍有待于進(jìn)一步研究。同時,在通過CFD研究流彈失穩(wěn)時,尚缺乏有效的確定主導(dǎo)激勵機(jī)理的手段和判定管束流彈系統(tǒng)穩(wěn)定性的指標(biāo)。

本文針對壓水堆蒸汽發(fā)生器中容易發(fā)生的阻尼控制型流彈失穩(wěn),基于開源CFD工具OpenFOAM建立流固耦合數(shù)值模型,從宏觀響應(yīng)特性、能量輸入與耗散、升力與位移間的相關(guān)性、升力的頻譜特性等層面研究流彈失穩(wěn)行為,提出更好反映管束流彈系統(tǒng)穩(wěn)定性和確定主導(dǎo)激勵機(jī)理的分析方法。

1 數(shù)值模型

1.1 流場模型

平行三角形管束及其計算域,如圖1所示。節(jié)徑比Pr=1.5。計算域包含了21行7列管,其中深色的管可在橫向自由振動(稱為運(yùn)動管),上下兩側(cè)為半管。進(jìn)口位于第一排管子的15D處(D為管子外徑),已有研究表明該入口長度是合適的。

圖1 計算域

湍流模擬采用兩方程的剪切應(yīng)力輸運(yùn)湍流模型(k-ωSST模型,以下簡稱SST模型),其在管束的流動預(yù)測中與試驗吻合較好[6]。管壁面附近采用邊界層,使y+≤1以保證黏性亞表層被合理求解。邊界層數(shù)量為15層,并使用1.01的擴(kuò)展系數(shù)從邊界層單元過渡到四面體單元,網(wǎng)格細(xì)節(jié)如圖2所示。

入口設(shè)置為均勻流動,入口湍流強(qiáng)度取為1%;出口設(shè)為inletOutlet條件,OpenFOAM中inletOutlet邊界條件允許流體流出或流入?yún)^(qū)域,因此流體不會被迫流向出口方向;計算域沿z向的兩個端面賦給empty邊界條件。除了第11行中心的管為運(yùn)動壁面外,流場上下側(cè)和其它管表面為固定無滑移壁面。根據(jù)平均間隙流速,雷諾數(shù)范圍為Re=2×103～4×104。

圖2 網(wǎng)格細(xì)節(jié)

1.2 運(yùn)動管的模擬

將運(yùn)動管視為單自由度系統(tǒng),運(yùn)動方程為

(1)

1.3 離散求解

流場模擬基于開源CFD工具OpenFOAM。本文所用的離散化方案和求解器誤差的詳細(xì)信息,如表1所示。時間導(dǎo)數(shù)為隱式一階Euler格式、對流項為線性迎風(fēng)格式、擴(kuò)散項為非正交修正的中心差分格式;動量方程的求解采用幾何代數(shù)多重網(wǎng)格求解器并配合Gauss-Seidel型平滑器;其它量的求解采用光滑求解器并配合symGaussSeidel型平滑器;動網(wǎng)格求解器為dynamicMotionSolverFvMesh,運(yùn)動邊界求解器為sixDoFRigidBodyMotion。

在瞬態(tài)求解過程中自動調(diào)整時間步長,保障最大局部Courant數(shù)保持在預(yù)設(shè)極限以下,在目前的模擬中,Courant數(shù)被限制在Co<0.8。非結(jié)構(gòu)三維網(wǎng)格上的Courant數(shù)計算如下

(2)

式中,Φf=Af(uf·nf)是垂直于f面的速度通量,f面是體積為V的單元的面Af。

每個模擬分為兩個步驟：采用穩(wěn)態(tài)計算求解流場,此時所有的管子保持靜止;以穩(wěn)態(tài)結(jié)果作為初始條件,進(jìn)行具有管子運(yùn)動的瞬態(tài)模擬。

2 模型驗證

為驗證本文建立的數(shù)值模型,首先開展網(wǎng)格收斂研究,然后將預(yù)測的臨界流速與已有試驗數(shù)據(jù)對比,驗證本文建立的數(shù)值模型。

表1 離散格式和求解器

網(wǎng)格收斂是CFD分析的一個基本步驟。本文在靜止繞流、強(qiáng)迫運(yùn)動兩個維度上進(jìn)行網(wǎng)格無關(guān)性分析,在建立管子周圍的邊界層網(wǎng)格時,采取了兩種方式：

方式一：僅在運(yùn)動管及其周圍的7根管子(稱為Cluster)周圍建立邊界層,即mesh1～mesh6,如圖3所示。

方式二：在所有管子周圍建立邊界層,即mesh21～mesh25。

圖3 僅Cluster管帶有邊界層的情況

網(wǎng)格情況如表2所示。其中size表示每個管周向網(wǎng)格的大小,NTotal表示網(wǎng)格單元數(shù)。

定義升力系數(shù)和阻力系數(shù)分別為

(3)

式中：Fl和Fd分別為作用在管子上的升力和阻力;ρ為流體密度;U0為入口流速;Ax和Ay分別為管在阻力方向和升力方向的投影面積。

表2 網(wǎng)格情況

2.1 靜止繞流下的網(wǎng)格驗證

首先基于靜止繞流進(jìn)行網(wǎng)格敏感性研究,典型的升力系數(shù)、阻力系數(shù)的時程,如圖4所示。實(shí)際計算得到的y+,如圖5所示。由圖5可知,計算得到的最大y+<1,滿足SST模型對網(wǎng)格的要求。流體力的統(tǒng)計分布,如圖6所示。與Mahon[7]得出的規(guī)律一致,即升力系數(shù)呈雙峰分布、阻力系數(shù)呈單峰分布。

圖4 升阻力系數(shù)時程

圖5 y+時程

將各網(wǎng)格的升力系數(shù)均方根、阻力系數(shù)均方根列于圖7,圖中的誤差線表示各網(wǎng)格相對于最精細(xì)網(wǎng)格(mesh6)的相對誤差,虛線表示僅Cluster管子具有邊界層,實(shí)線表示所有管子具有邊界層。由圖7可知,隨著網(wǎng)格越精細(xì),誤差越小,升阻力系數(shù)趨于收斂,mesh24的相對誤差小于5%?？傮w來說,網(wǎng)格對阻力系數(shù)的影響相對較小,最大相對誤差在10%以內(nèi),對升力系數(shù)的影響較大,最大相對誤差達(dá)到37%。將各網(wǎng)格相應(yīng)的所有管子具有邊界層與僅Cluster管子具有邊界層的結(jié)果進(jìn)行對比,如圖8所示。由圖8可知,所有管子帶邊界層的升阻力系數(shù)均小于僅Cluster帶邊界層的情況,兩兩之間的差異隨著網(wǎng)格尺寸的增大而增大,升力系數(shù)的差異遠(yuǎn)大于阻力系數(shù)的差異,且所有管子帶邊界層時,升力系數(shù)的統(tǒng)計分布特征更加合理,因此,在管束的流場計算需考慮所有管子的邊界層。

(a) 升力系數(shù)

(b) 阻力系數(shù)

圖7 流體力系數(shù)隨網(wǎng)格的變化情況

圖8 兩種邊界層網(wǎng)格的對比

2.2 強(qiáng)迫振動下的網(wǎng)格驗證

為進(jìn)一步驗證網(wǎng)格獨(dú)立性,使運(yùn)動管做強(qiáng)迫振動,按照振幅y0和振蕩頻率f指定運(yùn)動

y(t)=y0sin(2πft)

(4)

運(yùn)動管做強(qiáng)迫振動時的流體力系數(shù)隨網(wǎng)格的變化情況,如圖9所示。同樣地,誤差線表示各網(wǎng)格相對于最精細(xì)網(wǎng)格(mesh6)的相對誤差,虛線表示僅Cluster管子具有邊界層,實(shí)線表示所有管子具有邊界層?？梢悦黠@看到,在強(qiáng)迫振動情況下,全部管子具有邊界層時的誤差最小。另一方面,邊界層尤其對升力系數(shù)影響顯著,所有管子具有邊界層的升力系數(shù)比僅Cluster管子具有邊界層的顯著減小,且差異總體在20%以上。進(jìn)一步,通過表2可知,所有管子具有邊界層,網(wǎng)格總數(shù)的增加并不顯著。

圖9 強(qiáng)迫振動時,流體力系數(shù)隨網(wǎng)格的變化情況

綜合靜止繞流和強(qiáng)迫振動的結(jié)果,選取mesh24(330 000個單元)開展后續(xù)研究。

2.3 臨界流速預(yù)測值對比

進(jìn)一步,選取較低質(zhì)量阻尼參數(shù)和較高質(zhì)量阻尼參數(shù)兩種典型工況進(jìn)行研究(PMD=0.07和PMD=56.99),其中,PMD=mlδ/ρD2,δ=2πζ,ml為單位管長質(zhì)量,ρ為流體密度,ζ為阻尼比。入口流速用管子直徑和振動頻率無量綱化,Ur=Ug/fD,Ug為間隙流速,定義為Ug=U0Pr/(Pr-1),Pr為管束的節(jié)徑比,U0為入口流速。

采用本文建立的數(shù)值模型,使運(yùn)動管在流體作用下作自由振動,確定臨界流速,然后將預(yù)測結(jié)果繪入穩(wěn)定性圖,如圖10所示。并與文獻(xiàn)中相似管束的試驗結(jié)果[8]以及工程驗證性試驗結(jié)果進(jìn)行對比,可以看到預(yù)測結(jié)果與試驗結(jié)果在較大的質(zhì)量阻尼參數(shù)(mass damping parameter, MDP)范圍內(nèi)吻合良好。

圖10 臨界流速預(yù)測結(jié)果與試驗結(jié)果的對比

3 結(jié)果與討論

管的振動響應(yīng)根據(jù)不同的Ur可以衰減或放大,三種典型的管子響應(yīng)形式如圖11所示。由圖11可知,在第一種情況下,系統(tǒng)的凈阻尼為正,因此系統(tǒng)是穩(wěn)定的(見圖11(a));在第二種情況下,系統(tǒng)的凈阻尼為負(fù),是不穩(wěn)定的,此時對應(yīng)于阻尼控制的流彈失穩(wěn)(見圖11(b));第三種情況介于第一種情況和第二種情況,難以僅通過管的響應(yīng)來判定系統(tǒng)是否失穩(wěn)(見圖11(c))。

(a) Ur=15

(b) Ur=57

(c) Ur=45

3.1 宏觀響應(yīng)特性

振動位移、流體力等宏觀響應(yīng)是表征管束流彈失穩(wěn)行為最直觀、最常用的方法,如流速超過臨界值后,管的振幅會迅速增長,發(fā)生所謂的流彈失穩(wěn)。

升阻力系數(shù)、橫向位移隨速度的變化情況,如圖12所示。

(a) PMD=56.99

(b) PMD=0.07

由圖12可知,總體上,升阻力系數(shù)隨著流速增大而呈減小趨勢,橫向位移隨著流速增大而呈增加趨勢,橫向位移并非隨著流速變化呈單調(diào)關(guān)系,這與Weaver早期在水洞中測得的響應(yīng)形式類似[9]。對于PMD=56.99,在Ur=25～40的速度范圍內(nèi),出現(xiàn)了一個橫向位移的峰值,但通過檢查升力的頻譜及管的響應(yīng)頻率,發(fā)現(xiàn)其并不是漩渦脫落共振,而在Ur=40處,響應(yīng)曲線的斜率發(fā)生了顯著變化,相應(yīng)的升力系數(shù)曲線也有一個小幅的增加。對于較小的MDP(PMD=0.07),升力系數(shù)波動幅度較大,沒有明確的變化規(guī)律,橫向位移在Ur=1.1時顯著增加,在Ur=1.25處出現(xiàn)了一個小的尖峰,但此時管束并沒有發(fā)生流彈失穩(wěn),同樣也不是漩渦脫落共振,Ur>1.5后,位移呈持續(xù)增長趨勢。

由圖12可知,流體力系數(shù)與橫向位移之間明確的相關(guān)性,在橫向位移變化劇烈的速度處,未發(fā)現(xiàn)流體力系數(shù)的明顯分界點(diǎn)。難以通過振幅-速度曲線的斜率突變點(diǎn)來確定臨界流速。

3.2 升力與位移的相關(guān)性

為了進(jìn)一步定量分析升力與位移間的相關(guān)性,圖13列出了各Ur下的相關(guān)系數(shù)。在MDP較高的情況下(PMD=56.99),相關(guān)系數(shù)非常小,說明管子的振動位移不是由于升力直接導(dǎo)致的,說明管的響應(yīng)是系統(tǒng)的行為。而在MDP較低的情況下(PMD=0.07),升力與位移完全相關(guān),說明振動位移主要來自于升力的強(qiáng)迫振動,也反映出這些工況下系統(tǒng)沒有發(fā)生流彈失穩(wěn)。

圖13 升力與橫向位移間的相關(guān)系數(shù)

升力與橫向位移之間的相干性與相位差,如圖14所示。由圖14可知,在10 Hz(固有頻率)以下,升力與橫向位移呈顯著相干,相位差趨于0,說明兩者同相。在10 Hz附近,相位差發(fā)生了從同相到反相的跳躍。在10 Hz以上,不同的MDP顯示了不同的相干性與相位差,MDP較大時(PMD=56.99),相干系數(shù)逐漸減小到0.5以下,相位差在0～180°之間變化;MDP較小時(PMD=0.07),在約10～20 Hz的頻率范圍內(nèi)呈顯著相干,相位差約為180°,在約20～100 Hz的頻率范圍內(nèi),相干性先減小再增大,相位差逐漸減小,在約100 Hz以上,呈顯著相干且相位差又保持在0°附近。通過分析其它速度下的相干性和相位差,發(fā)現(xiàn)差別不大,從相干函數(shù)和相位差中,看不出系統(tǒng)發(fā)生失穩(wěn)的明顯臨界點(diǎn)。

(a) PMD=56.99,Ur=57

(b) PMD=0.07,Ur=3

3.3 振動能量

進(jìn)一步通過計算管束流彈系統(tǒng)的能量輸入與能量損失,以期得到一些有價值的見解。定義運(yùn)動管的振動能量為

(6)

振動能量的時程,如圖15所示。其中：W為系統(tǒng)的振動能量;WFf為所有流體力所做的功;WFd為所有耗散力所做的功?？梢钥吹?流速較高時,系統(tǒng)的振動能量為正,說明系統(tǒng)失穩(wěn);Ur=45時,系統(tǒng)的振動能量總體上為負(fù),但在某些時刻會變?yōu)檎?說明此時系統(tǒng)已經(jīng)具有失穩(wěn)趨勢。系統(tǒng)的振動能量隨Ur的變化情況,如圖16所示。從圖16可知,明顯的臨界點(diǎn)。與振幅-流速曲線或流體力-流速曲線(見圖12)相比,振動能量由于同時考慮了流體力與振動響應(yīng),可以更好地反映穩(wěn)定性。

(a) Ur=45

(b) Ur=57

(a) PMD=56.99

(b) PMD=0.07

3.4 升力的頻譜特性

升力系數(shù)的頻譜如圖17～圖18所示。當(dāng)流速和雷諾數(shù)較低時,升力系數(shù)的頻譜為窄帶,具有明顯的單峰特性(見圖17(a));當(dāng)流速和雷諾數(shù)較高時(水中Ur>1.5,基于入口流速的Re=4 500;空氣中Ur>25,基于入口流速的Re=5 134),升力系數(shù)的頻譜開始變?yōu)閷拵问?湍流激勵或其他激勵源是主要激勵機(jī)理。升力系數(shù)頻譜特性發(fā)生變化的雷諾數(shù)與單圓柱繞流時發(fā)生轉(zhuǎn)捩的雷諾數(shù)(Re=5×103)[10]相當(dāng),說明此時流場形式發(fā)生了變化,湍流激勵增強(qiáng)。

(a) Ur=0.5

(b) Ur=3

(a) Ur=45

(b) Ur=57

升力的頻譜隨流速的變化呈現(xiàn)一定規(guī)律,如圖18所示。對于流體介質(zhì)為空氣的情況(PMD=56.99),在管束發(fā)生失穩(wěn)前存在一個頻率區(qū)間,在這個區(qū)間內(nèi),升力系數(shù)頻譜呈現(xiàn)寬頻特性,具有兩個較明顯的峰,一個對應(yīng)管子的固有頻率,另一個更高的頻率對應(yīng)流體激勵頻率,隨著流速增大,較高頻率的峰值也逐漸增大,且逐漸超過管子固有頻率對應(yīng)的峰值;而當(dāng)系統(tǒng)發(fā)生了流彈失穩(wěn)后,升力系數(shù)的頻譜從比較寬的頻帶變到窄帶,具有明顯的單峰特性(見圖18(b))。與空氣中不同,水中(PMD=0.07)升力系數(shù)的寬頻特性主要出現(xiàn)在10 Hz以下(見圖17(b))。

3.5 主導(dǎo)激勵機(jī)理

響應(yīng)特性在Ur=25～40的速度范圍內(nèi)出現(xiàn)了一個橫向位移的峰值(見圖12(a)),在文中前面的分析中指出,其并不是漩渦脫落共振,同時結(jié)合升力頻譜的寬帶特性可以認(rèn)為是湍流激勵而非漩渦脫落和流彈失穩(wěn)占主導(dǎo)地位,為了進(jìn)一步驗證這種結(jié)論,對比分析了靜止繞流和流固耦合的結(jié)果,如圖19～圖21所示。其中FSI表示流固耦合,Static表示靜止繞流,靜止繞流時管子的響應(yīng)通過將靜止繞流的流體力作用于管的動力學(xué)方程,再采用4階Runge-Kutta法求解而得到。

(a) PMD=56.99

(b) PMD=0.07

從圖19可知,PMD=56.99時FSI和Static的流體力幾乎相同,說明運(yùn)動對流體力的影響較小,這也體現(xiàn)在升力與位移間非常小的相關(guān)系數(shù)上(見圖13);而PMD=0.07時,FSI的升力遠(yuǎn)大于Static的,說明流體-振動耦合對流體力有明顯的放大作用,這與彈性管發(fā)生渦致振動時振動對流體力的放大作用類似[11]。

橫向位移的比較(見圖20),對于PMD=56.99,當(dāng)Ur較小時,主要機(jī)理是湍流激勵,因此FSI與Static兩者相差不大;隨著Ur的增大(Ur>40),FSI的位移顯著增大,說明流彈失穩(wěn)逐漸起主導(dǎo)作用。對于PMD=0.07,盡管最大的橫向位移已經(jīng)達(dá)到了12%D,但FSI與Static的差別并不大,僅從位移的對比來看,仍然是湍流激勵占主導(dǎo)地位,從振動能量的對比更可以清晰地反映這些特征(見圖21)。通過對比靜止繞流下的強(qiáng)迫振動與流固耦合振動,同時結(jié)合升力的頻譜特性、系統(tǒng)的振動能量,可以更全面地反映出占主導(dǎo)地位的激勵機(jī)理。

(a) PMD=56.99

(b) PMD=0.07

(a) PMD=56.99

(b) PMD=0.07

4 結(jié) 論

本文基于開源CFD工具OpenFOAM,建立了平行三角形管束的阻尼控制型流彈失穩(wěn)預(yù)測模型,分析了影響阻尼控制型不穩(wěn)定機(jī)制的關(guān)鍵參數(shù)和流動響應(yīng),主要結(jié)論如下：

(1) 難以通過振幅-速度曲線的斜率突變點(diǎn)、升力與位移的相關(guān)性來確定臨界流速,而振動能量綜合考慮了流體力與振動響應(yīng),可以更好地反映管束流彈系統(tǒng)的穩(wěn)定性特征。

(2) 在MDP較高的情況下,升力與橫向位移間的相關(guān)系數(shù)非常小,說明管子的振動位移不是由升力直接導(dǎo)致,管的響應(yīng)為系統(tǒng)行為;而在質(zhì)量阻尼參數(shù)較低的情況下,升力與位移完全相關(guān),說明振動位移主要是升力的強(qiáng)迫振動。

(3) 當(dāng)雷諾數(shù)較低時,升力頻譜為窄帶,具有明顯的單峰特性;當(dāng)雷諾數(shù)較高時,升力頻譜變?yōu)閷拵问?主導(dǎo)激勵機(jī)理是湍流激勵;升力頻譜特性發(fā)生變化的雷諾數(shù)與單圓柱繞流時發(fā)生轉(zhuǎn)捩的雷諾數(shù)相當(dāng);當(dāng)系統(tǒng)發(fā)生了流彈失穩(wěn)后,升力頻譜從寬帶變到窄帶。

(4) 通過對比靜止繞流下的強(qiáng)迫振動與流固耦合振動,同時結(jié)合升力頻譜特性與系統(tǒng)的振動能量,可以更為全面地反映出系統(tǒng)中占主導(dǎo)地位的激勵機(jī)理。