傅亨仁，王靈芝，晏致濤，孫 毅

(1.重慶大學(xué) a.土木工程學(xué)院;b.山地城市建設(shè)新技術(shù)教育部重點實驗室，重慶 400044；2.重慶科技學(xué)院 建筑工程學(xué)院，重慶 401331)

除了定義為尾流致振，也有學(xué)者將上述振動定義為尾流馳振。Brika 等[9]對固定的上游圓柱和自由運動的下游圓柱進行實驗(雷諾數(shù)Re=5 000～27 000，兩圓柱流向間距比L/D=7～25，折減風(fēng)速Vr=4～11)，在L/D=7.5和8時，發(fā)現(xiàn)下游圓柱經(jīng)歷渦激共振與尾流馳振的組合振動現(xiàn)象。King 等[10]通過對2個自由振動的串列圓柱進行實驗，在流向間距L/D=2.5，Vr>11時，觀察到尾流馳振現(xiàn)象。Hover 等[11]對處于固定圓柱后的彈性支撐的下游圓柱(僅允許橫向振動)進行風(fēng)洞實驗，在串聯(lián)距離L/D=4.75時，觀察到高振幅的尾流馳振現(xiàn)象，且該振幅隨折減風(fēng)速不斷增加，在實驗的折減風(fēng)速范圍不出現(xiàn)峰值。Tokoroa 等[12]利用風(fēng)洞對串列布置(L/D=4.3～8.7，Vr=0～25)的雙圓柱進行全尺度實驗，在L/D=4.3時，觀察到下游圓柱的尾流馳振現(xiàn)象，對應(yīng)振幅隨著折減風(fēng)速無限制增大，該現(xiàn)象于L/D=6.5時消失。

上述研究中，雖然對雙圓柱尾流馳振振動響應(yīng)展開了實驗，但是并未獲得尾流馳振發(fā)生時的氣動荷載，數(shù)值模擬可以方便地得到尾流馳振全過程的氣動荷載并加以分析。根據(jù)Zdravkovich等[13]的研究，兩圓柱流向間距比L/D=2和橫流向間距比T/D=1時，處于尾流干擾的區(qū)域，文中采用基于RANS的SSTk-ω湍流模型，假定上游圓柱固定，對彈性支撐兩自由度的下游圓柱振動特性和氣動荷載進行研究。利用ICEM對流域進行結(jié)構(gòu)化網(wǎng)格劃分，結(jié)合動網(wǎng)格、滑移網(wǎng)格技術(shù)以及用戶自定義接口編程，將計算結(jié)構(gòu)響應(yīng)的Newmark-β代碼嵌入Fluent軟件進行數(shù)值模擬，折減風(fēng)速范圍為Vr=5～60，觀察下游圓柱的尾流馳振特性，并將對應(yīng)的氣動荷載模擬結(jié)果與準(zhǔn)定常數(shù)值計算結(jié)果進行對比分析。

1 控制方程

1.1 流體控制方程

二維不可壓縮均勻粘性牛頓流體運動的基本控制方程為連續(xù)性方程和Navier-Stokes方程，密度為ρ、動力粘度為μ的流體域的控制方程為

?u=0 ，

(1)

(2)

式中：p為壓力；u為流速矢量，包括流向x方向與橫向y方向的流速分量；ρ表示空氣的密度；μ表示空氣的動力粘度。

1.2 結(jié)構(gòu)控制方程

上游圓柱固定不動，尾流下運動圓柱不考慮扭轉(zhuǎn)自由度時，圓柱在2個方向的運動為往復(fù)運動，接近簡諧振動，在其運動過程中會受到阻尼力與彈性恢復(fù)力，可將其簡化成雙自由度的彈簧振子模型，該模型在運動時，除了要滿足上述的流體控制方程外，還需滿足如下運動方程：

(3)

(4)

(5)

(6)

式中：由于是二維圓柱，圓柱長度l為單位長度；V為來流速度；ρ為流體密度；mx,cx和kx為圓柱流向單位長度的質(zhì)量、阻尼和剛度；my,cy和ky為圓柱橫流向單位長度的質(zhì)量、阻尼和剛度；D為圓柱的直徑；x(t)和y(t)為圓柱順流向與橫流向在t時刻的位移；V為來流速度；FD(t)和FL(t)為圓柱在t時刻受到的阻力和升力；CD(t)和CL(t)為在t時刻對應(yīng)的阻力系數(shù)與升力系數(shù)。

