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

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

蒸汽發(fā)生器是核島內(nèi)的關鍵設備之一,傳熱管束是其核心部件,容易發(fā)生流致振動問題。在已發(fā)現(xiàn)的流致振動機理中,流彈失穩(wěn)(流體彈性不穩(wěn)定性)能在短時間內(nèi)造成結(jié)構(gòu)損壞,在工程設計中必須避免發(fā)生這種現(xiàn)象。

目前,工程中采用基于準靜態(tài)理論得出的Connors公式計算臨界流速[1]

(1)

式中:K為Connors系數(shù);fi、ζi、D、m分別為傳熱管的第i階頻率、第i階模態(tài)阻尼比、外徑、單位長度質(zhì)量;ρ為流體密度;無量綱量2πζim/ρD2稱為質(zhì)量阻尼參數(shù)(mass damping parameters,MDP)。隨著對流彈失穩(wěn)研究的深入,人們發(fā)現(xiàn)Connors公式較難描述流彈失穩(wěn)機理,而其他理論模型,如準穩(wěn)態(tài)模型、非穩(wěn)態(tài)模型和流管模型,由于需要大量試驗參數(shù),制約了工程應用與推廣。

在沒有經(jīng)驗數(shù)據(jù)的情況下詳細描述流動的最佳途經(jīng)是采用計算流體力學(computational fluid dynamics,CFD),但考慮到真實管束結(jié)構(gòu)及其眾多可能的排列形狀,使得仿真計算量大[2],時間成本高,目前這種方法僅限于非常簡單幾何結(jié)構(gòu)的數(shù)值模擬。一種更實用的方法是利用CFD來預測可以集成到理論模型中的參數(shù),如Meskell等[3]利用數(shù)值模型估算了正三角形管束的穩(wěn)態(tài)流體力,然后將這些力參數(shù)與準穩(wěn)態(tài)模型結(jié)合使用;Hassan等[4]利用CFD來模擬非穩(wěn)態(tài)流體力和由管道運動產(chǎn)生的相應力系數(shù),并與非穩(wěn)態(tài)模型結(jié)合使用;趙燮霖等[5]采用CFD獲取半解析模型中的相位延遲函數(shù),預測了橫向流作用下間距比為1.375的平行三角形與正三角形管束結(jié)構(gòu)的流彈失穩(wěn)閾值。然而,目前還沒有一個全面的數(shù)值研究來預測這些流彈失穩(wěn)理論模型中的全套流體力相關參數(shù),也缺乏對這些理論模型的對比研究。

因此,本文采用理論建模和CFD計算相結(jié)合的方式,推導出準穩(wěn)態(tài)模型、非穩(wěn)態(tài)模型和流管模型的動力學方程,發(fā)展理論模型中眾多參數(shù)的辨識方法,開展管束結(jié)構(gòu)流彈失穩(wěn)的數(shù)值預測,最終確定臨界流速、對比分析各預測模型,為蒸汽發(fā)生器的流致振動分析與設計提供技術支撐。

1 動力學方程

壓水堆核電站采用立式自然循環(huán)式蒸汽發(fā)生器,直立式倒U形傳熱管束是其核心部件,圖 1為U形彎管段的結(jié)構(gòu)和流場示意圖。本文的研究對象為華龍一號的ZH65型蒸汽發(fā)生器的管束,主要幾何參數(shù)如表1所示。

考慮由N根管組成的管束,將管束中的管子看成多跨梁結(jié)構(gòu),利用梁的振動理論,第i根管在x向和y向的運動方程可寫成

(2)

式中:EI為彎曲剛度;z為管子軸向坐標;m為單位管長質(zhì)量;c為結(jié)構(gòu)阻尼系數(shù);g和h分別為x向和y向的流體彈性力;下標i為第i根管。

根據(jù)式(2)中流體彈性力的選取,將形成不同的流彈失穩(wěn)理論模型。首先通過理論推導,全面建立起目前研究最廣泛的三種流彈失穩(wěn)理論模型的控制方程和關鍵參數(shù)的數(shù)學模型,以開展流體力相關參數(shù)辨識和數(shù)值預測方法建立。

1.1 準穩(wěn)態(tài)模型

(3)

橫流作用下的管束的運動方程可寫為

(4)

(5)

式中:[I2]、[A]、[B]、[K]分別為單位矩陣、附加質(zhì)量矩陣、流體阻尼矩陣、流體剛度矩陣;t為時間;ω為圓頻率;δ為對數(shù)衰減率;a=U/U∞。

