管束結構的阻尼控制型流彈失穩特性研究

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

(中國核動力研究設計院 核反應堆系統設計技術重點實驗室,成都 610213)

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

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

長期以來,管束流彈失穩研究主要依賴于試驗,提出的許多半經驗模型忽略了湍流強度、雷諾數等因素的影響,但在實際工程中,管束中的流動通常是黏性和湍流的,在沒有經驗數據的情況下對流動進行詳細描述的最佳途徑是采用計算流體力學(computational fluid dynamics, CFD),如De Pedro等[3]研究了阻尼控制型流彈失穩的數值模型,Khalifa等[4]利用ANSYS-CFX建立了管束的CFD模型,預測了流彈失穩的穩定性,馮志鵬等[5]通過理論建模和CFD相結合的方式,發展了流彈失穩行為的數值預測方法。

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

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

1 數值模型

1.1 流場模型

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

圖1 計算域

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

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

圖2 網格細節

1.2 運動管的模擬

將運動管視為單自由度系統,運動方程為

(1)

1.3 離散求解

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

在瞬態求解過程中自動調整時間步長,保障最大局部Courant數保持在預設極限以下,在目前的模擬中,Courant數被限制在Co<0.8。非結構三維網格上的Courant數計算如下

(2)

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

每個模擬分為兩個步驟：采用穩態計算求解流場,此時所有的管子保持靜止;以穩態結果作為初始條件,進行具有管子運動的瞬態模擬。

2 模型驗證

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

表1 離散格式和求解器

網格收斂是CFD分析的一個基本步驟。本文在靜止繞流、強迫運動兩個維度上進行網格無關性分析,在建立管子周圍的邊界層網格時,采取了兩種方式：

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

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

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

網格情況如表2所示。其中size表示每個管周向網格的大小,NTotal表示網格單元數。

定義升力系數和阻力系數分別為

(3)

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

表2 網格情況

2.1 靜止繞流下的網格驗證

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

圖4 升阻力系數時程

圖5 y+時程

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

(a) 升力系數

(b) 阻力系數

圖7 流體力系數隨網格的變化情況

圖8 兩種邊界層網格的對比

2.2 強迫振動下的網格驗證

為進一步驗證網格獨立性,使運動管做強迫振動,按照振幅y0和振蕩頻率f指定運動

y(t)=y0sin(2πft)

(4)

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

圖9 強迫振動時,流體力系數隨網格的變化情況

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

2.3 臨界流速預測值對比

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

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

圖10 臨界流速預測結果與試驗結果的對比

3 結果與討論

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

(a) Ur=15

(b) Ur=57

(c) Ur=45

3.1 宏觀響應特性

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

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

(a) PMD=56.99

(b) PMD=0.07

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

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

3.2 升力與位移的相關性

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

圖13 升力與橫向位移間的相關系數

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

(a) PMD=56.99,Ur=57

(b) PMD=0.07,Ur=3

3.3 振動能量

進一步通過計算管束流彈系統的能量輸入與能量損失,以期得到一些有價值的見解。定義運動管的振動能量為

(6)

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

(a) Ur=45

(b) Ur=57

(a) PMD=56.99

(b) PMD=0.07

3.4 升力的頻譜特性

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

(a) Ur=0.5

(b) Ur=3

(a) Ur=45

(b) Ur=57

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

3.5 主導激勵機理

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

(a) PMD=56.99

(b) PMD=0.07

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

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

(a) PMD=56.99

(b) PMD=0.07

(a) PMD=56.99

(b) PMD=0.07

4 結 論

本文基于開源CFD工具OpenFOAM,建立了平行三角形管束的阻尼控制型流彈失穩預測模型,分析了影響阻尼控制型不穩定機制的關鍵參數和流動響應,主要結論如下：

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

(2) 在MDP較高的情況下,升力與橫向位移間的相關系數非常小,說明管子的振動位移不是由升力直接導致,管的響應為系統行為;而在質量阻尼參數較低的情況下,升力與位移完全相關,說明振動位移主要是升力的強迫振動。

(3) 當雷諾數較低時,升力頻譜為窄帶,具有明顯的單峰特性;當雷諾數較高時,升力頻譜變為寬帶形式,主導激勵機理是湍流激勵;升力頻譜特性發生變化的雷諾數與單圓柱繞流時發生轉捩的雷諾數相當;當系統發生了流彈失穩后,升力頻譜從寬帶變到窄帶。

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