2 計算模型及設(shè)定

2.1 模型計算域以及邊界條件

計算域阻塞率是設(shè)計流域計算范圍的主要參數(shù)，其定義為結(jié)構(gòu)總迎風(fēng)面積與計算域風(fēng)場寬度的比值，圓柱繞流主要以側(cè)面繞流為主，阻塞率過大會對圓柱側(cè)面的流場湍流度以及圓柱表面風(fēng)壓產(chǎn)生不利影響，導(dǎo)致計算結(jié)果不精確，根據(jù)方平治等[14]研究的阻塞率對風(fēng)場數(shù)值模擬的影響結(jié)果，阻塞率取值2.5%～7%的計算域在計算風(fēng)工程中有重要參考意義。文中選定的阻塞率為5%，考慮到湍流的充分發(fā)展，采用的計算域為60D×40D，如圖1所示，上游圓柱中心距離入口邊界為25D，下游圓柱中心距離出口邊界為33D。計算域的邊界條件設(shè)為：流域入口邊界設(shè)為速度入口邊界條件(velocity-inlet)，流域出口邊界設(shè)為壓力出口邊界條件(pressure-outlet)，上下邊界定義為對稱邊界條件(symmetry)，圓柱周圍采用無滑移壁面(wall)。

圖1 計算域和邊界條件Fig. 1 Computational domain and boundary conditions

2.2 網(wǎng)格劃分

文中利用ICEM CFD對計算域進行網(wǎng)格劃分，如圖2所示。網(wǎng)格為非均勻四邊形結(jié)構(gòu)化網(wǎng)格，由于導(dǎo)線附近的流場變化劇烈，對圓柱近壁面的網(wǎng)格進行加密處理，通過控制壁面網(wǎng)格高度減小網(wǎng)格對數(shù)值計算的影響，第1層網(wǎng)格高度通過y+=1計算。

圖2 計算域網(wǎng)格和圓柱壁面網(wǎng)格Fig. 2 Global mesh and grid near the rigid wall

2.3 網(wǎng)格精確性驗證

為了驗證文中數(shù)值模擬方法的可行性與最優(yōu)的網(wǎng)格劃分策略，需要對不同網(wǎng)格數(shù)下圓柱的氣動特性進行模擬，模擬的圓柱間距選為L/D=6，T/D=0～4，雷諾數(shù)Re=3.48×104，處于亞臨界范圍，湍流度為1%，通過y+計算得到的壁面第一層網(wǎng)格厚度為0.013 mm，選用了4種不同數(shù)量的網(wǎng)格進行計算，并將模擬的平均升阻力系數(shù)與Wu等[15]的實驗結(jié)果以及肖春云[16]的數(shù)值計算結(jié)果進行對比。結(jié)果表明，隨著流域的網(wǎng)格數(shù)量從1.3×105增加到3.7×105時，圖3(a)的平均阻力系數(shù)曲線表現(xiàn)為不斷接近肖春云和Wu的曲線結(jié)果，圖3(b)的阻力系數(shù)時程曲線表現(xiàn)為上升，但是當(dāng)網(wǎng)格數(shù)量由3.7×105增加至4.5×105時，圖3(a)的2種網(wǎng)格對應(yīng)的平均阻力系數(shù)曲線非常接近，圖3(b)的阻力系數(shù)時程曲線在穩(wěn)定段幾乎重合，由此可見，繼續(xù)加密網(wǎng)格已經(jīng)對結(jié)果的精確度無明顯提升，后續(xù)動態(tài)圓柱的模擬選用與370 000網(wǎng)格數(shù)相似的網(wǎng)格劃分策略可達到精度要求，即圓柱壁面第一層網(wǎng)格利用y+=1計算，圓柱壁面的網(wǎng)格增長率為1.05，除了壁面以外的網(wǎng)格增長率為1.08，結(jié)構(gòu)化長方形網(wǎng)格的長邊尺寸與短邊尺寸最大比例控制在5以內(nèi)，相鄰block之間的網(wǎng)格尺寸應(yīng)盡量保持一致。

圖3 網(wǎng)格依賴性研究Fig. 3 Grid sensitivity research

