雙彈性管流固耦合振動(dòng)的數(shù)值模擬

馮志鵬，臧峰剛，張毅雄

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

管束類結(jié)構(gòu)的流固耦合振動(dòng)在海洋工程、航天航空、核工程等領(lǐng)域中非常普遍。運(yùn)行經(jīng)驗(yàn)和科學(xué)研究結(jié)果表明，流致振動(dòng)是導(dǎo)致蒸汽發(fā)生器傳熱管破裂失效的關(guān)鍵因素，由其導(dǎo)致的傳熱管破裂將會(huì)影響核電站運(yùn)行的可靠性。因此，對(duì)流體誘發(fā)的彈性管振動(dòng)進(jìn)行深入研究十分必要。

不同排列的雙管是實(shí)際應(yīng)用中最基本的組合單元，Sumner等[1]、Price等[2]通過實(shí)驗(yàn)研究了串列與交錯(cuò)排列雙管的渦激振動(dòng)；Antoine等[3]采用梁方程模擬直管，聯(lián)合Star-CD研究了單圓柱與流體的相互作用；Liu等[4]采用有限元方法和可變形網(wǎng)格技術(shù)對(duì)Re=200的并列彈性支撐剛性雙圓柱渦激振動(dòng)進(jìn)行了數(shù)值模擬；Lam等[5]采用表面渦方法對(duì)Re=2.67×104下小節(jié)徑比一列和錯(cuò)列的陣列圓柱的流致振動(dòng)進(jìn)行模擬，指出在大幅值振動(dòng)時(shí)流致振動(dòng)對(duì)橫向振動(dòng)的影響不可忽略；Longattea等[6]提出了一種采用ALE方法的流致振動(dòng)數(shù)值模型，但該模型只能預(yù)測(cè)橫流中彈性支撐剛性管的振動(dòng)頻率，無法考慮湍流與管的變形。已有的這些研究，主要針對(duì)彈性支撐剛性管，而對(duì)三維彈性管束與流體相互作用的研究相對(duì)較少。當(dāng)考慮管束的彈性變形時(shí)，由于流體流動(dòng)和管道運(yùn)動(dòng)間的相互作用、流場(chǎng)與結(jié)構(gòu)的相互作用變得更復(fù)雜，振動(dòng)特性與流場(chǎng)特征通常與節(jié)徑比、排列方式、流動(dòng)速度等相關(guān)。

本文采用雙向流固耦合方法，結(jié)合動(dòng)網(wǎng)格控制技術(shù)，同時(shí)考慮流體-結(jié)構(gòu)兩個(gè)物理場(chǎng)間的相互作用，聯(lián)合計(jì)算流體力學(xué)和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)建立三維彈性管束的流固耦合振動(dòng)模型，研究雙彈性管的流固耦合振動(dòng)，旨在為自主設(shè)計(jì)蒸汽發(fā)生器傳熱管束的流致振動(dòng)分析提供技術(shù)支持。

1 數(shù)值模型

1.1 CFD模型

在預(yù)測(cè)管束的橫掠問題中，較先進(jìn)的RANS(雷諾平均N-S方程)湍流模型，如可實(shí)現(xiàn)k-ε和RNGk-ε模型均極大地低估了密集管束中的湍動(dòng)能強(qiáng)度，而LES(大渦模擬)的結(jié)果與DNS(直接模擬)和實(shí)驗(yàn)結(jié)果吻合較好[7]，可得到RANS無法獲得的湍流結(jié)構(gòu)，因此本文采用大渦模擬方法來計(jì)算流體區(qū)域。連續(xù)性方程及動(dòng)量守恒方程在物理空間經(jīng)濾波后可得到大渦模擬的控制方程：

(1)

(2)

流場(chǎng)區(qū)域及網(wǎng)格如圖1所示，其中D為管子外徑，流體域網(wǎng)格為六面體結(jié)構(gòu)化網(wǎng)格。圖1中左邊入口采用速度入口邊界條件，右端出口采用壓力出口邊界條件，其他外邊界按對(duì)稱邊界和固定壁面處理；管壁面為流固耦合交界面，設(shè)為動(dòng)網(wǎng)格邊界。

1.2 結(jié)構(gòu)分析模型

管的結(jié)構(gòu)參數(shù)取為：管長(zhǎng)0.5 m；外徑D=0.01 m；內(nèi)徑0.009 5 m；彈性模量E=1×105MPa；泊松比υ=0.3；密度ρs=6 500 kg/m3；阻尼系數(shù)α=5.098，β=0.000 215。

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

