換熱管受湍流分布力隨機振動的虛擬激勵法

張琦， 張新玉， 張文平

(1.哈爾濱工程大學(xué) 動力與能源工程學(xué)院，黑龍江 哈爾濱 150001； 2.哈爾濱工程大學(xué) 動力裝置工程技術(shù)研究所，黑龍江 哈爾濱 150001)

湍流抖振存在于換熱器管束的整個運行過程中，以微幅振動的形式出現(xiàn)，導(dǎo)致結(jié)構(gòu)某些部位與支撐板的碰摩和疲勞損傷等累積失效。管束的微幅振動對強化換熱有積極的影響，可以激發(fā)管程工質(zhì)的脈動以提高換熱效率[1]。管束流致振動的早期研究是以漩渦脫落和流體彈性激振為主，湍流激勵為輔。針對湍流抖振的工作多是研究其對周期性激勵的影響[2-4]。計算和實驗手段的限制導(dǎo)致經(jīng)常將管束繞流場視為2D，垂直于管軸線方向的平面內(nèi)進(jìn)行研究[5-6]，包括漩渦脫落和流彈激勵。然而由于湍流激勵復(fù)雜的三維特性，在沿管軸線方向上的分布隨機，因此換熱管表面的每個空間位置都存在一個時間域上的隨機湍流力，數(shù)學(xué)上可以描述為多個隨機過程。為了研究湍流所引起的隨機振動，文獻(xiàn)[7-8]對激勵力的頻譜標(biāo)準(zhǔn)化方法做過不同的探究嘗試，并提出了受單位長度上已知功率譜的隨機載荷作用下?lián)Q熱管以一階振型和固有頻率進(jìn)行湍流抖振的撓度振幅預(yù)報經(jīng)驗公式。其中載荷的功率譜需要應(yīng)用Axisa[9]提出的等功率譜密度法來設(shè)定，該方法假設(shè)分布湍流隨機載荷是服從正態(tài)分布的各態(tài)歷經(jīng)隨機過程，進(jìn)而對其進(jìn)行頻譜的標(biāo)準(zhǔn)化，王定標(biāo)等[10]基于該方法對大型管殼式換熱器進(jìn)行了隨機響應(yīng)分析。隨著計算機技術(shù)的發(fā)展，通過流場仿真技術(shù)得到換熱管上的流體力分布，進(jìn)而利用隨機振動方法對其響應(yīng)做出預(yù)報。林家浩等[11]提出并發(fā)展的“虛擬激勵法”對于多點激勵平穩(wěn)隨機振動方程給出了高效精確的求解方法。對于復(fù)雜結(jié)構(gòu)的非平穩(wěn)隨機振動，其計算效率也得到大幅提高。其應(yīng)用領(lǐng)域目前主要是工程抗震應(yīng)用[12-14]、風(fēng)工程領(lǐng)域[15-19]、車輛工程領(lǐng)域[20]以及其他發(fā)展的領(lǐng)域[21-28]。本文以管束殼程流場的CFD仿真結(jié)果中的表面壓力作為載荷，首先驗證載荷的隨機性，并利用基于梁模型的傳統(tǒng)隨機振動計算方法和虛擬激勵法分別分析了換熱管不同位置的位移響應(yīng)功率譜。

1 換熱管隨機振動理論

1.1 換熱管的受迫振動

換熱管結(jié)構(gòu)振動分析基于將換熱管抽象為細(xì)長比較大的帶有中間支撐的兩端固支多跨Bernoulli-Euler梁模型計算，分布載荷作用下的受迫振動微分方程為：

(1)

式中：EI為梁的抗彎剛度；m為梁單位長度的質(zhì)量；c為介質(zhì)阻尼系數(shù)；f(x,t)表示時空變化的分布載荷。方程(1)在零初始條件下基于模態(tài)疊加法的時域解為：

(2)

式中：φr(x)是橫向振動的第r階模態(tài)振型；ωr為其對應(yīng)的固有圓頻率；hr(t)是第r階脈沖響應(yīng)函數(shù)。

1.2 平穩(wěn)隨機響應(yīng)的計算

Syy(f)=H*(f)H(f)SFF(f)=|H(f)|2SFF(f)

(3)

式中H(f)表示輸入點與輸出點之間的導(dǎo)納函數(shù)，可以由頻響函數(shù)求解得到。略去表達(dá)式中的自變量頻率f，輸入與輸出之間的互譜可以表示為：