3 尾流下彈性支撐圓柱馳振特性分析

3.1 計算參數(shù)

文中模擬的流場與圓柱參數(shù)為：空氣密度ρ=1.225 kg/m3，導(dǎo)線直徑D=30 mm，質(zhì)量比m*=m/(1/4πρD2L)=2.18，導(dǎo)線自振頻率為fn=9.33 Hz，阻尼比ξ=0.964%。整個數(shù)值模擬的雷諾數(shù)范圍為2 400≤Re≤28 200，始終處于亞臨界范圍，折減風(fēng)速Vr=V/fnD=5～60，湍流度為1%，時間步Δt=0.000 4 s。

3.2 尾流下下游圓柱的馳振特性分析

圖4給出了下游圓柱的橫向無量綱振幅隨折減風(fēng)速的變化曲線，在尾流馳振發(fā)生之前，最大無量綱振幅發(fā)生在Vr=7.5時，振幅值可達到0.13D，再結(jié)合圖5所示，Vr=7～9時，渦脫頻率與自然頻率比值fs/fn=0.97～0.99，可認為下游圓柱發(fā)生了渦激共振。本節(jié)模擬的Re與Socker等[17]的實驗保持一致，保證了流體的運動相似，為了更好地觀察到錯列布置下游圓柱的渦激振動特性，文中選用的圓柱質(zhì)量與Socker不同，所以在Vr=7～9時，Socker的渦激共振振幅較小，而文中渦激共振現(xiàn)象較為明顯。由圖4可以看出，在折減風(fēng)速Vr為38時，橫向的無量綱振幅突然以接近直線趨勢上升，說明在該折減風(fēng)速下游圓柱受上游圓柱的尾流影響開始發(fā)生尾流馳振，直到折減風(fēng)速Vr達到50時，直線上升斜率開始減小，本節(jié)的模擬結(jié)果與Socker的實驗結(jié)果吻合較好，驗證了所采用的數(shù)值模擬方法在研究尾流馳振響應(yīng)上的可行性。

圖4 下游圓柱橫向振幅隨Vr變化曲線Fig. 4 Variation of the dimensionless amplitude Ay/Dwith the reduced velocity Vr

此外，Qin等[18]和Wu等[19]的單圓柱渦激振動和文中模擬的下游圓柱振動的無量綱頻率隨折減風(fēng)速的變化曲如圖5所示，由于St=fsD/V,Vr=V/(fnD)，fs為渦脫主頻，D為圓柱直徑，V為來流風(fēng)速，Vr為折減風(fēng)速，fn為自然頻率。將兩式相乘可得St=fs/(fnVr)，曲線斜率的變化表示為斯托羅哈數(shù)的變化，由圖5可知，單圓柱的St幾乎為一個定常數(shù)0.2，而下游圓柱的St小于0.2，并且在渦激共振之后與尾流馳振起振之前的區(qū)域有較大變化，而文獻[18-19]中在渦激共振之后無明顯變化，St的值與圓柱的渦脫頻率息息相關(guān)，表明上游圓柱的尾流在下游圓柱產(chǎn)生的隨機渦脫弱于來流直接在單圓柱產(chǎn)生的渦脫，即在渦激共振區(qū)尾流抑制了下游圓柱表面的隨機渦脫。

圖5 無量綱頻率隨Vr變化曲線Fig. 5 Variation of the frequency ratio fs/fnwith the reduced velocity Vr

為了更清楚地獲得尾流馳振的振動特性，將Vr為42、45、46、50和60對應(yīng)下游圓柱的運動軌跡繪制在圖6中，可以明顯發(fā)現(xiàn)，橢圓長軸隨折減風(fēng)速的增加會朝流向傾斜再逐漸趨于穩(wěn)定，值得注意的是，所有的運動軌跡都為逆時針，說明尾流馳振具有明確方向性和自限性的風(fēng)致響應(yīng)的振動；定義θ為橢圓長軸與流向的夾角，圖7表明在發(fā)生尾流馳振后，θ隨折減風(fēng)速的增加一直在減小，當(dāng)Vr到達60時，θ為27°，與兩圓柱圓心連線與流向的夾角26.5°幾乎相等，這些現(xiàn)象表明，下游圓柱尾流馳振的最終橢圓運動軌跡長軸會與兩圓柱圓心連線重合。

