管束結(jié)構(gòu)流彈失穩(wěn)的數(shù)值預(yù)測(cè)方法研究

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

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

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

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

(1)

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

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

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

1 動(dòng)力學(xué)方程

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

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

(2)

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

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

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

(3)

橫流作用下的管束的運(yùn)動(dòng)方程可寫(xiě)為

(4)

(5)

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

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

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

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

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

(8)

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

1.3 流管模型

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

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)來(lái)考慮擾動(dòng)延遲效應(yīng),面積擾動(dòng)函數(shù)a(s,t)表示為

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

(11)

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

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

(12)

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

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

(13)

采用類(lèi)似推導(dǎo)過(guò)程,可得到控制體的動(dòng)量方程

(14)

(15)

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

(16)

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

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

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

2.1 CFD模型

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

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

(17)

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

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

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

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

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

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

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

3 臨界流速預(yù)測(cè)

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

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

(18)

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

(19)

(20)

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

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

圖2 作用于振動(dòng)管的合速度

圖3 作用于振動(dòng)管的流體力

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

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

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

圖8 流管邊界提取

圖9 相位滯后預(yù)測(cè)結(jié)果

圖10 有限元模型

圖11 典型振型

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

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

4 結(jié) 論

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

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

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

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