李岳霖，孫 超,3，艾詩欽，劉月嬋,2

(1.哈爾濱理工大學(xué) 測控技術(shù)與通信工程學(xué)院，哈爾濱 150080；2.哈爾濱工程大學(xué) 水聲工程學(xué)院,哈爾濱 150001；3.哈爾濱工程大學(xué) 動力工程及工程熱物理博士后科研流動站，哈爾濱 150001)

隨著節(jié)能技術(shù)的發(fā)展，換熱器在化工領(lǐng)域的應(yīng)用愈加廣泛，帶來了顯著的經(jīng)濟(jì)收益。然而為了滿足生產(chǎn)需要，不斷追求高熱效率，管束撓性越來越大，振動越來越敏感，易發(fā)生疲勞損壞[1-3]。因此對換熱器設(shè)備結(jié)構(gòu)進(jìn)行設(shè)計、評估和檢修時，必須對換熱管流致振動的相關(guān)特性問題進(jìn)行研究分析，以保證設(shè)備的可靠性。

殼程流體流過管束時，管束后部會形成交替脫落的漩渦，對圓柱產(chǎn)生周期性的交變應(yīng)力，打破管束原有安全狀態(tài)。長時期受到流體沖刷的影響，能量不斷積累，超過安全臨界值，易產(chǎn)生劇烈爆炸。軸向流繞流結(jié)構(gòu)體存在流致振動這種現(xiàn)象于1950年末就被發(fā)現(xiàn)，可惜并未得到足夠的重視[4]。

賴永星[5]研究指出只有在流速遠(yuǎn)高于正常流速場合，縱向流致振動才需要考慮。正常情況下，橫向流速可能引起較大振動，對結(jié)構(gòu)產(chǎn)生危害最大，人們更為關(guān)注橫向流致振動產(chǎn)生的危害；晉文娟[6]通過水洞實驗指出單柔性管進(jìn)行研究存在不足之處；吳倩倩等[7]對內(nèi)外充液管束流致振動進(jìn)行雙向流固耦合數(shù)值模擬研究，發(fā)現(xiàn)單向流固耦合數(shù)值模擬結(jié)果并不十分可靠；Zhao等[8]對相鄰兩根管束的相互影響進(jìn)行分析，得出兩相鄰管的節(jié)徑比與位移之間的關(guān)系；Katinas等[9]通過實驗分析，劉景偉等[10]通過數(shù)值計算與實驗對比分析，研究了線性排列與交錯排列管束兩種不同排列方式的振動特性；馮志鵬等[11]針對正方形順排列彈性管束流固耦合系統(tǒng)三維數(shù)值模型，實現(xiàn)了不同彈性管束模型的流體力及振動響應(yīng)特性分析。

目前，通過大量的理論與實驗研究，學(xué)術(shù)界比較認(rèn)同的管束流致振動機(jī)理有以下四種：漩渦脫落，流體彈性不穩(wěn)定，紊流抖振與聲共振[12]。現(xiàn)有研究基本是以這些理論為支撐進(jìn)行開展，大多數(shù)都是針對外流域不同形狀管束、單管及排列的不同方式對流致振動特性進(jìn)行單向耦合分析，且針對不同節(jié)徑比開展管束集群研究較少。因此，筆者就換熱器管束流致振動問題展開研究，針對包含不同節(jié)徑比管束的振動特性進(jìn)行雙向流固耦合分析，為換熱器安全運行提供重要的設(shè)計依據(jù)。

1 換熱管的固有頻率

1.1 單跨管的固有頻率

首先進(jìn)行理論值計算，計算模型如圖1，考慮到換熱管支撐條件的影響，對其結(jié)構(gòu)進(jìn)行合理簡化，即換熱管視為直管，兩端管板視為固定支撐，管內(nèi)部流體介質(zhì)為水，求解換熱管振動問題可以等效為連續(xù)梁振動問題。換熱管束固有頻率計算式[13]為

