低雷諾數(shù)中等間距串列雙方柱渦激振動的數(shù)值模擬

杜曉慶，邱 濤，鄭德乾，趙 燕

(1.上海大學(xué) 土木工程系，上海 200444; 2.上海大學(xué) 風(fēng)工程和氣動控制研究中心，上海 200444; 3.河南工業(yè)大學(xué) 土木建筑學(xué)院，鄭州 450001)

渦激振動現(xiàn)象時常發(fā)生，并易對結(jié)構(gòu)造成疲勞破壞，因而受到廣泛關(guān)注[1].與單方柱相比，雙方柱渦激振動的影響因素眾多[2]，其干擾機理尚不明確.

對于固定串列雙方柱繞流，Sohankar[3]將串列雙方柱的柱間流態(tài)分為單一鈍體流態(tài)、剪切層再附流態(tài)(L/B=2～4，其中L為柱心間距，B為方柱邊長)和雙渦脫流態(tài)三類.學(xué)者將雙方柱從剪切層再附流態(tài)向雙渦脫流態(tài)轉(zhuǎn)變的柱心間距比稱為臨界間距比，Sohankar[3]得出的臨界間距比為L/B=4.Yen等[4]發(fā)現(xiàn)雙方柱在剪切層再附流態(tài)時的繞流場特性非常復(fù)雜，并得出雙方柱的臨界間距比為L/B=3～5.

對于單方柱渦激振動，Zhao等[5]發(fā)現(xiàn)當(dāng)流向角為0°時，單方柱在低雷諾數(shù)下并無明顯的“鎖定”現(xiàn)象.Sen等[6]得出低質(zhì)量比的單方柱會出現(xiàn)“弱鎖定”現(xiàn)象，即方柱在振動鎖定區(qū)內(nèi)的振動頻率與自振頻率的比值fy/fn遠小于1.0；“弱鎖定”現(xiàn)象發(fā)生時，單方柱常有大幅橫流向振動；隨著質(zhì)量比增加，方柱的“弱鎖定”現(xiàn)象會消失，而出現(xiàn)“鎖定”現(xiàn)象，即fy/fn≈1.0.

針對雙方柱振動問題，往年研究并不充分.Kumar等[7-8]發(fā)現(xiàn)在上游方柱運動、下游方柱固定的情況下，當(dāng)兩方柱的柱心間距較小(L/B=1.25)且處于串列布置時，上游方柱出現(xiàn)橫流向最大振幅.Mithun等[9]則在雷諾數(shù)Re=100時，通過固定方柱振幅、改變橫流向振動頻率，研究了串列方柱的氣動力和流場結(jié)構(gòu)特征，發(fā)現(xiàn)當(dāng)L/B=2和3時，下游方柱的升力隨時間推移變化紊亂.Pillalamarri等[10-11]則在Re=200、m*=m/(ρB2)=5(m為方柱質(zhì)量，ρ為流體密度)和L/B=4情況下，對上游方柱固定、下游方柱運動的串列雙方柱進行了數(shù)值模擬研究.Guan等[12]則在Re=200、m*=10時，研究了剛性連接并列雙方柱的渦激振動特性.

學(xué)者對雙方柱在中小間距和臨界間距時的振動問題均較為關(guān)注，未見到在此間距下，研究上、下游方柱均可自由振動的雙方柱渦激振動的相關(guān)文獻.本文在Re=150時，對L=2B和4B的串列雙方柱渦激振動進行了數(shù)值研究.

1 數(shù)值方法

1.1 計算模型和控制方程

將雙方柱渦激振動簡化為質(zhì)量-彈簧-阻尼系統(tǒng)(見圖1)，圖中UC和DC分別表示上、下游方柱，V為均勻來流速度.入口邊界與上下邊界距原點O均為30B；出口邊界距原點O為40B.

對于流體域，假定方柱周圍流場為二維不可壓縮流體，其連續(xù)方程和N-S方程分別為：

(1)

(2)

(3)

式中：p為流體壓力，μ為黏性系數(shù)，u、v分別為x、y方向的速度，F(xiàn)x、Fy分別為流體在x、y方向的體力.

對于固體域，兩個方柱均考慮順流向和橫流向運動，其控制方程為：

(4)