圖6 尾流馳振風(fēng)速范圍下游圓柱運動軌跡Fig. 6 Trajectory of the downstream cylinder under wake galloping region

圖7 θ隨折減風(fēng)速的變化曲線Fig. 7 Variation of θ with Vr under wake galloping region

圖8和圖9為渦激共振區(qū)Vr=7.5與尾流馳振區(qū)Vr=50所對應(yīng)的運動軌跡和升力頻譜結(jié)果。圖8(a)表明，受上游圓柱的渦脫和自身的渦脫影響，下游圓柱在“鎖定區(qū)”的振動軌跡呈現(xiàn)為橢圓環(huán)的形式，已有研究對二維單圓柱兩自由度的渦激共振進行模擬，發(fā)現(xiàn)單圓柱的渦激共振響應(yīng)軌跡為8字形的Lissajou圖，說明尾流會改變下游圓柱的運動軌跡。由圖9(a)可知，發(fā)生尾流馳振后，下游圓柱的運動軌跡仍然保持著橢圓環(huán)的形式，說明下游圓柱從流體中吸取的能量與自身耗散的能量最終會達到一個動態(tài)平衡的過程，尾流馳振同渦激共振一樣是一種具有自限性的風(fēng)致響應(yīng)。圖8(b)和9(b)表明，在渦激共振區(qū)，漩渦脫落主頻除了有接近于自振頻率的一階頻率9.3 Hz，還有一個二階頻率18.6 Hz，是非線性振動產(chǎn)生的倍頻，頻率成分較為單一；而在尾流馳振區(qū)，除了漩渦脫落主頻還有大量的其他頻率成分，包括自振頻率的倍頻和其他渦脫頻率，說明尾流馳振是一種渦脫模式復(fù)雜并且?guī)в袕姺蔷€性的振動響應(yīng)。所有發(fā)生了尾流馳振的其他風(fēng)速也都有著類似的頻譜成分，由圖5可知，這些渦脫頻率值的大小隨折減風(fēng)速的增加呈線性增加，說明尾流馳振區(qū)域圓柱的St為一個定值。值得注意的是，尾流馳振高頻成分占比高于低頻成分，這些頻率會對運動響應(yīng)造成影響。

圖8 Vr=7.5下游圓柱運動軌跡和升力系數(shù)頻譜Fig. 8 Trajectory and frequency spectrum of CL at vortex-induced vibration reduced velocity Vr=7.5

圖9 Vr=50下游圓柱運動軌跡和升力系數(shù)頻譜Fig. 9 Trajectory and frequency spectrum of CL at wake galloping oscillation reduced velocity Vr=50

4 流固耦合模擬與準(zhǔn)定常數(shù)值方法對比

準(zhǔn)定常數(shù)值計算方法是求解尾流馳振振動響應(yīng)的傳統(tǒng)方法，準(zhǔn)定常理論是將流經(jīng)微振動細結(jié)構(gòu)的氣流假設(shè)為定常流，依據(jù)相對運動原理，視物體為靜止?fàn)顟B(tài)，在建立尾流馳振力學(xué)模型時忽略了流體與結(jié)構(gòu)運動的相互反饋作用。利用準(zhǔn)定常方法求解下游圓柱的動力響應(yīng)時，先擬合得到下游圓柱靜態(tài)的氣動力系數(shù)CD,CL與位置的分布函數(shù)，如式(7)所示，是x和y的多項式函數(shù)，根據(jù)文獻[15]以及考慮到擬合曲面的整體光順性要求，文中選用擬合的最高次數(shù)為6，代入式(8)求解獲得圓柱的運動位移x(t)和y(t)，再將x(t)和y(t)重新代回式(7)，最終得到圓柱運動過程中氣動力系數(shù)CD,CL隨時間的變化曲線。以上過程求解式(8)時氣動力系數(shù)CD,CL是圓柱處于靜態(tài)時(其周圍的流場為定常流)得到的，故稱之為準(zhǔn)定常方法。

(7)

(8)

4.1 氣動荷載對比

