馮志鵬，臧峰剛，張毅雄

(中國核動力研究設計院 核反應堆系統(tǒng)設計技術重點實驗室，四川 成都 610041)

在核反應堆結構中，圓柱和圓柱束結構是流體誘發(fā)振動的主要部件，例如在橫流作用下的蒸汽發(fā)生器傳熱管、燃料棒等。根據(jù)國內(nèi)外核電站中對蒸汽發(fā)生器部件失效的多例事故的統(tǒng)計分析可知，傳熱管的流致振動和與其相關的沖擊磨損、微動接觸疲勞加上介質(zhì)腐蝕使得管壁逐漸變薄，是導致傳熱管承壓能力降低而破裂的主要原因。因此，對其流體誘導振動特性進行深入研究是十分必要的。

傳熱管的流體誘導振動問題得到了許多學者的關注[1]，他們在很大雷諾數(shù)范圍內(nèi)對圓柱體的流致振動進行了一系列的研究，包括用有限元方法[2]或計算流體力學(CFD)方法[3]。Gabbai等[4]對過去二十幾年中，經(jīng)受穩(wěn)態(tài)均勻流的大多數(shù)圓柱結構的流致振動所取得的研究進展進行了全面綜述，詳細討論了目前流體誘導振動問題研究現(xiàn)狀的要點與不足。隨著工業(yè)設備的過程參數(shù)越來越高，流體與結構之間存在著強烈的交互作用[5]，因此，需要構建更精確的物理模型來分析結構和流體的相互作用。而現(xiàn)有的研究對象大多是二維彈性支撐剛性柱體，不能反映真實的三維流動狀態(tài)，且將彈性管簡化成彈性支撐的剛性管，忽略了管自身的彈性，不能考慮結構的彈性變形與流體流動的相互影響。

本文基于有限體積法和有限元法，利用動網(wǎng)格控制技術模擬結構的運動，聯(lián)合CFD和計算結構動力學(CSD)建立傳熱管的流固耦合振動模型，旨在通過數(shù)值方法來研究傳熱管和流體間的相互作用，為設計者提供可靠、合理的振動分析，為自主設計蒸汽發(fā)生器傳熱管束的流致振動分析提供技術支持。

1 數(shù)值模型

對于流體誘導振動問題，任何數(shù)值模擬都需考慮以下4個基本問題：流場模擬、結構振動模擬、流體-結構交互作用模擬、數(shù)據(jù)處理。對這樣一個系統(tǒng)來說，控制方程包括流體的質(zhì)量和動量守恒及描述結構運動的動力學方程。因此，流體誘發(fā)振動問題的研究涉及到流體控制方程與結構動力的雙向耦合，近年來已有許多關于流體誘發(fā)振動的綜述性文獻[6]報道。

本文采用有限體積法離散三維、黏性、瞬態(tài)、不可壓縮Navier-Stokes方程，并聯(lián)合大渦模擬方法求解流場區(qū)域，利用有限元方法離散傳熱管結構體，采用Newmark積分方法求解瞬態(tài)動力平衡方程來獲得結構的位移、速度等響應?？紤]結構大變形以及由大變形帶來的流場網(wǎng)格的變形問題，采用基于擴散光順的Diffusion方法控制運動邊界的網(wǎng)格更新，通過流固耦合交界面進行固體域和流體域間的數(shù)據(jù)傳遞，文獻[7]對解耦和求解過程進行了詳細說明，本文不再贅敘。

2 物理模型及計算參數(shù)

傳熱管的物理參數(shù)：管長0.5 m，外徑D=0.01 m，內(nèi)徑Di=0.009 5 m，彈性模量E=10 GPa，泊松比υ=0.3，密度ρs=6 500 kg/m3，阻尼比ζ=0.047。

流體參數(shù)：流體為水，密度ρ=998.2 kg/m3，動力黏度μ=1.003 mPa·s，無量綱入口流速Ur=U/fnD=0.5～17，其中U為來流速度(m/s)，fn為管的固有頻率(Hz)。

流場區(qū)域及邊界條件：流場計算區(qū)域及網(wǎng)格如圖1所示。選用ICEM CFD作為網(wǎng)格劃分工具，采用六面體結構化網(wǎng)格。圖1中左邊入口采用速度入口邊界條件，右端出口采用壓力出口邊界條件，其他外邊界按對稱邊界和固定壁面處理；管壁面為流固交界面，設為動網(wǎng)格邊界。