(5)

式中：X和Y分別為方柱的順流向與橫流向位移，c和k分別為結(jié)構(gòu)阻尼與剛度，F(xiàn)D(t)和FL(t)分別為作用在方柱上的順流向與橫流向流動力.

數(shù)值計算中，對于流體域，采用層流模型，通過 SIMPLEC法求解速度與壓力耦合，通過二階精度離散格式求解式(2)、(3)；對于固體域，式(4)、(5)的求解采用四階龍格庫塔方法.網(wǎng)格更新采用彈簧光滑法.

數(shù)值模擬參數(shù)Re=150，L=2B和4B，m*=10；為節(jié)省計算資源，c取0；折減速度Vr=V/(fnB)=1～12，其中結(jié)構(gòu)自振頻率fn在x和y兩個方向保持一致，且在計算中保持V不變，而通過改變fn實現(xiàn)不同Vr的計算.由于L=2B的雙方柱在Vr>13時出現(xiàn)相互碰撞，因此本文僅分析Vr=1～12的工況.

1.2 邊界條件及網(wǎng)格劃分

圖1還給出了網(wǎng)格模型的邊界條件，本文將入口邊界、兩側(cè)壁面、出口邊界和方柱表面分別設(shè)為速度入口、對稱邊界、自由出流和無滑移壁面.

圖1 計算模型和計算域

雙方柱網(wǎng)格方案見圖2，對計算域網(wǎng)格進行分塊處理，在方柱近壁面內(nèi)采用結(jié)構(gòu)化網(wǎng)格；圓形框內(nèi)的網(wǎng)格采用非結(jié)構(gòu)化，框外則為結(jié)構(gòu)化.通過加密距離方柱表面較近處的網(wǎng)格來捕捉其周圍復(fù)雜的流場，本文也考慮了近壁面最小網(wǎng)格厚度的影響(見表1).為減小網(wǎng)格更新計算量，僅將圓形框區(qū)域設(shè)定為動網(wǎng)格區(qū)域.

圖2 網(wǎng)格方案

2 網(wǎng)格檢驗和結(jié)果驗證

為驗證計算模型的合理性，分別對單方柱繞流和渦激振動進行網(wǎng)格獨立性檢驗和結(jié)果驗證.

2.1 單方柱繞流的結(jié)果驗證

以雷諾數(shù)為150的靜止單方柱繞流為算例，進行了網(wǎng)格和步長獨立性驗證，并與文獻結(jié)果進行對比(表1)，表中St、CL,rms、CD,mean分別為斯特羅哈數(shù)、脈動升力系數(shù)和平均阻力系數(shù)，A1～A4和B1～B3工況分別為周向和徑向網(wǎng)格獨立性檢驗，Z1～Z3工況為無量綱時間步(Δt*=ΔtV/B，Δt為計算時間步長)獨立性驗證.可見，本文所有工況的計算結(jié)果較為接近，也與文獻結(jié)果吻合較好.考慮到模擬精度和計算效率，后續(xù)網(wǎng)格模型均采用A3工況的參數(shù).

表1 結(jié)果驗證與網(wǎng)格參數(shù)

2.2 單方柱渦激振動的結(jié)果驗證

采用A3工況的網(wǎng)格，將Re=100、m*=3的單方柱渦激振動的計算結(jié)果與文獻進行比較(見圖3)，圖中Ay/B(Ay=(Ymax-Ymin)/2)、fx/fn和fy/fn分別為橫流向無量綱振幅、順流向振動頻率比和橫流向振動頻率比，其中fx、fy分別為方柱的順、橫流向振動頻率；Ymin、Ymax分別為橫流向最小位移和最大位移.由圖3可知本文結(jié)果非常吻合文獻結(jié)果.

圖3 方柱渦激振動結(jié)果驗證

3 計算結(jié)果及分析

3.1 振動響應(yīng)特性

圖4、5為雙方柱無量綱振幅Ay/B、Ax/B(Ax=(Xmax-Xmin)/2)和振動頻率比fy/fn隨Vr的變化曲線，為便于比較，圖中還給出了單方柱的計算結(jié)果，其中Xmin和Xmax為順流向最小位移和最大位移；上游柱(UC)參數(shù)腳標(biāo)設(shè)為1，下游柱(DC)設(shè)為2.可見：

