考慮雙向流固耦合的換熱器彈性管束流致振動特性三維數值模擬研究

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

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

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

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

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

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

1 換熱管的固有頻率

1.1 單跨管的固有頻率

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

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

(1)

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

(2)

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

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

m=mt+mi+mo

(3)

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

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

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

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

(a) 單跨管網格

1.2 四跨管的固有頻率

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

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

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

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

1.3 非等跨管的固有頻率

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

(a) “011”型

采用數值仿真模擬對三種失效型前四階固有頻率進行計算，得出相應失效狀態下振型圖，圖5為“011”振型，圖6為“101”振型，圖7為“100”振型，固有頻率計算結果如表3所示。

(a) 一階振型

(a) 一階振型

(a) 一階振型

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

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

2 換熱管的流致振動

2.1 網格獨立性及數值方法驗證

在進行數值模擬計算之前，通常要進行網格的獨立性檢驗。網格的疏密及質量直接影響著計算結果的精度。通常認為當網格密度達到一定程度后，繼續增加網格單元數量對于計算結果的影響非常小，此時可認為網格疏密對于計算結果的影響可以忽略。

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

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

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

表4 不同網格數量下的各項仿真數據對比Tab.4 Comparison of various simulation data under different grid numbers

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

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

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

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

外流速度作用下換熱器管束將產生變形，漩渦脫落對管束產生力的作用，從而管束位移呈現一定的變化規律。Khushnood等[17]通過實驗研究指出在一定殼程流速范圍內，湍流是引起管內振動的主要激勵機制。

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

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

(a) 管結構與管外流場端面網格

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

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

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

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

(a) t=0.441 s

據圖14與圖15可知，升、阻力曲線與x、y向位移曲線在時域上的變化趨勢具有相似性，達到一定穩定狀態后呈現簡諧輸出。對其穩定部分分別作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，符合阻力與升力頻率成兩倍關系理論[18]。據x、y向位移頻域曲線顯示，y和x向位移響應頻率分別與之升、阻力頻率一一對應，x向位移主要是由阻力引起，y向位移主要是由升力引起。

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

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

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

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

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

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

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

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

(a)Ⅰ跨

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

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

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

3 管束的流致振動

3.1 網格獨立性驗證

本節針對3節徑比的管束進行網格驗證。選取不同數量的四種網格進行評估。

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

表7 不同網格數量下的各項仿真數據對比Tab.7 Comparison of various simulation data under different grid numbers

3.2 節徑比1.5時管束的流致振動

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

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

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

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

圖23 xy切面網格Fig.23 xy section mesh

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

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

(a) 速度云圖

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

圖25 監測點位置圖Fig.25 Location map of monitoring points

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

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

取穩定部分，將橫坐標對應x向位移，縱坐標對應y向位移，得到管束的振動位移曲線，如圖28所示。

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

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

3.3 節徑比3時管束的流致振動

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

(a) 速度云圖

由圖29可以看出，與節徑比1.5相比，節徑比3時的管束尾流出現明顯的渦旋脫落，在5號管上取同樣的一監測點，得到該點處管子的振動位移曲線，如圖30所示。將時域較為穩定部分做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號管的振動更趨于穩定。x向與y向振動位移頻率響應呈現兩倍關系。

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

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

4 結 論

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

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

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

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

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