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

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

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

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

1 控制方程

1.1 流體控制方程

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

?u=0 ，

(1)

(2)

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

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

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

(3)

(4)

(5)

(6)

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

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

2.1 模型計(jì)算域以及邊界條件

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

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

2.2 網(wǎng)格劃分

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

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

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

為了驗(yàn)證文中數(shù)值模擬方法的可行性與最優(yōu)的網(wǎng)格劃分策略，需要對(duì)不同網(wǎng)格數(shù)下圓柱的氣動(dòng)特性進(jìn)行模擬，模擬的圓柱間距選為L/D=6，T/D=0～4，雷諾數(shù)Re=3.48×104，處于亞臨界范圍，湍流度為1%，通過y+計(jì)算得到的壁面第一層網(wǎng)格厚度為0.013 mm，選用了4種不同數(shù)量的網(wǎng)格進(jìn)行計(jì)算，并將模擬的平均升阻力系數(shù)與Wu等[15]的實(shí)驗(yàn)結(jié)果以及肖春云[16]的數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比。結(jié)果表明，隨著流域的網(wǎng)格數(shù)量從1.3×105增加到3.7×105時(shí)，圖3(a)的平均阻力系數(shù)曲線表現(xiàn)為不斷接近肖春云和Wu的曲線結(jié)果，圖3(b)的阻力系數(shù)時(shí)程曲線表現(xiàn)為上升，但是當(dāng)網(wǎng)格數(shù)量由3.7×105增加至4.5×105時(shí)，圖3(a)的2種網(wǎng)格對(duì)應(yīng)的平均阻力系數(shù)曲線非常接近，圖3(b)的阻力系數(shù)時(shí)程曲線在穩(wěn)定段幾乎重合，由此可見，繼續(xù)加密網(wǎng)格已經(jīng)對(duì)結(jié)果的精確度無明顯提升，后續(xù)動(dòng)態(tài)圓柱的模擬選用與370 000網(wǎng)格數(shù)相似的網(wǎng)格劃分策略可達(dá)到精度要求，即圓柱壁面第一層網(wǎng)格利用y+=1計(jì)算，圓柱壁面的網(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 計(jì)算參數(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%。整個(gè)數(shù)值模擬的雷諾數(shù)范圍為2 400≤Re≤28 200，始終處于亞臨界范圍，折減風(fēng)速Vr=V/fnD=5～60，湍流度為1%，時(shí)間步Δt=0.000 4 s。

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

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

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

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

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

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

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

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

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

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

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

(7)

(8)

4.1 氣動(dòng)荷載對(duì)比

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

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

圖11 氣動(dòng)力時(shí)程曲線傅里葉變換 Fig. 11 Frequency of the aerodynamic coefficients obtained by the two methods

4.2 運(yùn)動(dòng)結(jié)果對(duì)比

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

5 結(jié)束語

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

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

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

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