SFy=HSFF

(4)

SyF=H*SFF

(5)

可以推導(dǎo)多輸入F(t)=[F1(t),F2(t),…,FM(t)]、多輸出y(t)=[y1(t),y2(t),…,yM(t)]的響應(yīng)功率譜矩陣為：

Syy=H*SFFHT

(6)

SFy=SFFHT

(7)

SyF=H*SFF

(8)

Ryy(τ)=E[[yi(t)][yj(t+τ)]T]=

hr(τ1)hs(τ2)dτ1dτ2

(9)

對上式應(yīng)用維納-辛欽關(guān)系進(jìn)行Fourier變換可得系統(tǒng)響應(yīng)的功率譜矩陣為：

(10)

式(10)稱為CQC方法，計入了所有的參振模態(tài)與耦合項。對于M自由度系統(tǒng)，使用N階模態(tài)疊加法對計算進(jìn)行截斷，CQC的計算量為N2M2次實數(shù)乘法，在工程中推薦使用一個簡化的近似方法SRSS，忽略式(10)中所有的模態(tài)耦合項：

(11)

式(11)的計算量為NM2次實數(shù)乘法。SRSS僅對于參振頻率全部為稀疏分布，且各階阻尼比很小的均質(zhì)材料結(jié)構(gòu)才可用。大多數(shù)三維結(jié)構(gòu)的參振頻率很少為稀疏分布，因此SRSS的近似誤差也引起廣泛關(guān)注和研究。

1.3 虛擬激勵法

2.抽屜中有10只相同的黑襪子和10只相同的白襪子。若你在黑暗中打開抽屜，伸手拿出襪子，請問至少要拿出幾只襪子，才能確定拿到一雙？

(12)

(13)

式(13)與式(6)的形式相同，因此與CQC具有相同的精度，且計算量僅為M2次實數(shù)乘法。此方法被稱為虛擬激勵法(pseudo excitation method，PEM)。

2 管束繞流場仿真

2.1 結(jié)構(gòu)模型、網(wǎng)格劃分與仿真設(shè)置

建立19根轉(zhuǎn)置三角排列管束繞流的三維模型，管子外徑16 mm，有效長度0.62 m，分為2跨，節(jié)徑比為1.5?；诖鬁u模擬(large eddy simulation，LES)湍流模型對管束繞流的計算[24-26]已經(jīng)成為目前管束繞流仿真的應(yīng)用趨勢，本文采用LES方法對如圖1的管束繞流模型進(jìn)行仿真。針對此計算模型共進(jìn)行了2套網(wǎng)格劃分方案：方案1，流場區(qū)域采用六面體的結(jié)構(gòu)化網(wǎng)格離散，管束表面邊界層厚度0.04 mm，并以每層4%的增長率遞增至3 mm，管束區(qū)和進(jìn)、排氣道網(wǎng)格細(xì)節(jié)如圖2所示；方案2，為了提高進(jìn)、排氣道區(qū)域網(wǎng)格質(zhì)量，在方案1的基礎(chǔ)上，管束區(qū)仍采用全結(jié)構(gòu)化劃分，進(jìn)、出口過渡區(qū)域采用非結(jié)構(gòu)化劃分。

圖1 流域三視圖和管束區(qū)橫截面示意Fig.1 3-view drawing of the flow domain and section of the tube bundle cross-flow domain

圖2 流域結(jié)構(gòu)化網(wǎng)格細(xì)節(jié)Fig.2 Structural grid details in cross-flow domains

方案1經(jīng)過計算，管束表面最大壁面y+值為3.83，為保證LES的壁面y+≤1的要求，邊界層厚度降為0.01 mm，網(wǎng)格厚度增長率與方案1相同，方案2的網(wǎng)格質(zhì)量更高，數(shù)量略有增加。表1為2種網(wǎng)格方案的對比。

表1 2種網(wǎng)格劃分方案對比Table 1 Comparison of 2 schemes of mesh

流場介質(zhì)為常溫空氣，采用有限體積法求解控制方程，并使用3-D Smagorinsky-Lilly亞格子尺度LES模型，PISO壓力速度耦合格式，PRESTO!壓力離散格式，二階迎風(fēng)動量離散格式，有界二階隱式離散化的方案應(yīng)用于此模型。在入口速度為3 m/s的亞臨界雷諾數(shù)條件下，Re為1.6×104，出口為0 Pa，其余邊界包括換熱管束外壁、殼體內(nèi)壁和擋板壁均為無滑移壁面條件，瞬態(tài)計算時間步長為10-4s，總計算時長為2 s。