利用有限元方法對(duì)管結(jié)構(gòu)進(jìn)行離散，即可得到每個(gè)管的質(zhì)量矩陣與剛度矩陣[8]。采用Newmark數(shù)值積分方法求解結(jié)構(gòu)動(dòng)力學(xué)方程，進(jìn)行瞬態(tài)動(dòng)力學(xué)分析：

(3)

1.3 動(dòng)網(wǎng)格模型

采用基于擴(kuò)散光順(Diffusion-Based Smooth)方法的動(dòng)網(wǎng)格模型來控制運(yùn)動(dòng)邊界的網(wǎng)格更新，其網(wǎng)格運(yùn)動(dòng)由式(4)的擴(kuò)散方程控制。

·(γ

(4)

式中：u為網(wǎng)格移動(dòng)的速度，在變形邊界上網(wǎng)格的運(yùn)動(dòng)與邊界相切；γ為擴(kuò)散系數(shù)，表示為單元體積的函數(shù)，γ=1/Vα，用于控制邊界網(wǎng)格對(duì)內(nèi)部網(wǎng)格的影響，α為控制參數(shù)，本文取2。式(4)用有限體積法離散并迭代求解，單元中心的移動(dòng)速度u通過距離加權(quán)平均插值到各網(wǎng)格節(jié)點(diǎn)，根據(jù)式(5)即可更新節(jié)點(diǎn)位置：

xnew=xold+uΔt

(5)

1.4 流體-結(jié)構(gòu)間的耦合

通過流固耦合交界面進(jìn)行固體域和流體域間的數(shù)據(jù)傳遞，將計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(CSD)和計(jì)算流體動(dòng)力學(xué)(CFD)聯(lián)合。CSD通過計(jì)算結(jié)構(gòu)位移來指定流體域的變形，而用CFD計(jì)算作用在結(jié)構(gòu)上的載荷，結(jié)構(gòu)域和流體域間的雙向耦合通過采用相同時(shí)間步的流體計(jì)算和結(jié)構(gòu)動(dòng)力學(xué)計(jì)算依次迭代實(shí)現(xiàn)。

1) 用Newmark算法求解結(jié)構(gòu)動(dòng)力學(xué)方程，得到結(jié)構(gòu)的位移、速度等參數(shù)，將結(jié)構(gòu)的位移與速度傳遞給動(dòng)網(wǎng)格求解器。

2) 動(dòng)網(wǎng)格求解器根據(jù)結(jié)構(gòu)位移更新流體域網(wǎng)格。

3) CFD求解，得到流場(chǎng)的u、v、w和壓力場(chǎng)p等流場(chǎng)參數(shù)，并將作用在結(jié)構(gòu)上的流體載荷傳遞給結(jié)構(gòu)求解器。

4) 回到1，進(jìn)行下一時(shí)間步求解。

流體力雖不是一連續(xù)函數(shù)，但可在每個(gè)時(shí)間步結(jié)束時(shí)，通過流體求解器得到一瞬時(shí)值，從而更新當(dāng)前結(jié)構(gòu)的位置，并得到下一時(shí)間步的流體計(jì)算網(wǎng)格。

2 結(jié)果分析

為表述方便，定義以下無量綱變量：升力系數(shù)均方根ClRMS=FlRMS/(0.5ρAU2)，阻力系數(shù)均方根CdRMS=FdRMS/(0.5ρAU2)，其中FdRMS、FlRMS分別為阻力Fd與升力Fl的均方根值；流向位移x/D、橫向位移y/D；流向振幅Ax/D=xRMS/D，橫向振幅Ay/D=yRMS/D，其中xRMS、yRMS分別為流向振動(dòng)位移與橫向振動(dòng)位移的均方根值；渦脫頻率St=fstD/U；間隙流速Upr=Up/(fnD)，Up=UP/(P-D)，其中P/(P-D)為節(jié)徑比，U為來流速度，D為管外徑，P為兩管的間距，fn為管的固有頻率，fst為靜止圓柱體繞流的瀉渦頻率；A為所要計(jì)算方向上的投影面積。

2.1 模型驗(yàn)證

為驗(yàn)證數(shù)值模型，計(jì)算了單根彈性管的流致振動(dòng)。圖2中分別就數(shù)值計(jì)算得到的管的振動(dòng)頻率fex與固有頻率fn之比fex/fn及橫向振幅Ay/D隨來流速度(無量綱)Ur(Ur=U/fnD)的變化情況與已有實(shí)驗(yàn)結(jié)果[3]進(jìn)行了比較，通過比較可看到，兩者在趨勢(shì)和數(shù)值上吻合較好，說明本文的流體-結(jié)構(gòu)交互作用模型是可行的和正確的。