圖1 單跨管模型Fig.1 Single-span-tube model

(1)

式中：fi為i階固有頻率，Hz；l為跨距，m；EI為管的抗彎剛度；λi為由下面式(2)決定

(2)

式中：N為跨數(shù)；n=0,1,2,……,(頻帶為基數(shù)時)；n=1,2,3,……,(頻帶為偶數(shù)時)。

m為管的單位長度質(zhì)量，kg/m。計算方法如(3)所示

m=mt+mi+mo

(3)

式中：mt為單位長度空管質(zhì)量，kg/m；mi為換熱管內(nèi)單位長度流體(本文為水)的質(zhì)量，kg/m；mo為單位長度管所排開管外流體(本文為空氣)的質(zhì)量，kg/m，相較于前兩項質(zhì)量可忽略不計。換熱管材料為0Cr19Ni9鋼，密度：7 900 kg/m3，泊松比：0.247，彈性模量：195 GPa。換熱管長：l=1 200 mm，內(nèi)徑：d=15 mm，外徑：D=16 mm。將數(shù)值代入公式，計算結(jié)果列于表1中。

表1 單跨管固有頻率理論值與數(shù)值計算值Tab.1 Theoretical and numerical values of natural frequency of single span pipe

張杰等[14]對U型充液管道進(jìn)行了流固耦合模態(tài)分析，結(jié)果表明隨著模態(tài)階數(shù)的增加，管道的固有頻率逐漸增大，無耦合計算的管道固有頻率大于流固耦合狀態(tài)時的結(jié)果。本文運用數(shù)值模擬仿真方法對單跨管固有頻率進(jìn)行分析，計算出前八階固有頻率。為保證計算精度提高計算準(zhǔn)確性，流體域與單跨管結(jié)構(gòu)劃分為六面體結(jié)構(gòu)網(wǎng)格，結(jié)構(gòu)化網(wǎng)格可以很容易實現(xiàn)區(qū)域邊界擬合，最大單元：24 mm，最小單元：0.24 mm，最大單元增長率：1.3，網(wǎng)格劃分總數(shù)為44 980個。網(wǎng)格劃分結(jié)果如圖2所示。

通過表1可知，單跨管固有頻率的理論值和數(shù)值計算值隨著階次的增加誤差逐漸增大，前八階頻率誤差基本低于3%。結(jié)果表明，利用數(shù)值仿真對單跨管固有頻率計算不會引起較大誤差。

(a) 單跨管網(wǎng)格

1.2 四跨管的固有頻率

實際生產(chǎn)設(shè)計中考慮到管子熱膨脹和方便安裝等問題，換熱管與支撐板間存在間隙，管束因流體繞流后會有支撐沒有接觸上，被稱為非有效支撐。首先假設(shè)全部支撐都是有效支撐，換熱管材料不變，四跨管折流板簡化為簡支支撐，兩端管板簡化為固定支撐如圖3所示。四跨管總長：L=1 200 mm，每跨長：l=300 mm，內(nèi)徑：d=15 mm，外徑：D=16 mm。固有頻率理論計算結(jié)果和數(shù)值仿真計算結(jié)果如表2所示。

圖3 四跨管模型Fig.3 Four-span tube model

表2 四跨管固有頻率理論值與數(shù)值計算值Tab.2 Theoretical and numerical values of natural frequency of four span pipe

通過表2可以看出兩種計算方法所得固有頻率誤差仍基本低于3%，可以采用該方法對多跨管固有頻率進(jìn)行計算。相同階數(shù)下，多跨管固有頻率普遍高于單跨管固有頻率，說明管束跨數(shù)會對固有頻率計算有影響。

1.3 非等跨管的固有頻率