圖1 流場區(qū)域及局部網(wǎng)格示意圖

時間參數(shù)：結構計算與流體計算的時間步長取為0.000 25 s。

3 結果分析

3.1 響應分支

圖2 振幅隨Ur的變化

圖3 Griffin圖

圖4 升力、阻力、振幅隨Ur的變化

在下分支階段Ⅳ(5

在非相干階段Ⅴ(Ur≥9)，流動速度的進一步增加導致Ay/D與ClRMS的下落，兩者均趨于零，由此也可看出漩渦強度是自限定的；阻力系數(shù)維持在一個比較恒定的值。

從圖4還可看到，最大的阻力系數(shù)和升力系數(shù)并不是同時出現(xiàn)的，升力峰值先于阻力峰值。ClRMS曲線的顯著特征是存在1個尖峰，其恰好出現(xiàn)在初始分支與擬上端分支的轉換處，在轉換前，升力急劇增加，而當響應模式發(fā)生轉換后，升力迅速下落，因此就導致了這樣1個尖峰。

3.2 頻率鎖定

圖5示出傳熱管的旋渦脫落頻率fvs和振動響應頻率fex與固有頻率fn之比fvs/fn、fex/fn隨Ur的變化情況。圖5b中對角線代表靜止管的旋渦脫落頻率，即St=0.2，水平線代表固有頻率，即fex=fn，實點指功率譜中的主頻峰值。

圖5 旋渦脫落頻率和響應頻率隨Ur的變化

響應頻率fex隨Ur的變化規(guī)律為：在初始分支，升力與橫向位移的頻率與靜止管的旋渦脫落頻率fst相同，遵循St=0.2的線性關系；當進入擬上端分支和下端分支階段時，旋渦脫落被升力方向的運動所控制，從圖5可看出，旋渦脫落頻率與振動響應頻率偏離了St=0.2這條線，基本上等于管的固有頻率，即fex/fn≈1.0，即發(fā)生了“頻率鎖定”，同步區(qū)域為3

3.3 相位差

傳熱管的流致振動特性還可通過升力與位移間的相位角來表征有關響應信息，相位角φ定義為升力相對于橫向位移的相位差，φ與時間的函數(shù)可通過Hilbert變換得到。Lissajou圖表示了傳熱管運動與流體間的能量傳遞，其傾角則給出了位移與升力間相位角的估計。

圖6 相位角φ隨Ur的變化

圖6為相位角φ隨Ur的變化情況?？擅黠@看到，當Ur=9～13時，相位差經(jīng)歷了從同相到反相的改變，相位差的這種跳躍稱為“相位開關”，是一種典型的非線性現(xiàn)象。圖7示出了各Ur下的相位角時程，通過綜合分析各響應階段的相位角時程，發(fā)現(xiàn)當Ur<9時，相位角時程的變化趨勢類似，相位角φ接近于0°，其時程也非常恒定(如Ur=1.75)；在相位差的轉變階段(如Ur=10)，相位角φ的時程開始變得紊亂，也不再接近于0°；而當流速超出“頻率鎖定”的臨界速度，處于高Ur時，相位角φ在360°范圍內(nèi)變化，時程變得非常紊亂(如Ur=16)。

圖7 不同Ur下的相位差時程

圖8為不同Ur下管的運動軌跡，圖9為不同Ur下傳熱管的Lissajou圖。從圖8、9可看出，當Ur≤3時(圖9a)，相位差接近于0，升力與位移同相；隨著流速的增大，平衡位置向下游移動，橫向變形逐漸增大，管的運動軌跡是阻力方向和升力方向頻率比為2的“8”字形圖，當Ur=6時達到最大，相應的流體力也明顯增大，Lissajou圖為類似橢圓的平行四邊形(圖9b)；在非相干階段(Ur≥9)，運動軌跡變得紊亂，不再是“8”字形，升力系數(shù)與位移間的相位也由同相變?yōu)榉聪?圖9c)，升力、振幅隨流速的增大而逐漸減小。

圖8 不同Ur下的運動軌跡

3.4 相圖

圖9 不同Ur下的Lissajou圖

圖10 不同Ur下位移的相圖

圖11 不同Ur下升力的相圖

圖12 不同Ur下位移的Poincare截面

圖13 不同Ur下升力的Poincare截面

當2