2.2 臨界間距

1) 兩串列管

為研究?jī)纱泄艿慕换プ饔?，建立了如圖3所示的兩串列管，其中管2位于管1尾部的下游，流向間距Px的變化范圍為1.2D～4D。

圖4示出了阻力系數(shù)均方根CdRMS、升力系數(shù)均方根ClRMS與振幅隨節(jié)徑比的變化情況。管1、管2與單管間的相互比較分析可發(fā)現(xiàn)：(1) 振幅隨節(jié)徑比的變化趨勢(shì)與流體力隨節(jié)徑比的變化趨勢(shì)相同，兩串列管的臨界間距為2D。(2) 當(dāng)Px/D>2時(shí)，管1幾乎不受管2的影響，其流體力、振幅與單管接近，管2的CdRMS遠(yuǎn)小于管1與單管的，但隨節(jié)徑比的增加，管2的CdRMS迅速增加到約為單管的0.5倍；由于管1尾流的影響，管2的ClRMS與Ay/D均大于單管的，且隨節(jié)徑比的增大而增大，當(dāng)Px/D=4時(shí)，管2的ClRMS約是管1的4倍，振幅約是管1的5倍。(3) 當(dāng)Px/D≤2時(shí)，兩串列管的ClRMS與振幅均接近單管，而兩串列管的CdRMS則小于單管的，管1約為單管的0.8倍，管2僅為單管的0.1倍。

圖2 振動(dòng)響應(yīng)隨Ur的變化情況

圖3 串列管示意圖

圖4 阻力系數(shù)、升力系數(shù)均方根和振幅隨節(jié)徑比的變化

顯然，流體力、振幅與節(jié)徑比密切相關(guān)，當(dāng)節(jié)徑比小于臨界值(Px/D≤2)時(shí)，管1與管2的橫向變量(升力、橫向振幅)均較小，接近于單管，而當(dāng)節(jié)徑比大于臨界值時(shí)，管2的橫向變量遠(yuǎn)大于管1與單管的；下游管2與上游管1間的Ay/D隨節(jié)徑比的增大而增大。

2) 兩并列管

為研究橫流中兩并列管間的相互作用，建立了如圖5所示的參數(shù)完全相同的兩并列管，橫向間距記為Py。由于并列管中管1與管2的升力與橫向振幅相等、相位相反，阻力與流向振幅基本相同，因此只取管1進(jìn)行分析。

圖5 并列管示意圖

圖6a示出兩并列管CdRMS、ClRMS隨節(jié)徑比的變化，圖6b顯示了振動(dòng)幅值隨節(jié)徑比的變化，并分別與單管的進(jìn)行了比較。可看到，CdRMS和ClRMS隨Py/D的減小而增加，振幅隨間距的增加而呈非線性趨勢(shì)減小，同時(shí)還發(fā)現(xiàn)當(dāng)間距大于一臨界間距(3D)后，兩管間基本不再相互影響，ClRMS與Ay/D與單管的接近。當(dāng)間距較小時(shí)，兩管的流場(chǎng)強(qiáng)烈耦合，漩渦激發(fā)振動(dòng)和流體彈性不穩(wěn)定性均可能導(dǎo)致兩管相互碰撞。

2.3 臨界流速

通過對(duì)不同節(jié)徑比的兩串列管、兩并列管的流致振動(dòng)分析發(fā)現(xiàn)，并列管、串列管的振動(dòng)特性會(huì)在其臨界間距發(fā)生較大的轉(zhuǎn)變，以下分別以P/D=1.6與P/D=3為例分析流速與兩彈性管的流固耦合振動(dòng)間的關(guān)系，其中P/D=1.6的兩并列管會(huì)在Upr>4時(shí)發(fā)生碰撞。

圖7為串列管與并列管的CdRMS和ClRMS隨Upr的變化。從圖7a可看出，CdRMS隨Upr的增大先減小再增加到一較恒定值。從圖7b可看出，P/D=3的兩串列管、P/D=1.6與P/D=3的兩并列管ClRMS隨Upr的增大先減小再增加，達(dá)到最大值后迅速下落；而對(duì)P/D=1.6的兩串列管，其ClRMS隨Upr的增大而增大，達(dá)到最大值后也迅速下落。

圖8為串列管與并列管的流向振幅Ax/D與橫向振幅Ay/D隨Upr的變化。通過綜合比較圖7、8可得出，P/D=3的串列管1、并列管的CdRMS、ClRMS、Ax/D、Ay/D與串列管2的ClRMS、Ay/D開始急劇增加的速度為Upr=3，只