1)不同柱心間距雙方柱和單方柱的振動均以橫流向為主.與雙方柱相比，單方柱Ay較小且無明顯的“弱鎖定”和“鎖定”現(xiàn)象，這與文獻[5]結(jié)論相同.

2)當(dāng)L=2B時，上下游方柱的振動鎖定區(qū)為Vr=9.75～12，并伴有“弱鎖定”現(xiàn)象發(fā)生，其對應(yīng)的fy/fn≈0.87；兩方柱在振動鎖定區(qū)內(nèi)出現(xiàn)Ay最大值，其值分別為Ay1=1.04B和Ay2=0.70B.當(dāng)L增大至4B時，上游方柱無大幅橫流向振動，因此未發(fā)生“鎖定”和“弱鎖定”現(xiàn)象，而下游方柱的振動鎖定區(qū)為Vr=5.5～6.5，且發(fā)生“鎖定”現(xiàn)象，對應(yīng)的fy/fn≈0.94；兩方柱Ay最大值均出現(xiàn)在振動鎖定區(qū)外，其值分別為Ay1=0.21B和Ay2=0.67B.

圖4 雙方柱振幅隨折減速度的變化

圖5 雙方柱橫流向振動的頻率比

值得注意的是，不同柱心間距的下游方柱在振動鎖定區(qū)外也存在大幅橫流向振動，其振動區(qū)域分別為Vr=6.5～7(L=2B)和Vr=7～9(L=4B)，下文將對此進一步分析.此外，本文將下游方柱的橫流向振幅曲線與文獻[11]進行比較.發(fā)現(xiàn)隨著Vr的增加，文獻[11]的下游方柱橫流向振幅也會出現(xiàn)兩個峰值，但振幅與本文差異較大，這可能是上游方柱運動狀態(tài)的不同造成的.

為進一步研究雙方柱的振動特性，圖6給出了雙方柱在振動鎖定區(qū)內(nèi)的位移時程.可見，對于順流向，L=2B和4B的上下游方柱的振動相位隨時間推移而不斷變化.對于橫流向，L=2B的兩方柱振動相位為同相位，其相位角ψ≈38°；L=4B的兩方柱振動相位則為反相位，其ψ≈152°.

圖6 振動鎖定區(qū)內(nèi)上下游方柱位移時程

3.2 氣動力特性

3.2.1 平均阻力系數(shù)

當(dāng)柱心間距隨Vr變化發(fā)生改變時，將對雙方柱的振動響應(yīng)和流場結(jié)構(gòu)產(chǎn)生影響，因此，圖7、8給出了雙方柱CD,mean和平均柱心間距Lmean隨Vr的變化規(guī)律.

圖7 平均阻力系數(shù)隨折減速度的變化

當(dāng)L=2B時，在振動鎖定區(qū)(Vr=9.75～12)內(nèi)，CD,mean1值為正且數(shù)值急劇增加，致使上游方柱沿順流向產(chǎn)生較大偏移；此時，上游方柱對下游方柱的“遮擋”作用[16]造成CD,mean2<0，致使下游方柱逆流向上產(chǎn)生偏移，最終導(dǎo)致Lmean的急劇減小(見圖8).

圖8 平均柱心間距隨折減速度的變化

當(dāng)L=4B時，CD,mean1隨Vr增加基本不變，而CD,mean2隨Vr的增加逐漸接近CD,mean1，致使上下游方柱順流向偏移量相近，從而導(dǎo)致Lmean基本不變.

3.2.2 脈動升力系數(shù)

圖9為CL,rms隨Vr的變化曲線，可見，當(dāng)L=2B時，在振動鎖定區(qū)內(nèi)(Vr=9.75～12)，兩方柱CL,rms迅速增加，且CL,rms1>CL,rms2，因而兩方柱出現(xiàn)Ay最大值，且上游方柱Ay最大值大于下游方柱；下游方柱在振動鎖定區(qū)外(Vr=6.5～7)也會有較大的CL,rms2，此時下游方柱的fy與fn較為接近(見圖5)，其升力更易激勵下游方柱振動，因而出現(xiàn)大幅橫流向振動.當(dāng)L=4B時，上游方柱的CL,rms1隨Vr增加變化幅度較??；下游方柱在振動鎖定區(qū)內(nèi)、外均會出現(xiàn)較大的CL,rms2，其中在振動鎖定區(qū)外(Vr=7～9)，下游方柱的fy與fn較為接近(見圖5)，因此出現(xiàn)較大橫流向振幅.