Lissajou圖、運動軌跡、相圖、Poincare截面映射都表現(xiàn)了三維彈性管流致振動的周期性，對所有流速范圍內(nèi)各種響應的分析發(fā)現(xiàn)，根據(jù)不同階段的響應分支特征，Lissajou圖、運動軌跡、相圖、Poincare截面映射均呈現(xiàn)不同的形式，聯(lián)合這4種圖與響應分支圖對三維彈性管的流致振動進行分析，便可較容易地確定其運動行為及振幅響應特性。

4 結論

本文基于雙向流固耦合方法，建立了三維橫向流體誘發(fā)傳熱管振動的數(shù)值模型，對不同流速下單管的流致振動特性進行了數(shù)值計算，得出以下主要結論：

1) 高m*ζ低m*流體-彈性管耦合系統(tǒng)的響應與高m*ζ高m*系統(tǒng)和低m*ζ低m*系統(tǒng)相比，會出現(xiàn)擬上端分支，振幅響應具有更寬的同步區(qū)域。

2) 當Ur≤2時，管的振幅很小，流致振動機理為湍流激勵；升力系數(shù)均方根隨Ur的增加而增加，阻力系數(shù)均方根隨Ur的增加先增大后減小，在Ur約為2時達到最小值；相位角φ接近于0°。

3) 當3

4) 當Ur≥9時，升力系數(shù)均方根趨于0，橫向振幅隨流動速度的增加而下落；當Ur=9～13時，相位差經(jīng)歷了從同相到反相的改變，相位角φ在360°范圍內(nèi)變化，其時程變得非常紊亂。

5) 在均勻湍流流動作用下，三維彈性管的升力與橫向位移在Ur=0.5～18的范圍內(nèi)未出現(xiàn)周期解的分叉。

參考文獻：

[1] 高李霞，蔣自龍，馬建中，等. 螺旋式傳熱管束流致振動試驗研究[J]. 原子能科學技術，2008，42(增刊)：468-471.

GAO Lixia, JIANG Zilong, MA Jianzhong, et al. Experimental study on flow induced vibration of heat transfer helical tube bundle[J]. Atomic Energy Science and Technology, 2008, 42(Suppl.): 468-471(in Chinese).

[2] ANAGNOSTOPOULOS P. Numerical study of the flow past a cylinder excited transversely to the incident stream, Part 1: Lock-in zone, hydrodynamic forces and wake geometry[J]. Journal of Fluids and Structures, 2000, 14(6): 819-851.

[3] PLACZEK A, SIGRIST J F. Numerical simulation of an oscillating cylinder in a cross-flow at low Reynolds number: Forced and free oscillations[J]. Computers & Fluids, 2009, 38: 80-100.

[4] GABBAI R D, BENAROYA H. An overview of modeling and experiments of vortex-induced vibration of circular cylinders[J]. Journal of Sound and Vibration, 2005, 282: 575-616.

[5] 馮志鵬，張毅雄，臧峰剛. 直管束流固耦合振動的數(shù)值模擬[J]. 應用數(shù)學和力學，2013，34(11):1 165-1 172.

FENG Zhipeng, ZHANG Yixiong, ZANG Fenggang. Numerical simulation of fluid-structure interaction for tube bundles[J]. Applied Mathematics and Mechanics, 2013, 34(11): 1 165-1 172(in Chinese).

[6] SARPKAYA T. A critical review of the intrinsic nature of vortex-induced vibrations[J]. Journal of Fluids and Structures, 2004, 19: 389-447.

[7] 馮志鵬，臧峰剛，張毅雄. 雙彈性管流固耦合振動的數(shù)值模擬[J]. 原子能科學技術，2014，48(8)：1 428-1 434.

FENG Zhipeng, ZANG Fenggang, ZHANG Yi-xiong. Numerical simulation of fluid structure interaction in two flexible tubes[J]. Atomic Energy Science and Technology, 2014, 48(8): 1 428-1 434(in Chinese).

[8] KHALAK A, WILLIAMSON C H K. Motions, forces and mode transitions in vortex-induced vibrations at low mass-damping[J]. Journal of Fluids and Structures, 1999, 13(7-8): 813-851.

[9] SKOP R A, BALASUBRAMANIAN S. A new twist on an old model for vortex-excited vibrations[J]. Journal of Fluids and Structures, 1997, 11: 395-412.

[10] LI Tian, ZHANG Jiye, ZHANG Weihua. Nonlinear characteristics of vortex-induced vibration at low Reynolds number[J]. Commun Nonlinear Sci Numer Simulat, 2011, 16: 2 753-2 771.