有處于下游串列管2的流向變量——CdRMS與Ax/D開始大幅增加的速度發(fā)生在Upr=4。P/D=1.6的串列管與并列管的CdRMS最小值均在Upr=4，ClRMS、Ax/D、Ay/D開始急劇增加的速度也在Upr=4。因此，P/D=3的雙管的臨界流速為Upr=3，P/D=1.6的雙管的臨界流速為Upr=4。通過仔細(xì)比較串列管1、串列管2與并列管在相同條件下的流體力系數(shù)及振幅，可發(fā)現(xiàn)，雙管的臨界速度與節(jié)徑比相關(guān)，而排列方式對(duì)其影響不大。相同節(jié)徑比情況下，并列管的振動(dòng)較串列管強(qiáng)烈，受到的流體力更大，振幅也大于串列管的。

綜上所述，當(dāng)節(jié)徑比大于臨界節(jié)徑比(如P/D=3)時(shí)，兩彈性管的流體力系數(shù)隨Upr的變化曲線存在一最小值，其恰好對(duì)應(yīng)振幅開始急劇增加的流速，此時(shí)的流速即為臨界流速。在臨界流速以下，流體力系數(shù)隨間距的減小而增大，而在臨界流速以上，CdRMS則隨Upr的增加先減小再增大到一較恒定值，ClRMS隨Upr的增加先減小再增大，在達(dá)到最大值后迅速下落。

當(dāng)節(jié)徑比小于臨界節(jié)徑比時(shí)(如P/D=1.6)，CdRMS隨Upr的變化曲線存在一最小值，其隨Upr的增加先減小再增大到一較恒定值，而ClRMS隨Upr的增加而增大，達(dá)到最大值后迅速下落。

圖6 阻力系數(shù)、升力系數(shù)均方根(a)和振幅(b)與節(jié)徑比的關(guān)系

圖7 CdRMS(a)和ClRMS(b)隨Upr的變化

圖8 Ax/D和Ay/D隨Upr的變化

流向與橫向振幅均隨Upr的增加而增加。

2.4 運(yùn)動(dòng)軌跡

通過流體力系數(shù)均方根及振幅響應(yīng)隨Upr的變化可知，P/D=3和P/D=1.6兩種節(jié)徑比管束的臨界流速分別為Upr=3和Upr=4，因此以下分別取低于臨界速度和高于臨界速度兩種情況進(jìn)行分析。對(duì)P/D=3，取Upr=2.2與Upr=3.8；而對(duì)P/D=1.6，由于兩并列管在Upr>4時(shí)發(fā)生碰撞，因此取Upr=2.7與Upr=4。

1) 兩串列管

對(duì)不同間距的兩串列管來說，當(dāng)流速低于臨界流速時(shí)，串列管的振動(dòng)很小，運(yùn)動(dòng)軌跡紊亂；當(dāng)流速高于臨界流速時(shí)，兩管開始在流向和橫向發(fā)生大幅振動(dòng)，運(yùn)動(dòng)軌跡為“8”字形圖，如圖9所示，阻力使運(yùn)動(dòng)軌跡向下游方向彎曲，上游管1的振幅較下游管2的大。另外，從圖9還可發(fā)現(xiàn)，當(dāng)間距大于臨界間距后，兩管的運(yùn)動(dòng)軌跡與相同情況下單管的類似，兩管間的相互影響較弱，而當(dāng)間距小于臨界間距時(shí)，兩管間的相互作用明顯較強(qiáng)烈，如圖9a、b所示。

2) 兩并列管

不同間距兩并列管在臨界流速內(nèi)的響應(yīng)與串列管類似，但在臨界流速外的運(yùn)動(dòng)軌跡與串列管明顯不同，且與節(jié)徑比相關(guān)。臨界間距內(nèi)(如P/D=1.6)的兩并列管在較高流速時(shí)的運(yùn)動(dòng)軌跡為如圖9c所示的橢圓形；而臨界間距外的兩并列管(如P/D=3)，其在高于臨界流速的流體作用下的運(yùn)動(dòng)軌跡如圖9d所示，與相同情況下兩串列管的不同，其軌跡不是“8”字形，而為類似“小雨滴”的形狀。從圖9c、d中也可看到，兩并列管的橫向運(yùn)動(dòng)相位相反，幅值相等。

2.5 尾渦結(jié)構(gòu)

1) 兩串列管