圖10給出了使用2種方法得到的下游圓柱的升力系數(shù)時程和阻力系數(shù)時程。由圖可見，2種方法的氣動力周期都幾乎相同，在一個特定周期t=10.20～10.31 s里，2種方法的阻力系數(shù)在10.20～10.22 s以及10.28～10.31 s很接近，而2種方法的升力系數(shù)在整個周期具有明顯的差異。對氣動力進行FFT傅里葉變換，如圖11所示，由頻譜可知，2種方法都存在較為明顯的自然頻率的倍頻2fn-4fn，氣動力的差異主要表現(xiàn)在準(zhǔn)定常數(shù)值計算方法對較高次倍頻和渦脫頻率的貢獻考慮不足。

圖10 兩種方法得到的氣動力時程曲線Fig. 10 Time history of the aerodynamic coefficients obtained by the two methods

圖11 氣動力時程曲線傅里葉變換 Fig. 11 Frequency of the aerodynamic coefficients obtained by the two methods

4.2 運動結(jié)果對比

圖12為2種求解方法順流向與橫流向的位移時程，可以知道，準(zhǔn)定常數(shù)值計算結(jié)果與文中模擬結(jié)果在波峰處達到最大誤差，x方向位移相對誤差為2%，y方向位移相對誤差為6%，其余部分曲線幾乎完全重合。因此，即使氣動力在較高次倍頻和渦脫頻率有非常大的差別時，運動響應(yīng)也保持著較高的吻合性，說明高倍頻的自激力和渦激力對運動響應(yīng)的貢獻非常小，這些力在下游圓柱的整個運動過程中做的總功非常小，自振頻率的前四階倍頻的自激力對尾流馳振的位移和速度起主要控制作用，在一定程度上說明尾流馳振是一種自激振動，高倍頻的自激力占比較小，所以影響較小，而渦激力占比較大，但是由于力與位移的相位關(guān)系，只在尾流馳振起振階段提供初始動力，在整個馳振過程對圓柱做的總功非常小，存在氣動力差異較大但是運動響應(yīng)結(jié)果相差不大的現(xiàn)象。如果需要準(zhǔn)確得到下游圓柱的氣動力，文中的流固耦合模擬方法可以提供較好的途徑。

5 結(jié)束語

文中利用ICEM對流域進行結(jié)構(gòu)化網(wǎng)格劃分，結(jié)合動網(wǎng)格及滑移網(wǎng)格技術(shù)以及用戶自定義接口編程，將計算結(jié)構(gòu)響應(yīng)的Newmark-β代碼嵌入Fluent軟件，利用流固耦合方法對雙圓柱間距L/D=2、T/D=1，在折減風(fēng)速Vr=5～60的范圍內(nèi)模擬研究了圓柱尾流振動特性，并將準(zhǔn)定常數(shù)值方法的結(jié)果與模擬結(jié)果進行對比，可以得到主要結(jié)論：

1)文中的流固耦合數(shù)值模擬結(jié)果與Socker的實驗結(jié)果有著較高的吻合度，驗證了文中采用的數(shù)值模擬方法在研究尾流馳振響應(yīng)上的可行性。尾流馳振階段對應(yīng)的下游圓柱自振頻率與渦脫主頻相差較大，與已有結(jié)論一致。相比于單圓柱，尾流擴大了斯托羅哈數(shù)的不穩(wěn)定風(fēng)速區(qū)域，在渦激共振區(qū)尾流抑制了下游圓柱表面的隨機渦脫。

2)尾流馳振發(fā)生后，下游圓柱的橫向振幅隨折減風(fēng)速的增大以接近直線上升，并且都呈現(xiàn)為橢圓形的逆時針運動軌跡，橢圓長軸與流向的夾角會先減小再穩(wěn)定，穩(wěn)定角θ與兩圓柱圓心連線與流向的夾角幾乎相等，說明尾流馳振是一種有自限性和明確方向性的風(fēng)致失穩(wěn)響應(yīng)。

3)2種方法得到氣動力系數(shù)的周期幾乎相同，都存在較為明顯的自然頻率的倍頻2fn-4fn，氣動力的差異主要表現(xiàn)在準(zhǔn)定常數(shù)值計算方法對較高次倍頻和渦脫頻率的貢獻考慮不足。高倍頻的自激力和渦激力對運動響應(yīng)做的貢獻非常小，自振頻率的前四階倍頻的自激力對尾流馳振的位移和速度起主要控制作用，在一定程度上說明尾流馳振是一種自激振動。