2.2 流場仿真結(jié)果討論

圖3所示為2個網(wǎng)格方案的管束壁面y+值分布，方案1由于管束表面流體邊界層網(wǎng)格厚度較大，導(dǎo)致壁面y+較高，最大值為3.83；改進(jìn)后的方案2邊界層網(wǎng)格厚度為方案1的1/4，因此壁面y+可以達(dá)到0.73以下，空間分布也更加明顯。證明方案2的網(wǎng)格細(xì)化使得計算精度有了很大的提高，為后續(xù)能夠精確計算管束表面流體激勵力提供了先決條件。2個網(wǎng)格方案計算得到進(jìn)出口壓降分別為168.5 Pa和171.7 Pa。

圖3 2種網(wǎng)格方案管束表面y+結(jié)果云圖Fig.3 The contours of wall y+ on the tube bundle using the 2 schemes of mesh

考察某一跨管束上壓力的分布，圖4展示的是沿流速方向視角，迎流第1排管(1～5號管)、中間排管(11～15號管)與尾管(19號管)上的壓力云圖和壁面y+云圖，按圖1(d)中對管子的編號排列。如圖4(a)所示，迎流的第1排管束上所受的壓力分布固定，在靠近右側(cè)位置可看出壓力沒有左側(cè)高，是因為進(jìn)氣道與管束區(qū)相連的位置沒有覆蓋到整個管長，空氣繞流在下游管子上會沿軸線方向擴散分布。圖4(b)可以看出在第1排管束上y+的分布沿管軸線方向上相當(dāng)均勻，說明第1排管束的繞流流速分布基本相同，而且圖4(a)和(b)中對稱的2對管，2號與3號、4號與5號管上的壓力和流速分布對稱性很強；中心排管束如圖4(c)、(d)，對稱性完全消失，管束表面壓力分布混亂，繞流進(jìn)入充分發(fā)展階段。流場下游管束周期性載荷不再主導(dǎo)，隨機激勵成為流致振動的主要因素。

圖4 管束中第1排、中間排和尾管迎流面總壓和y+云圖Fig.4 Total pressure and wall y+ contours of the 1st array, middle array and the last tube

為了驗證計算模型的準(zhǔn)確性，本文對僅有1號管的單管繞流場也進(jìn)行了模擬，圖5(a)和(b)分別為單管繞流場速度和19管繞流場的速度，在單管尾流中有足夠的空間可以形成完整的卡曼旋渦并脫落在單管上產(chǎn)生周期性的脈動升阻力，密排管束則不然，尾流中難以產(chǎn)生固定頻率的漩渦脫落。單管仿真得到的升阻力系數(shù)時域與頻域圖如圖6所示，升力的脈動頻率為74 Hz，阻力的頻率為148 Hz，滿足文獻(xiàn)[6]中所述的阻力頻率是升力頻率的2倍關(guān)系，此時St數(shù)為0.394 67。

圖5 繞流場速度云圖Fig.5 The contours of cross flow velocity magnitude

圖6 單管仿真升阻力系數(shù)曲線Fig.6 The Cl and Cd curves of single tube simulation

在與單管仿真相同的位置和邊界條件下，密排管束中的迎流1號管的升阻力系數(shù)結(jié)果如圖7(a)～(c)所示。密排管束中已不能形成完整的旋渦，流體激勵力的能量頻帶變寬，旋渦脫落受限。圖7(d)～(f)和(g)～(i)表明處于位置對稱7、8號管，在阻力系數(shù)相近的情況下，升力系數(shù)反向，頻譜近似。圖7(j)～(l)11號中間管與圖7(m)～(o)19號尾管的頻譜中也出現(xiàn)若干主導(dǎo)頻率，然而高頻中的分量變得越來越多。

圖7 密排管束中若干典型位置圓管的升阻力系數(shù)曲線Fig.7 Curves of lift and drag coefficients of tubes at typical positions in densely packed tube bundles

考察第19號尾管，每個橫截面上升力或阻力分量合成一個合力就是該橫截面上的2個垂直的橫向作用力。尾管上的阻力和升力在時域內(nèi)的分布如圖8所示，阻力隨時間變化不大，在計算范圍內(nèi)的相關(guān)系數(shù)不低于0.99，相當(dāng)于時不變分布載荷；而每個截面的升力在時間上則是隨機的。