換熱器設(shè)備正常運行時，設(shè)備各個構(gòu)件達(dá)到最佳狀態(tài)，具有高度可靠性和穩(wěn)定性，隨著運行時間提升，設(shè)備出現(xiàn)不同程度損壞，可靠性和穩(wěn)定性得不到保證，其換熱管束中多跨管某一支撐換熱管的折流板失效時，就會出現(xiàn)非等跨的有效支撐現(xiàn)象，出現(xiàn)較大安全生產(chǎn)隱患；以四跨管為例，將失效支撐折流板用數(shù)字0代替，有效支撐折流板用數(shù)字1代替，對非等跨直管固有頻率計算分析，有三種不同位置處失效，即“011”型；“101”型；“100”型如圖4所示。

(a) “011”型

采用數(shù)值仿真模擬對三種失效型前四階固有頻率進(jìn)行計算，得出相應(yīng)失效狀態(tài)下振型圖，圖5為“011”振型，圖6為“101”振型，圖7為“100”振型，固有頻率計算結(jié)果如表3所示。

(a) 一階振型

(a) 一階振型

(a) 一階振型

表3 多跨管不同失效型固有頻率計算結(jié)果Tab.3 Calculation results of natural frequencies of different failure modes of multi span pipes

據(jù)表3計算結(jié)果顯示，相同失效支撐數(shù)量下失效位置不同的多跨管固有頻率相接近，而隨著支撐失效數(shù)的增加，同階固有頻率會降低[15]。

2 換熱管的流致振動

2.1 網(wǎng)格獨立性及數(shù)值方法驗證

在進(jìn)行數(shù)值模擬計算之前，通常要進(jìn)行網(wǎng)格的獨立性檢驗。網(wǎng)格的疏密及質(zhì)量直接影響著計算結(jié)果的精度。通常認(rèn)為當(dāng)網(wǎng)格密度達(dá)到一定程度后，繼續(xù)增加網(wǎng)格單元數(shù)量對于計算結(jié)果的影響非常小，此時可認(rèn)為網(wǎng)格疏密對于計算結(jié)果的影響可以忽略。

流場與單跨管結(jié)構(gòu)都采用六面體結(jié)構(gòu)網(wǎng)格進(jìn)行劃分。由于在近壁面速度變化梯度較大，所以本文對近壁面處網(wǎng)格進(jìn)行細(xì)化，為達(dá)到計算精度要求，本文設(shè)置了5層邊界層，第一層網(wǎng)格高度應(yīng)滿足y+≈1，因此在Re=64 000時，第一層網(wǎng)格高度δ≤0.04 mm，近壁面網(wǎng)格劃分如圖8所示。

圖8 圓柱xy面局部加密網(wǎng)格Fig.8 Local encryption grid of cylindrical xy plane

由表4可知，179 496網(wǎng)格單元與232 956網(wǎng)格單元的各項數(shù)據(jù)十分接近。因此，為了保證計算及后處理的效率和求解精度之間的平衡，采用179 496網(wǎng)格單元進(jìn)行數(shù)值模擬。

表4 不同網(wǎng)格數(shù)量下的各項仿真數(shù)據(jù)對比Tab.4 Comparison of various simulation data under different grid numbers

為了驗證數(shù)據(jù)的準(zhǔn)確性，本文與文獻(xiàn)[16]中的試驗數(shù)據(jù)進(jìn)行對比，結(jié)果如圖9所示。計算了不同外流速度沖擊下的St數(shù)。由表5可知，仿真數(shù)據(jù)與試驗數(shù)據(jù)誤差較小。因此，可以采用該數(shù)值模擬方法進(jìn)行換熱管的流致振動分析。

表5 不同外流速度下St數(shù)的對比分析Tab.5 Comparative analysis of St number under different outflow velocities

圖9 圓柱繞流的St數(shù)與雷諾數(shù)Re的關(guān)系曲線Fig.9 Relationship curve between St number and Reynolds number Re of circular cylinder flow

2.2 水作用下單跨管的雙向耦合