圖9 脈動升力系數(shù)隨折減速度的變化

3.3 能量輸入分析

本節(jié)將結(jié)合瞬時能量輸入過程，分析下游方柱Ay2在各振動區(qū)達到極值時的工況.圖10為下游方柱能量輸入、升力系數(shù)與橫流向位移的時程.其中P*=CL(t)v(t)/V為無量綱功率，v(t)為方柱橫流向運動速度.

對于振動鎖定區(qū)外的情況(圖10(a)、(d))，L=2B和4B時，下游方柱的CL與Y隨時間推移均呈穩(wěn)定的相位變化，致使升力對下游方柱能量輸入穩(wěn)定.且此時下游方柱的fy與fn相近，其升力極易激勵下游方柱振動，因而下游方柱橫流向振幅在振動鎖定區(qū)外出現(xiàn)明顯極值.

對于振動鎖定區(qū)內(nèi)的情況(圖10(b)、(c))，當(dāng)L=2B時，雖然下游方柱CL具有較大的幅值，但由于上游方柱尾渦的影響，下游方柱CL的周期性較弱，造成CL與Y的相位不斷變化，從而擾亂并減弱了升力對下游方柱的能量輸入，因此，下游方柱橫流向振幅雖達到最大值，但小于上游方柱.當(dāng)L=4B時，下游方柱的CL與Y相位穩(wěn)定，且CL出現(xiàn)較大幅值，造成其對下游方柱穩(wěn)定輸入強勁能量，因而下游方柱在振動鎖定區(qū)內(nèi)出現(xiàn)大幅橫流向振動.

圖10 下游方柱的升力系數(shù)、橫流向位移和能量輸入時程

3.4 流場結(jié)構(gòu)特征

圖11為雙方柱尾流模態(tài)及柱間流態(tài)隨Vr增加的演變過程.可見，共出現(xiàn)柱間流態(tài)2種、尾流模態(tài)6種.表2歸納了雙方柱尾流模態(tài)與柱間流態(tài)的關(guān)系.

在剪切層再附流態(tài)內(nèi)，L=2B和4B時的雙方柱尾流模態(tài)均為模態(tài)Ⅰ～Ⅲ.其中模態(tài)Ⅰ下游方柱背風(fēng)面存在較長的尾流；模態(tài)Ⅱ和模態(tài)Ⅲ分別為剪切層再附流態(tài)內(nèi)的平行渦街模態(tài)和“2S”模態(tài).當(dāng)方柱處于雙渦脫流態(tài)時，L=2B時僅出現(xiàn)模態(tài)Ⅴ，此時下游方柱背風(fēng)面尾渦極不穩(wěn)定；L=4B的則為模態(tài)Ⅳ和Ⅵ，其分別為雙渦脫流態(tài)內(nèi)的平行渦街模態(tài)和“2S”模態(tài).

結(jié)合上文知，L=2B時的雙方柱振動鎖定區(qū)(Vr=9.75～12)出現(xiàn)在雙渦脫流態(tài)內(nèi)，此時的尾流為模態(tài)Ⅴ；L=4B時的下游方柱振動鎖定區(qū)(Vr=5.5～6.5)出現(xiàn)在剪切層再附流態(tài)內(nèi)，此時的尾流為模態(tài)Ⅱ.兩種柱心間距的下游方柱橫流向振幅在振動鎖定區(qū)外也有明顯極值，其尾流分別為模態(tài)Ⅱ(L=2B，Vr=6.5)和Ⅳ(L=4B，Vr=8.5).

表2 不同柱心間距的柱間流態(tài)和尾流模態(tài)

圖11 雙方柱柱間流態(tài)和尾流模態(tài)

3.5 流固耦合機制