圖8 19號管分布阻力和升力的時空場Fig.8 The spatial-temporal field of the drag force and lift force distribution on tube No.19

3個典型位置第1跨中點、全管中點以及第2跨中點升力數(shù)據(jù)的概率分布如圖9所示，虛線為標(biāo)準(zhǔn)正態(tài)分布概率線，基本服從均勻分布的概率曲線。每個截面升力都是一個隨機過程，將升力的計算結(jié)果截取無關(guān)的10次時序樣本，每個樣本內(nèi)包含100步連續(xù)的瞬態(tài)數(shù)據(jù)，經(jīng)過計算上述3個位置隨機過程的相關(guān)函數(shù)系數(shù)如圖10所示。在時間間隔相等的情況下，相關(guān)系數(shù)處于同一水平且與計算起始時間點無關(guān)，即數(shù)據(jù)的相關(guān)函數(shù)保持平穩(wěn)。各個隨機過程的時間平均并不相等，且不等于統(tǒng)計平均，因此本文中研究的管子上M個空間載荷是廣義平穩(wěn)隨機過程，屬于多點輸入平穩(wěn)隨機問題，不同于文獻(xiàn)[9]假設(shè)的單位長度載荷為各態(tài)歷經(jīng)平穩(wěn)隨機過程。

圖9 19號管典型位置升力正態(tài)概率Fig.9 Normal probability diagrams of the lift force at typical positions of tube No.19

圖10 19號管典型位置升力相關(guān)系數(shù)等值線Fig.10 Contours of the correlation coefficient of the lift force at typical positions of tube No.19

3 換熱管隨機振動響應(yīng)計算結(jié)果分析

3.1 換熱管的固有特性

由1.2和1.3節(jié)可知，在計算多跨換熱管結(jié)構(gòu)的隨機振動響應(yīng)時，離散分布激勵的功率譜矩陣SFF和結(jié)構(gòu)自身跨點頻響矩陣H是必要條件。計算符合國標(biāo)[24]的兩等跨304不銹鋼換熱管，其有效管長0.62 m，材料密度為7 930 kg/m3，材料彈性模量為200 GPa，管外徑為16 mm，管壁厚8 mm，管內(nèi)流體密度為998 kg/m3，阻尼比為0.02。利用有限元方法計算其固有頻率和模態(tài)，計算過程中充液圓管的抗彎剛度即為金屬圓管本身的剛度，充液圓管的線質(zhì)量為單位長度的金屬管質(zhì)量、管內(nèi)流體質(zhì)量與排開的管外流體質(zhì)量之和，管外流體為空氣，質(zhì)量相對于前2項可以忽略。前10階固有頻率結(jié)果為580、841、1 484、1 878、3 919、4 545、5 936、6 701、10 226和11 223 Hz。

3.2 隨機振動響應(yīng)算例

3.2.1 載荷完全相干算例

載荷在空間上完全相干，則rank(SFF)=1。令功率譜為SFF(x1,x2,f)=I，載荷為空間上完全相干的白噪聲，對此載荷分別應(yīng)用CQC、SRSS和PEM方法截斷于10階模態(tài)進(jìn)行隨機響應(yīng)求解，3種方法分別用時441.1、42.4和4.7 s，基本滿足計算速度與模態(tài)疊加階數(shù)的倍率關(guān)系。

由上文中的仿真結(jié)果可知，換熱管受到湍流隨機力作用時的能量主要集中于500 Hz以下的低頻段，高頻模態(tài)無法受激起振，本節(jié)中應(yīng)用空間完全相干的白噪聲作為載荷，理論上的響應(yīng)即為系統(tǒng)的傳遞函數(shù)，且高階模態(tài)的響應(yīng)幅值遠(yuǎn)遠(yuǎn)低于低階模態(tài)，因此計算過程中截斷于高達(dá)11 223 Hz的十階模態(tài)僅是為了驗證計算效率。

圖11(a)、(b)分別為一跨中點位置的位移響應(yīng)自譜和1/3與2/3位置的位移響應(yīng)互譜，計算范圍為0～5 000 Hz。曲線中各有3個峰值，分別對應(yīng)于充液管第2、4、6階固有頻率，圖11(b)中CQC和PEM的結(jié)果中還存在3 730 Hz附近的反共振峰，SRSS的計算結(jié)果中則不存在。