外流速度作用下?lián)Q熱器管束將產(chǎn)生變形，漩渦脫落對管束產(chǎn)生力的作用，從而管束位移呈現(xiàn)一定的變化規(guī)律。Khushnood等[17]通過實驗研究指出在一定殼程流速范圍內(nèi)，湍流是引起管內(nèi)振動的主要激勵機(jī)制。

為研究三維管束的流致振動規(guī)律，本文設(shè)置外流速度為4 m/s，研究單跨管在該流速沖擊下的振動特性。

管束兩端保持固定支撐，換熱管材料不變。換熱管長：L=500 mm，內(nèi)徑：d=15 mm，外徑：D=16 mm。為保證流場流動得到充分的發(fā)展，流場進(jìn)口與圓柱中心距離為2D，流場寬度為10D，流場出口距離圓柱中心8D。網(wǎng)格劃分如圖10所示。

(a) 管結(jié)構(gòu)與管外流場端面網(wǎng)格

根據(jù)雷諾數(shù)計算公式得出4 m/s速度下雷諾數(shù)為Re=64 000，流域特征處于亞臨界區(qū)，選用k-ε湍流模型可以達(dá)到計算要求，考慮到水作用下的雙向流固耦合，選取時間步長為5×10-4s，進(jìn)行瞬態(tài)的數(shù)值模擬。

圖11 速度xy截面云圖Fig.11 xy section velocity nephogram

圖12 壓力xy截面云圖Fig.12 xy section pressure nephogram

由圖13可以看出，換熱管相對于zx平面進(jìn)行往復(fù)擺動，管束中部位移量最大，運動軌跡與后文圖18相吻合，呈現(xiàn)“8”字型。

(a) t=0.441 s

據(jù)圖14與圖15可知，升、阻力曲線與x、y向位移曲線在時域上的變化趨勢具有相似性，達(dá)到一定穩(wěn)定狀態(tài)后呈現(xiàn)簡諧輸出。對其穩(wěn)定部分分別作FFT變換，即可得兩者各自頻域曲線，如圖16、圖17所示。

圖14 升、阻力時域曲線圖Fig.14 Time domain curves of lift and drag

圖15 x、y向位移時域曲線圖Fig.15 Time domain curve of displacement in x and y direction

圖16 升、阻力頻率曲線圖Fig.16 Lift and drag frequency diagram

從圖16可知，升力頻率為49.75 Hz，阻力頻率為104 Hz，符合阻力與升力頻率成兩倍關(guān)系理論[18]。據(jù)x、y向位移頻域曲線顯示，y和x向位移響應(yīng)頻率分別與之升、阻力頻率一一對應(yīng)，x向位移主要是由阻力引起，y向位移主要是由升力引起。

管束振動達(dá)到穩(wěn)定狀態(tài)后，將橫坐標(biāo)對應(yīng)x向位移，縱坐標(biāo)對應(yīng)y向位移，繪制三維管束振動位移曲線，得到振動位移曲線圖，如圖18所示?？梢园l(fā)現(xiàn)位移曲線呈現(xiàn)“8”字型，這是因為x向位移響應(yīng)頻率是y向位移響應(yīng)頻率的兩倍，且振幅較大，出現(xiàn)典型的Strouhal振型，位移曲線圖呈現(xiàn)上下對稱的兩個部分。

圖18 振動位移曲線Fig.18 Vibration displacement curve

改變外流域流速大小，得到不同速度下的振動位移曲線如圖19所示，可以看出三種速度下同樣出現(xiàn)Strouhal振型，振動位移曲線幅值與外流速成正比關(guān)系，隨流速增大而增大。

圖19 不同速度下振動位移曲線Fig.19 Vibration displacement curves at different speeds

2.3 水作用下四跨管的雙向耦合

跨度較長的換熱器設(shè)備在理想條件下折流板視為簡支有效支撐，換熱器管束兩端保持著固定有效支撐。