1.2 非穩(wěn)態(tài)模型

在非穩(wěn)態(tài)理論中[7],作用在第i根管上的非穩(wěn)態(tài)流體力取決于管子自身及相鄰管子在流場中的位移、速度和加速度,表示為

式中:下標i、j為第j根管運動對第i根管所受流體力的影響;α、σ、τ、β的各種形式為流體力相關系數(shù);R為管半徑;Ug為管束的間隙流速,其它變量與前文相同。

得到流體力后,與振動方程聯(lián)立,管束的動力學方程可表示為

(8)

式中:[Ms]、[Cs]、[Ks]為管的質(zhì)量矩陣、阻尼矩陣和剛度矩陣;[Ma]、[Ca]、[Ka]為流體的附加質(zhì)量矩陣、附加阻尼矩陣和附加剛度矩陣;z={x,y}′。

1.3 流管模型

流管模型基于一維非穩(wěn)態(tài)流動理論,用曲線坐標s表示流動路徑,瞬時流管面積、速度和壓力分解為平均項和波動項

A(s,t)=A0+a(s,t)

(9)

U(s,t)=U0+u(s,t),P(s,t)=P0+p(s,t)

(10)

采用相位滯后函數(shù)φ(s)來考慮擾動延遲效應,面積擾動函數(shù)a(s,t)表示為

a(s,t)=a(sm,t)f(s)eiφ(s)

(11)

式中:a(sm,t)為最小間隙位置的面積擾動;f(s)為人工衰減函數(shù)。

不可壓縮流體的連續(xù)性方程為

(12)

式中:U(s,t)為流體速度向量;n(s)為垂直于控制體表面的單位向量。

將式(9)和式(10)代入式(12),假設擾動是頻率為ω的諧波,從入口s=si到出口s=se沿曲線坐標s積分,消除穩(wěn)態(tài)項和高階項,并無量綱化

(13)

采用類似推導過程,可得到控制體的動量方程

(14)

(15)

得到流體力后,與振動方程聯(lián)立,即可建立基于流管模型的動力學方程

(16)

理論上,流管模型可以不需要任何試驗數(shù)據(jù),但在實際使用過程中,模型中許多參數(shù)來自于定性的流動可視化和假定,缺乏合理的定量表征方式,本文利用數(shù)值方法來確定流管模型中的兩個關鍵參數(shù):流管邊界和相位滯后。

2 流體力參數(shù)辨識

本文將理論模型與CFD計算相結(jié)合,通過仿真數(shù)據(jù)驅(qū)動的方式,發(fā)展計算這些參數(shù)的辨識方法,進而建立預測管束流彈失穩(wěn)行為的數(shù)值預測方法。

2.1 CFD模型

管束及計算域和局部網(wǎng)格示意圖,如圖4所示。其中,上下兩側(cè)為半管,出口在最后一行管的間隙處截斷;入口區(qū)域設置為速度入口(U∞),出口區(qū)域設為inletOutlet條件,振動管設置為運動壁面,流場上下兩側(cè)和其它管表面為固定無滑移壁面;管束周圍邊界層數(shù)量為15層,周向尺寸為0.01D,第一個網(wǎng)格距管束壁面為0.1%D,以保證y+≤1,使用1.01的擴展系數(shù)從邊界層單元過渡到四面體單元。流場模擬基于開源CFD工具OpenFOAM,湍流模擬采用SST模型。文獻[8]詳細介紹了數(shù)值模型,本文不再贅述。

為了定量描述數(shù)值預測結(jié)果與實驗結(jié)果的相對誤差情況,定義相對均方根誤差RMSE

(17)

式中:xi為計算值;yi為試驗值;N為樣本數(shù)目。

2.2 準穩(wěn)態(tài)模型的流體力參數(shù)辨識

2.3 非穩(wěn)態(tài)模型的流體力參數(shù)辨識

對于非穩(wěn)態(tài)模型,以圖4中1根管在升力和阻力方向上作強迫振蕩,其余管固定,來計算作用在各管上的流體力等信號。

計算的升力系數(shù)幅值和相位差與文獻中的試驗結(jié)果比較,如圖6所示[9]。計算的流體阻尼、流體剛度與文獻中的試驗結(jié)果比較[10],如圖7所示。從圖7可知,數(shù)值預測結(jié)果基本上反映了流體力系數(shù)的變化情況,預測結(jié)果與試驗數(shù)據(jù)吻合較好,升力系數(shù)幅值的相對均方根誤差為88%,在Ur=10處的誤差較大;相位差的相對均方根誤差為12.99%。采用相同的方式,可得到所有其它流體力參數(shù)。