結(jié)合圖10，為研究導(dǎo)致下游方柱振動的流固耦合機制，圖12、13給出了下游方柱出現(xiàn)大幅振動時的瞬時渦量圖及對應(yīng)的能量輸入時程.其中，圖12為剪切層再附流態(tài)下，兩種柱心間距的Ay2分別在振動鎖定區(qū)內(nèi)(L=4B)和區(qū)外(L=2B)出現(xiàn)明顯極值時的典型工況；圖13為雙渦脫流態(tài)下，兩種柱心間距的Ay2分別在振動鎖定區(qū)內(nèi)(L=2B)和區(qū)外(L=4B)出現(xiàn)明顯極值時的典型工況

見圖12，L=2B和4B的上游方柱均未見旋渦脫落，下游方柱的大幅橫流向振動主要由背風(fēng)面卷起的強烈旋渦引起；此外，下游方柱的大幅橫流向振動導(dǎo)致旋渦的橫流相間距變大，因而形成了模態(tài)Ⅱ.

見圖13，L=2B時的Lmean急劇減小，致使下游方柱流場更易受上游方柱尾渦影響，造成下游方柱升力對其運動進行不穩(wěn)定的能量輸入，同時也導(dǎo)致下游方柱的尾渦極不穩(wěn)定(模態(tài)Ⅴ).當(dāng)L=4B時，由于兩方柱間距較大，上游方柱尾渦會撞擊下游方柱迎風(fēng)面，并一分為二，其中一個旋渦的強度隨時間流逝而減弱(t1時刻)，而另一旋渦則和下游方柱脫落的渦發(fā)生融合(t1和t2時刻)，從而增強了升力對下游方柱的能量輸入，故而促進橫流向振動.由于結(jié)構(gòu)與流場的耦合作用，大幅橫流向振動拉大了旋渦的橫流向間距，最終形成模態(tài)Ⅳ.

圖12 剪切層再附流態(tài)內(nèi)典型工況的瞬時渦量圖

圖13 雙渦脫流態(tài)內(nèi)典型工況的瞬時渦量圖

4 結(jié) 論

通過對兩種典型柱心間距的串列雙方柱渦激振動進行數(shù)值研究，得到如下結(jié)論：

1)L=2B時，上、下游方柱均在雙渦脫流態(tài)內(nèi)發(fā)生“弱鎖定”現(xiàn)象，即振動鎖定區(qū)內(nèi)的振動頻率鎖定值fy/fn遠小于1.0.L=4B時，僅下游方柱在剪切層再附流態(tài)內(nèi)發(fā)生“鎖定”現(xiàn)象.

2)不同柱心間距雙方柱均以橫流向為主.L=2B時的上游方柱橫流向最大振幅大于下游方柱，且均出現(xiàn)在振動鎖定區(qū)內(nèi).L=4B時的上游方柱橫流向最大振幅遠小于下游方柱，且均出現(xiàn)在振動鎖定區(qū)外.各柱心間距下游方柱在振動鎖定區(qū)內(nèi)、外均有大幅橫流向振動.

3)L=2B時的上游方柱平均阻力系數(shù)在振動鎖定區(qū)內(nèi)出現(xiàn)激增，而下游方柱平均阻力系數(shù)在振動鎖定區(qū)內(nèi)為負值，導(dǎo)致兩方柱平均柱心間距急劇減小.各柱心間距的下游方柱在振動鎖定區(qū)內(nèi)、外均會有較大的脈動升力.

4)不同柱心間距的雙方柱隨折減速度變化共出現(xiàn)了2種柱間流態(tài)，即剪切層再附流態(tài)和雙渦脫流態(tài)，還出現(xiàn)了6種尾流模態(tài).其中，L=2B和4B時的下游方柱均在雙渦脫流態(tài)內(nèi)出現(xiàn)橫流向最大振幅，其中，L=2B時平均柱心間距的急劇減小致使上游方柱尾流對下游方柱的影響更為復(fù)雜，擾亂了下游方柱升力對其運動的能量輸入，使得下游方柱尾流中的旋渦極不穩(wěn)定；L=4B時上、下游方柱脫落的旋渦相融合，增強了升力對下游方柱的能量輸入，增大了橫流向振幅，導(dǎo)致旋渦的橫流向間距變大，也使得下游方柱的尾流為雙渦脫流態(tài)內(nèi)的平行渦街模態(tài).