P/D=1.6和P/D=3兩串列管的尾渦結(jié)構(gòu)示于圖10。從圖10a、b可看出，在臨界間距內(nèi)的兩串列管的流場(chǎng)場(chǎng)結(jié)構(gòu)與單管的類似，漩渦間距較小，強(qiáng)度也較低，尾渦尺度隨流速的增加而有所增大，上游管1無穩(wěn)定的渦脫落。從圖10c、d可看出，當(dāng)間距大于臨界間距時(shí)，兩管均有渦脫落，且尾渦結(jié)構(gòu)會(huì)在臨界流速發(fā)生顯著改變，在低速時(shí)，串列管2后的尾渦模態(tài)與單管的類似，但當(dāng)流速較高時(shí)，在串列管2的后面出現(xiàn)了兩排平行的渦，尾渦的橫向間距明顯增大。

a、c——P/D=1.6,Upr=4；b、d——P/D=3,Upr=3.8

2) 兩并列管

P/D=1.6和P/D=3兩并列管的尾渦結(jié)構(gòu)示于圖11。從圖11a、b可看出，當(dāng)Upr=2.7時(shí)，漩渦分別在管1與管2各自尾部的1.5D范圍內(nèi)形成和脫落，但沿著流向的發(fā)展，其尾渦相互干擾并合并，在1.5D后，尾渦合并為單排的、交替的漩渦。從圖11c、d可看出，流速較低時(shí)，在管1與管2的后面形成幾乎獨(dú)立的渦街，而當(dāng)流速較高時(shí)，管1與管2會(huì)在其尾部分別形成兩排平行的旋渦，且兩管后面的渦會(huì)在尾流中間相互擠壓和干擾，渦的橫向間距也明顯增大。從尾渦結(jié)構(gòu)中也可看到并列管尾渦相位相反的現(xiàn)象。

圖10 P/D=1.6(a、b)和P/D=3(c、d)兩串列管的尾渦結(jié)構(gòu)

圖11 P/D=1.6(a、b)和P/D=3(c、d)兩并列管的尾渦結(jié)構(gòu)

3 結(jié)論

1) 兩串列管的臨界間距為2D，兩并列管的臨界間距為3D。2) 在臨界間距內(nèi)，兩彈性管相互影響，而當(dāng)間距大于臨界間距后，兩管間的相互影響基本消失，響應(yīng)與單管的接近；兩串列管的下游管與上游管間的橫向振幅比值隨間距的增大而增大，兩并列管的流體力與振幅隨著間距的增加呈非線性趨勢(shì)減小。3) 當(dāng)流速高于臨界流速時(shí)，兩串列管的運(yùn)動(dòng)軌跡均為“8”字形圖，而間距大于臨界間距的兩并列管的運(yùn)動(dòng)軌跡為類似“小雨滴”形，間距小于臨界間距的兩并列的運(yùn)動(dòng)軌跡為橢圓形；當(dāng)流速低于臨界流速時(shí)，兩串列管與兩并列管的振動(dòng)均很小，運(yùn)動(dòng)軌跡紊亂。4) 在臨界流速或臨界間距內(nèi)，串列的下游管、并列的兩管的尾渦均為獨(dú)立的渦街；而在臨界流速與臨界間距以外，則會(huì)在串列的下游管與并列的兩管尾部分別形成兩排平行的渦街。

參考文獻(xiàn)：

[1] SUMNER D, RICHARDS M D, AKOSILE O O. Two staggered circular cylinders of equal diameter in cross-flow[J]. Journal of Fluids and Structures, 2005, 20: 255-276.

[2] PRICE S J, PAIDOUSSIS M P, KRISHNAMOORTHY S. Cross-flow past a pair of nearly in-line cylinders with the upstream cylinder subjected to a transverse harmonic oscillation[J]. Journal of Fluids and Structures, 2007, 23: 39-57.

[3] ANTOINE P, SIGIST J F, HAMDOUNI A. 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] LIU Y, SO R M C, LAU Y L, et al. Numerical studies of two side-by-side elastic cylinders in a cross-flow[J]. Journal of Fluids and Structures, 2001, 15: 1 009-1 030.

[5] LAM K, JIANG G D, LIU Y, et al. Simulation of cross-flow-induced vibration of cylinder arrays by surface vorticity method[J]. Journal of Fluids and Structures, 2006, 22: 1 113-1 311.

[6] LONGATTEA E, BENDJEDDOUB Z, SOULIB M. Methods for numerical study of tube bundle vibrations in cross-flows[J]. Journal of Fluids and Structures, 2003, 18: 513-528.

[7] BENHAMADOUCHE S, LAURENCE D. LES, coarse LES, and transient RANS comparisons on the flow across a tube bundle[J]. International Journal of Heat and Fluid Flow, 2003, 24: 470-479.

[8] 王勖成. 有限單元法[M]. 北京：清華大學(xué)出版社,2003:468-495.