2.4 流管模型的流體力參數(shù)辨識

聯(lián)合圖像處理技術和CFD流場計算,開發(fā)了流管特征提取程序,識別的流管邊界如圖8所示。然后在識別出的流管邊界中心線上布置監(jiān)測點,采集這些監(jiān)測點的速度時程,利用互相關函數(shù)法,計算出管運動和由此產(chǎn)生的流動擾動之間相位滯后。計算得到的相位滯后與試驗數(shù)據(jù)[11]的比較,如圖9所示。從圖9可知,本文建立的數(shù)值預測方法可以很好地預測相位滯后,相對均方根誤差為35%。

3 臨界流速預測

利用CFD計算可以辨識不同輸入?yún)?shù)的優(yōu)勢,在成功獲得準穩(wěn)態(tài)模型、非穩(wěn)態(tài)模型、流管模型中的全套流體力相關參數(shù)的基礎上,建立管束結(jié)構(gòu)流彈失穩(wěn)的數(shù)值預測方法,最終確定臨界流速,并對比研究三種理論模型的預測效果。

在工程中,流體密度、橫向流速沿傳熱管分布并不均勻,需要計算相應的等效值,等效流速定義為

(18)

等效質(zhì)量和等效密度定義為

(19)

(20)

根據(jù)工程經(jīng)驗,選取由4根防振條支承且彎曲半徑最大的傳熱管進行分析。將振型函數(shù)和熱工水力條件(橫向速度分布、流體密度分布)代入式(18)～式(20),其中管子的振型函數(shù)通過有限元方法得到,傳熱管有限元模型示意圖如圖10所示。傳熱管的典型振型圖如圖11所示。

圖1 蒸汽發(fā)生器U形彎管段的結(jié)構(gòu)和流場示意圖

圖2 作用于振動管的合速度

圖3 作用于振動管的流體力

圖4 計算域和局部網(wǎng)格細節(jié)

圖5 升力系數(shù)的空間變化

圖6 升力系數(shù)的幅值和相位差與試驗數(shù)據(jù)的比較

圖8 流管邊界提取

圖9 相位滯后預測結(jié)果

圖10 有限元模型

圖11 典型振型

利用本文建立的數(shù)值預測方法計算得到的穩(wěn)定性圖如圖12所示。并與文獻試驗結(jié)果[12]、工程驗證性試驗結(jié)果以及現(xiàn)有工程方法的結(jié)果進行對比。其中,工程驗證性試驗和工程方法均基于Connors公式,工程方法的結(jié)果采用自主研發(fā)的專用軟件計算獲得。選取了幾階典型模態(tài)對應的不穩(wěn)定比如表2所示。從表2可知,準穩(wěn)態(tài)理論計算的結(jié)果接近線性,且預測誤差較大;非定常模型與流管模型的計算結(jié)果相近,與試驗結(jié)果吻合較好;現(xiàn)有工程方法給出的結(jié)果最保守;對于所有的模態(tài),流彈不穩(wěn)定的裕量都是足夠的,工程算例的管束不會發(fā)生流彈失穩(wěn),與工程驗證性試驗的結(jié)論一致。

表2 各階模態(tài)對應的流彈不穩(wěn)定比

4 結(jié) 論

本文采用理論建模和CFD計算相結(jié)合的方式,發(fā)展適用于工程的管束結(jié)構(gòu)流彈失穩(wěn)預測的數(shù)值方法,得到以下結(jié)論:

(1) 形成了通過CFD計算流體力相關參數(shù)的辨識方法,提出了結(jié)合圖像處理技術提取流管特征的定量表征技術,通過仿真數(shù)據(jù)驅(qū)動的方式,獲得了全套流體力相關參數(shù)。

(2) 基于流彈失穩(wěn)的三種理論模型,建立了數(shù)值預測方法,完成了工程應用和對比研究,并用試驗數(shù)據(jù)進行了驗證,CFD計算與理論模型耦合的數(shù)值預測方法,具備一定的實用性。

(3) 數(shù)值預測方法結(jié)合了理論建模的嚴謹性和CFD計算可考慮實際結(jié)構(gòu)動力學特征的優(yōu)點,又避免了復雜管束流固耦合模擬對海量計算資源的需求,降低了理論模型對試驗數(shù)據(jù)的依賴,有利于在工程中推廣應用。