將2.2中單跨管簡化等分為四跨管，中間折流板視為簡支支撐，兩端管板視為固支支撐模型，從上到下每跨分別命名為Ⅰ跨、Ⅱ跨、Ⅲ跨、Ⅳ跨如圖20所示。四跨管同樣取來流速度為4 m/s的水對其進(jìn)行沖擊，不同跨振動位移曲線如圖21所示。

圖20 四跨管模型圖Fig.20 Four span pipe model diagram

(a)Ⅰ跨

據(jù)圖21所示，互為對稱跨的振動位移曲線“8”字型相似且位移幅度相等，位處中間兩跨的位移幅度要大于兩側(cè)跨。

表6與圖16、圖17所得數(shù)據(jù)對比可知，相同流速沖擊下，多跨管振動頻率與單跨管基本一致。四跨管升、阻力頻率滿足兩倍關(guān)系，位移響應(yīng)頻率與之升、阻力頻率存在一一對應(yīng)關(guān)系。

表6 各跨升、阻力頻率及軸向位移響應(yīng)頻率計算結(jié)果Tab.6 Calculation results of rise,resistance frequency and axial displacement response frequency of each span

3 管束的流致振動

3.1 網(wǎng)格獨立性驗證

本節(jié)針對3節(jié)徑比的管束進(jìn)行網(wǎng)格驗證。選取不同數(shù)量的四種網(wǎng)格進(jìn)行評估。

由表7可知，303 256網(wǎng)格單元與344 215網(wǎng)格單元的各項數(shù)據(jù)十分接近。因此，為了保證計算及后處理的效率和求解精度之間的平衡，數(shù)值模擬計算時采用303 256網(wǎng)格進(jìn)行。

表7 不同網(wǎng)格數(shù)量下的各項仿真數(shù)據(jù)對比Tab.7 Comparison of various simulation data under different grid numbers

3.2 節(jié)徑比1.5時管束的流致振動

三維單管能從一定程度上反映振動特性，但并未考慮其他管的影響，存在一定的局限性。受管與管之間的距離影響，每根管的流場不能得到充分發(fā)展，管間流體會相互影響，相互作用，導(dǎo)致?lián)Q熱器管束流場和受力與三維單管存在差異。

選取五根管組成的管束群，探討不同節(jié)徑比下管束的流致振動特性，分別對節(jié)徑比為1.5和3.0的兩種管束群進(jìn)行建模和網(wǎng)格劃分。每根換熱管模型和材料與2.2中相同，流域為水，流場寬度為8D，流場長度為20D，每兩根相鄰換熱管間距為1.5D。計算模型如圖22所示。

圖22 管束計算模型Fig.22 Tube bundle calculation model

整個流場中管束周圍存在較大的速度梯度變化，為使計算結(jié)果更加精確，對管束周圍進(jìn)行加密細(xì)化網(wǎng)格處理，流域與換熱管結(jié)構(gòu)選擇六面體結(jié)構(gòu)化網(wǎng)格劃分，換熱管網(wǎng)格最大單元：3.29×10-3m，最小單元：5.06×10-5m，最大增長率：1.05；流域網(wǎng)格最大單元：9.36×10-3m，最小單元：1.01×10-3m，最大增長率：1.1，網(wǎng)格總數(shù)為303 256。網(wǎng)格劃分如圖23所示。

圖23 xy切面網(wǎng)格Fig.23 xy section mesh

取來流速度為4 m/s，采用湍流k-ε模型，瞬態(tài)求解器，時間步長取5×10-4s。在T=1 s時刻，xy切面結(jié)果如圖24所示。

由圖24可知，當(dāng)外流橫向掠過管束時，管束之間由于距離太近，尾流之間相互影響，流場整體流動要比單管時更加混亂[19]。

(a) 速度云圖

在5號管中間段上取一監(jiān)測點，如圖25所示。監(jiān)測該點處管子的振動位移，并得到如圖26所示的管子振動位移時域曲線，將時域曲線較穩(wěn)定部分做FFT變換可知5號管振動位移響應(yīng)頻率，如圖27所示。