圖11 完全相干載荷作用下一跨幾個位置的位移功率譜密度Fig.11 The auto-PSD of middle-point deformation and the cross-PSD of deformation of 1/3 and 2/3 span length positions under fully coherent load

圖11中峰值頻率處的幅值相對誤差在表2中所列，相對誤差的定義為SRSS和PEM方法與CQC方法在固有頻率處幅值的差與CQC結(jié)果的比值。對于圓管上的單點響應(yīng)自譜，PEM可以表現(xiàn)出與CQC十分接近的精度。對于兩點響應(yīng)間的互譜，在低于1 000 Hz的頻段內(nèi)，PEM的精度依舊很好；在高于1 000 Hz的頻段內(nèi)，PEM與CQC的峰值出現(xiàn)分離，并且只出現(xiàn)在反共振峰處，但是誤差仍然遠(yuǎn)遠(yuǎn)小于SRSS。

表2 完全相干載荷3種方法的計算結(jié)果和相對誤差Table 2 Calculation results of three methods and relative error of the results of SRSS and PEM

3.2.2 載荷不相干算例

對2.2節(jié)中仿真得到的尾管升力進(jìn)行功率譜計算，其rank(SFF)=M，功率譜矩陣滿秩，即載荷在空間上不相干。應(yīng)用3種方法對仿真載荷作用于充液圓管模型進(jìn)行響應(yīng)計算，得到一跨中點位置的自譜如圖12(a)所示，500 Hz以內(nèi)為載荷控制頻段見圖12(b)，用相對范數(shù)誤差[29]來表征一個頻段內(nèi)SRSS或PEM的計算精度：

(14)

式中：V為SRSS或PEM方法計算得到的某點響應(yīng)自譜;V0為CQC方法計算的該點響應(yīng)自譜。載荷控制頻段內(nèi)，SRSS的相對范數(shù)誤差為43.5%，PEM的相對范數(shù)誤差為0.39%，可見PEM的精度遠(yuǎn)遠(yuǎn)高于SRSS。

圖12(a)的高頻段則回歸系統(tǒng)的傳遞函數(shù)，CQC與PEM的響應(yīng)自譜在2 170 Hz和4 100 Hz附近存在反共振峰，SRSS的響應(yīng)自譜只有在4 545 Hz處存在與CQC相近的峰值，且精度仍然低于PEM方法。

圖12 空間不相干載荷作用下一跨中點位移功率譜密度Fig.12 The auto-PSD of middle-point deformation under incoherent load

4 結(jié)論

1) 簡化換熱器中管束繞流的CFD模型瞬態(tài)數(shù)值模擬獲取管束表面湍流激勵的瞬態(tài)分布，通過本文計算在時間上是多個服從均勻分布的廣義平穩(wěn)隨機過程，而不是文獻(xiàn)[10]中假設(shè)理想的服從高斯分布的各態(tài)歷經(jīng)隨機過程。作為平穩(wěn)隨機過程，繞流場中密排管束表面，尤其是尾部管表面的流體力主要為能量集中于低頻段的寬頻作用力。

2) 利用空間完全相干的白噪聲作用于充液換熱管上以驗證3種方法計算橫向振動的位移響應(yīng)譜的效率，PEM的計算速度約等于SRSS的N倍、CQC的N2倍，其中N為CQC和SRSS方法中參與計算的模態(tài)截斷階數(shù)，且在系統(tǒng)傳遞函數(shù)結(jié)果中，PEM在峰值頻率處幅值精度與精確算法相當(dāng)。

3) 空間不相干分布載荷作用于換熱管時，造成受迫振動在低頻段內(nèi)載荷控制的響應(yīng)功率譜或高頻段內(nèi)傳函控制的響應(yīng)功率譜，以CQC結(jié)果為參考，應(yīng)用PEM相對于SRSS方法均高精且高效，因此在已知換熱管橫向載荷譜的情況下，虛擬激勵法在計算其隨機響應(yīng)方面具有重要的工程實際意義，可以為換熱器管束無支撐跨距的設(shè)計標(biāo)準(zhǔn)提供參考依據(jù)。

本文中尚未解決折流板換熱器內(nèi)多殼程流動時管束受力存在的空間相干性問題及其造成的隨機振動響應(yīng)預(yù)報方法。