圖25 監(jiān)測點位置圖Fig.25 Location map of monitoring points

圖27 x、y方向振動位移頻域曲線圖Fig.27 Frequency domain curves of vibration displacement

從時域及頻域曲線可以看出，管束距離太近，5號管尾流發(fā)展受到其他管影響，導(dǎo)致x、y方向的振動頻率并不存在兩倍關(guān)系，并出現(xiàn)多個峰值，73.17 Hz時x向的振幅要大于y向。

取穩(wěn)定部分，將橫坐標(biāo)對應(yīng)x向位移，縱坐標(biāo)對應(yīng)y向位移，得到管束的振動位移曲線，如圖28所示。

由圖28可以看出，由于管束間距離較近，各管束尾流間相互影響，管束的振動位移曲線并不規(guī)律，Strouhal振型并不明顯。

圖28 管束的振動位移曲線Fig.28 Vibration displacement curve of tube bundle

3.3 節(jié)徑比3時管束的流致振動

改變管束間的距離為3D，流域不變，計算條件與3.2相同，得到t=1 s時刻，xy切面速度和壓力云圖如圖29所示。

(a) 速度云圖

由圖29可以看出，與節(jié)徑比1.5相比，節(jié)徑比3時的管束尾流出現(xiàn)明顯的渦旋脫落，在5號管上取同樣的一監(jiān)測點，得到該點處管子的振動位移曲線，如圖30所示。將時域較為穩(wěn)定部分做FFT變換得出管子的頻域曲線，如圖31所示。

圖30 x、y方向振動位移時域曲線圖Fig.30 Time domain curves of vibration displacement

圖31 x、y方向振動位移頻域曲線圖Fig.31 Frequency domain curves of vibration displacemen

從時頻域曲線可知，由于管束間距變大，管束間影響變小，5號管的振動更趨于穩(wěn)定。x向與y向振動位移頻率響應(yīng)呈現(xiàn)兩倍關(guān)系。

取穩(wěn)定部分，將橫坐標(biāo)對應(yīng)x向位移，縱坐標(biāo)對應(yīng)y向位移得到振動位移曲線，如圖32所示?？梢园l(fā)現(xiàn)，5號管尾部形成完整的漩渦脫落[20]，尾流幾乎未受到影響，呈現(xiàn)明顯的Strouhal振型。其他管子由于受到相互影響并未呈現(xiàn)明顯的“8”字型，但仍比節(jié)徑比1.5時更接近“8”字型。

圖32 管束的振動位移曲線Fig.32 Vibration displacement curve of tube bundle

4 結(jié) 論

本文以換熱器管束為研究對象，針對換熱器管束固有頻率、單管及管束群流致振動分析，結(jié)合換熱管固有頻率理論計算值與仿真計算值進(jìn)行對比，誤差在3%左右。對單跨管、多跨管流致振動特性及五管束兩種不同節(jié)徑比(1.5和3.0)流致振動數(shù)值模擬分析，計算結(jié)果表明：

(1)多跨管的同階固有頻率隨著支撐失效個數(shù)增多而降低。

(2)相同數(shù)量但不同位置的失效支撐基本不會改變管的固有頻率大小。

(3)單跨管與多跨管振動位移曲線呈現(xiàn)“8”字型，出現(xiàn)典型的Strouhal振型，且幅值隨著流體速度增大而增大。多跨管的振動頻率與單跨管基本一致，且對稱跨“8”字型曲線相似，位移幅值相等。

(4)管束間的距離會對流場產(chǎn)生較大影響，當(dāng)節(jié)徑比為1.5時，管束距離較小，管子的振動頻率較為雜亂。當(dāng)節(jié)徑比為3.0時，管束間距離變大，各管間流體影響較小，振動更加穩(wěn)定。受影響最小的5號管振動位移頻率響應(yīng)呈現(xiàn)兩倍關(guān)系，振動位移曲線呈現(xiàn)明顯的Strouhal振型。