基于浸入邊界算法的振動(dòng)二維矩形柱繞流模擬研究

楊 青， 曹曙陽(yáng)， 周軍文

(1. 安徽建筑大學(xué) 土木工程學(xué)院，合肥 230601; 2. 同濟(jì)大學(xué) 土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092;3. 常州工學(xué)院 常州市建設(shè)工程結(jié)構(gòu)與材料性能研究重點(diǎn)實(shí)驗(yàn)室，江蘇 常州 213032)

振動(dòng)鈍體繞流是一類具有廣泛工程應(yīng)用背景的流固耦合問(wèn)題，其普遍存在于大跨橋梁結(jié)構(gòu)、高層建筑以及風(fēng)力發(fā)電設(shè)備中[1-3]。柱體作為上述各類結(jié)構(gòu)中常見的鈍體斷面形式，更是成為此類基礎(chǔ)研究領(lǐng)域的重點(diǎn)。

早期振動(dòng)柱體研究多借助于風(fēng)洞試驗(yàn)完成[4-5]，近年來(lái)，隨著計(jì)算理論和硬件設(shè)備的發(fā)展，越來(lái)越多的學(xué)者傾向于利用CFD技術(shù)研究此類繞流問(wèn)題。

Singh等[6]利用有限體積法模擬雷諾數(shù)Re=100和150時(shí)，不同強(qiáng)迫振動(dòng)頻率下的方柱繞流，詳細(xì)分析了方柱鎖定區(qū)內(nèi)的氣動(dòng)系數(shù)及流場(chǎng)特征；Placzek等[7]采用有限體積方法模擬出Re=100時(shí)強(qiáng)迫振動(dòng)的二維圓柱繞流，通過(guò)柱體升力系數(shù)隨振動(dòng)頻率和振幅的變化曲線分析柱后尾流特征；國(guó)內(nèi)徐楓、歐進(jìn)萍[8]則借助Fluent研究Re=150時(shí)不同截面形狀對(duì)柱體振動(dòng)形式的影響；丁林等[9]進(jìn)一步在高雷諾數(shù)范圍(30 000≤Re≤110 000)展開不同截面形式柱體流致振動(dòng)問(wèn)題的研究，發(fā)現(xiàn)柱體流致振動(dòng)振幅和頻率與尾跡旋渦形態(tài)緊密相關(guān)；劉志文等[10]采用“剛性運(yùn)動(dòng)區(qū)域+動(dòng)網(wǎng)格區(qū)域+靜止網(wǎng)格區(qū)域”的思路建立網(wǎng)格，利用Fluent對(duì)寬高比為4的矩形斷面渦激振動(dòng)響應(yīng)進(jìn)行了數(shù)值模擬，其研究結(jié)果與已有文獻(xiàn)風(fēng)洞試驗(yàn)結(jié)果的吻合性，證明了Newmark-β方法結(jié)合動(dòng)網(wǎng)格技術(shù)的可行性。

然而，傳統(tǒng)CFD技術(shù)在求解振動(dòng)鈍體繞流問(wèn)題時(shí)，多需借助動(dòng)網(wǎng)格技術(shù)，計(jì)算耗費(fèi)較高。相比之下，浸入邊界方法(Immersed Boundary Method，IBM)利用力源擬化邊界的思想，直接可在靜止網(wǎng)格體系中實(shí)現(xiàn)振動(dòng)鈍體的邊界再現(xiàn)，略去網(wǎng)格再生步驟，能夠顯著減輕計(jì)算負(fù)擔(dān)，且保證較好的計(jì)算精度。Luo等[11-12]從此類算法理論基礎(chǔ)入手，致力于提高此類方法在動(dòng)態(tài)鈍體數(shù)值求解的精度；Chern等[13]則利用IBM模擬圓柱在Re=100～300時(shí)兩種不同自由度(豎向及扭轉(zhuǎn))設(shè)置下的風(fēng)致振動(dòng)，研究發(fā)現(xiàn)：?jiǎn)我回Q向自由度工況時(shí)，圓柱在鎖定區(qū)初始時(shí)刻，振動(dòng)較大；雙向自由度工況時(shí)，圓柱尾部流場(chǎng)渦型由P型向2C型轉(zhuǎn)換；鐘國(guó)華等[14]也基于IBM研究了均勻流場(chǎng)中圓柱的振動(dòng)問(wèn)題，研究顯示：圓柱阻力在運(yùn)動(dòng)過(guò)程中一直存在振蕩，在最大振幅位置處，阻力與已有數(shù)據(jù)對(duì)比，呈現(xiàn)較大的數(shù)值差異。

綜上所述，振動(dòng)柱體繞流研究目前仍多集中于方柱、圓柱等基本外形，對(duì)于實(shí)際建筑結(jié)構(gòu)中同時(shí)出現(xiàn)“分離”和“再附”現(xiàn)象的鈍體斷面，則鮮有涉及。IBM雖然可以更高效地應(yīng)用此類結(jié)構(gòu)繞流的研究，但因其早期過(guò)多關(guān)注于算法改進(jìn)，現(xiàn)有文獻(xiàn)僅在輪機(jī)、活塞等工程簡(jiǎn)化研究問(wèn)題中作出嘗試[15-16]，少有針對(duì)土木結(jié)構(gòu)流固耦合繞流模擬的應(yīng)用。因此，系統(tǒng)地開展不同邊長(zhǎng)比下振動(dòng)柱體的繞流研究，具有較高的工程應(yīng)用價(jià)值。

本文基于IBM核心概念，編寫了數(shù)值計(jì)算模型，在Re=1×103時(shí)，展開不同邊長(zhǎng)比下豎向強(qiáng)迫振動(dòng)矩形柱的繞流研究，其中柱體運(yùn)動(dòng)軌跡為y(t)=Asinφ=Asin(2πfct)，式中y(t)為物體軌跡，A為振幅，fc為振動(dòng)頻率，t為時(shí)間。模擬中振幅恒定為柱體迎風(fēng)面高度的14%，以便同已有文獻(xiàn)數(shù)據(jù)結(jié)果對(duì)比。柱體邊長(zhǎng)比B/D=1，3，5，分別對(duì)應(yīng)柱體流動(dòng)完全分離階段及再附階段。通過(guò)渦脫鎖定區(qū)、流場(chǎng)特征、氣動(dòng)力系數(shù)三個(gè)方面的討論，研究不同邊長(zhǎng)比下振動(dòng)柱體的繞流，豐富IBM在工程領(lǐng)域中的基礎(chǔ)研究。

1 數(shù)值模型

本文數(shù)值模型以IBM為基礎(chǔ)，結(jié)合雙線性插值(Bi-linear interpolation method)和內(nèi)部流體處理方法(Internal flow treatment)，以實(shí)現(xiàn)不同邊長(zhǎng)比柱體的強(qiáng)迫振動(dòng)繞流模擬。下文即對(duì)對(duì)該數(shù)值模型展開詳細(xì)介紹。

1.1 IBM控制方程

IBM利用流體與物體邊界相互作用的力學(xué)概念，將物體邊界簡(jiǎn)化為作用于邊界點(diǎn)上的力源，并引入至Navier-Stokes方程中參加流場(chǎng)迭代計(jì)算，進(jìn)而演化出邊界的幾何形狀和物理特性(式(1)～(2))：

(1)

(2)

式中：ui、uj為流體速度；ρ為流體密度；P為壓力；ν為流體運(yùn)動(dòng)黏性系數(shù)；Fi即為邊界力源，施加于圖1中邊界點(diǎn)上。

圖1 浸入邊界示意圖Fig.1 Illustration of immersed boundary

邊界力源Fi采用反饋構(gòu)造形式，其力源概念更加清晰，反饋循環(huán)的應(yīng)用也使其能夠?qū)崿F(xiàn)較高的邊界數(shù)值精度，計(jì)算式為

式中：αf和βf為反饋系數(shù)項(xiàng)，其取值直接關(guān)系到邊界擬化效率；ui(sur)為邊界控制點(diǎn)處流體速度，由流體控制方程計(jì)算得出；Vdesired為物體邊界設(shè)定速度，靜止?fàn)顟B(tài)時(shí)其值為零，若物體處于強(qiáng)迫振動(dòng)狀態(tài)時(shí)，Vdesired可由振動(dòng)方程計(jì)算得出；ui(sur)-Vdesired為邊界處流體速度與邊界設(shè)定速度的差異，表達(dá)式(3)正是基于此差異，通過(guò)反饋系數(shù)回饋出邊界力源使兩者速度接近，以滿足無(wú)滑移剛性邊界條件。內(nèi)部流體則選用自由發(fā)展方法展開處理。

1.2 振動(dòng)邊界實(shí)現(xiàn)技術(shù)

IBM模型的計(jì)算核心在于獲取實(shí)時(shí)邊界控制點(diǎn)處流體速度ui(sur)，用以構(gòu)造邊界力源，進(jìn)而迭代出鈍體外形。本次模擬對(duì)象為振動(dòng)邊界，且計(jì)算網(wǎng)格采用交錯(cuò)網(wǎng)格，邊界控制點(diǎn)很難與速度網(wǎng)格點(diǎn)重合，如圖2所示。

圖2 振動(dòng)邊界示意圖Fig.2 Illustration of oscillating boundary

圖2中振動(dòng)狀態(tài)下的邊界控制點(diǎn)E位置隨時(shí)間發(fā)生改變(Etn—Etn+1)，此時(shí)邊界速度不能由流體控制方程直接求得。因此，為保證邊界控制點(diǎn)速度的實(shí)時(shí)獲取，需借助本節(jié)開始所介紹的雙線性插值方法完成周圍速度網(wǎng)格定義點(diǎn)與邊界控制點(diǎn)之間的數(shù)據(jù)傳遞。

雙線性插值方法，其過(guò)程分為兩個(gè)步驟：速度獲取和力源分配。

圖3 雙線性插值方法示意圖Fig.3 Illustration of bilinear interpolation method

其中速度獲取如圖3(a)所示，邊界控制點(diǎn)速度u(E)可由周圍速度定義網(wǎng)格點(diǎn)插值得出，具體插值公式為

(4)

式中：ui,j為插值所需速度網(wǎng)格定義點(diǎn)；φi,j為在插值過(guò)程中起到關(guān)鍵作用的雙線性函數(shù)，其一般形式為

(5)

式中:Δ(xi)和Δ(yi)為二維計(jì)算網(wǎng)格間距；(xE,yE) 為邊界控制點(diǎn)E所處坐標(biāo); (xi,yi) 則表示為各插值點(diǎn)坐標(biāo)值。

力源分配則是指力源在邊界點(diǎn)構(gòu)造完成后再分配到周圍速度定義網(wǎng)格點(diǎn)上參與計(jì)算的過(guò)程。圖3(b)為該步驟示意圖。

F(ui,j)=φi,jF(uE)

(6)

式中:F(uE)為邊界控制點(diǎn)力源項(xiàng)；F(ui,j) 為插值網(wǎng)格點(diǎn)處力源值;φi,j為式(5)所示，表示雙線性函數(shù)。

具體至本文振動(dòng)矩形柱繞流模擬，首先設(shè)計(jì)出定位模塊，根據(jù)物體振動(dòng)軌跡計(jì)算出某一瞬時(shí)邊界點(diǎn)的運(yùn)動(dòng)坐標(biāo)，其次依據(jù)插值方法迅速確定邊界點(diǎn)周圍的速度定義點(diǎn)，利用邊界插值函數(shù)獲得邊界點(diǎn)處流動(dòng)參數(shù)以構(gòu)造邊界力源，最后再將邊界力源分配到周圍速度網(wǎng)格定義點(diǎn)上，以參與該瞬時(shí)流場(chǎng)的整體計(jì)算，進(jìn)而實(shí)現(xiàn)對(duì)運(yùn)動(dòng)邊界的擬化。圖4為本文振動(dòng)邊界實(shí)現(xiàn)技術(shù)流程圖。

圖4 振動(dòng)邊界實(shí)現(xiàn)技術(shù)流程圖Fig.4 Sketch of realization technique for moving Boundary

2 繞流計(jì)算域說(shuō)明

本次振動(dòng)矩形柱繞流均采用二維非均勻計(jì)算網(wǎng)格。流體從計(jì)算域左端施加，流動(dòng)從左至右。邊界條件統(tǒng)一設(shè)定為：入口邊界條件u=1.0，v=0；出口邊界條件采用對(duì)流邊界條件：?u/?t+um▽u=0；即流動(dòng)變量除壓力項(xiàng)外其余的邊界梯度值為零。上下邊界條件設(shè)定為無(wú)摩擦邊界：?u/?y=0; ?P/?y=0;v=0。力源反饋系數(shù)取值為αf=-1.6×105；βf=-60。

柱體繞流計(jì)算域中，特征尺寸取為柱體迎風(fēng)面高度D，計(jì)算區(qū)域中入口距方柱中心距離為10D，出口邊界則設(shè)置為31D。上下計(jì)算域邊界長(zhǎng)度為H，D/H設(shè)定為1/28?？紤]到方柱運(yùn)動(dòng)軌跡為豎向強(qiáng)迫振動(dòng)，為增強(qiáng)柱體表面流動(dòng)細(xì)節(jié)的捕捉，在其寬度B范圍內(nèi)，沿零時(shí)刻柱體上下表面各設(shè)置尺寸為D×B的均勻網(wǎng)格區(qū)域，格子大小統(tǒng)一選為0.01D。超出此區(qū)域，則設(shè)置拉伸網(wǎng)格，網(wǎng)格總數(shù)為412×242；時(shí)間間隔Δt為1.5×10-4，計(jì)算域布置如圖5所示。

圖5 振動(dòng)方柱計(jì)算域示意圖Fig.5 Illustration of calculation domain for oscillating rectangular

3 數(shù)值驗(yàn)證

為驗(yàn)證本文程序可靠性，首先利用該動(dòng)態(tài)程序模擬Re=103，A=0，邊長(zhǎng)比B/D=1的靜止方柱繞流。

圖6 方柱時(shí)均化表面壓力分布曲線Fig.6 Distribution curve of mean pressure coefficient Cp for square cylinder at Re=103

圖6為Re=103時(shí)方柱表面壓力分布。圖中可以看出柱體前緣表面壓力分布與Cao 等[17]數(shù)據(jù)基本吻合，方柱其余表面壓力分布特征亦吻合較好，但其值略大。原因可能存在于本文模擬為二維模擬，模擬維度的不同導(dǎo)致了流動(dòng)分離時(shí)流場(chǎng)結(jié)構(gòu)的差異性變化，進(jìn)而影響到方柱表面壓力分布。

圖7 方柱尾流區(qū)時(shí)均流向速度分布Fig.7 Mean stream velocity in wake region of square

圖7為方柱尾流區(qū)平均流向速度分布，U為入口參考流速。圖7(a)中x/D=0處方柱平均流向速度分布，由于其在鈍體尖銳邊角的流動(dòng)分離會(huì)而造成的分離層內(nèi)壓力下降，導(dǎo)致尾流區(qū)流體的動(dòng)量虧損，尾流區(qū)內(nèi)速度減小，在偏于中心線的部分存在速度最小值；圖7(b)中x/D=1處，由于已經(jīng)遠(yuǎn)離方柱，此效應(yīng)減弱，速度經(jīng)過(guò)展現(xiàn)延伸后，趨于回歸到來(lái)流風(fēng)速。方柱此時(shí)模擬所得St為0.133，與Cao等所得結(jié)果0.129基本吻合。

在此基礎(chǔ)上，展開模擬Re=103，A=0，邊長(zhǎng)比B/D=5的方柱繞流。

圖8 矩形柱渦量圖(B/D=5)Fig.8 Vortex contour of rectangular (B/D=5)

圖8為邊長(zhǎng)比B/D=5的矩形柱渦量圖。此時(shí)，矩形柱繞流已完全進(jìn)入沖擊剪切層階段，柱體側(cè)面附著渦形狀不再穩(wěn)定，渦脫結(jié)構(gòu)更加復(fù)雜。柱后流場(chǎng)結(jié)構(gòu)完全由柱體前緣流動(dòng)分離所決定。本文模擬所得流場(chǎng)特征圖符合Nakamura等[18]的研究描述。

進(jìn)一步將本次模擬所得不同邊長(zhǎng)比矩形柱斯托羅哈數(shù)St(D)(Re=103)與Nakamura等對(duì)比，如表1所示，可以看出，本文數(shù)值結(jié)果吻合較好。

表1 不同邊長(zhǎng)比矩形柱St(Re=103)

上述各類計(jì)算數(shù)據(jù)的吻合驗(yàn)證了本文IBM算法程序的可行性和準(zhǔn)確性。

4 振動(dòng)柱體繞流模擬

振動(dòng)矩形柱繞流計(jì)算中，關(guān)鍵無(wú)量綱參數(shù)包含：Stc=fcD/U，Stv=fvD/U，Stn=fvnD/U；其中fc為強(qiáng)迫振動(dòng)頻率，fv為振動(dòng)柱體渦脫頻率，fvn則代表靜止渦脫頻率；D為柱體高度；U為來(lái)流速度。

本節(jié)著手展開Re=1×103，振動(dòng)頻率Stc=0～0.5、振幅A=0.14D的不同邊長(zhǎng)比(B/D=1，3，5)強(qiáng)迫振動(dòng)柱體的繞流模擬，分析柱后尾流渦脫特征，并結(jié)合柱后近尾跡流場(chǎng)顯示、氣動(dòng)力系數(shù)以及氣動(dòng)穩(wěn)定性三個(gè)方面，探索邊長(zhǎng)比效應(yīng)對(duì)振動(dòng)柱體繞流氣動(dòng)特征的影響。

4.1 渦脫頻率分析

本次模擬參照Nakamura等的試驗(yàn)設(shè)置，亦在柱體后部1.5D處設(shè)置速度觀測(cè)點(diǎn)，用以收集流場(chǎng)速度時(shí)程數(shù)據(jù)，通過(guò)對(duì)其作頻譜分析，繪制出振動(dòng)矩形柱尾流渦脫隨振動(dòng)頻率變化的曲線。

圖9 方柱渦脫隨振動(dòng)頻率變化曲線(B/D=1)Fig.9 Variation of frequency ratio for squares with varying oscillating frequencies(B/D=1)

圖9為方柱渦脫隨振動(dòng)頻率變化曲線，通過(guò)速度頻譜分析得出其鎖定區(qū)長(zhǎng)度(Lock-in region length)為0.06～0.2。當(dāng)外加強(qiáng)迫振動(dòng)頻率Stc超過(guò)0.2時(shí)，柱體尾流渦脫頻率回歸靜止渦脫的現(xiàn)象十分突出，并一直保持，使得振動(dòng)頻率與渦脫頻率的比值(Stc/Stv)呈線性分布特征。本次模擬主要特征與Singh等得出的現(xiàn)象吻合，驗(yàn)證了本次模擬的可靠性。Taylor等[19]定義鎖定區(qū)中Stc=Stv=Stn的點(diǎn)，為共振鎖定點(diǎn)(Resonant point)，此時(shí)外加振動(dòng)頻率與靜止渦脫頻率完全保持一致，本次模擬中方柱(B/D=1)共振鎖定點(diǎn)為0.133。

圖10 B/D=3時(shí)柱體渦脫隨振動(dòng)頻率變化曲線Fig.10 Variation of frequency ratio for rectangular(B/D=3)with varying oscillating frequencies

圖10為B/D=3時(shí)矩形柱渦脫隨振動(dòng)頻率變化曲線。圖中可以看出：與振動(dòng)方柱對(duì)比，矩形柱(B/D=3)渦脫亦在外界振動(dòng)頻率Stc=0.06時(shí)進(jìn)入鎖定區(qū)，柱體尾流渦脫頻率與外界振動(dòng)頻率保持一致；當(dāng)振動(dòng)頻率提升至Stc=0.2時(shí)，柱后尾流旋渦脫鎖，振動(dòng)柱體渦脫頻率回歸到靜止時(shí)渦脫頻率(Stn=0.178)；但不同于方柱(B/D=1)，此后繼續(xù)增大強(qiáng)迫振動(dòng)頻率，至0.35時(shí)柱體再次進(jìn)入鎖定區(qū)，出現(xiàn)了呼和敖德同類試驗(yàn)所指出的多級(jí)鎖定現(xiàn)象[20]，并一直持續(xù)到Stc=0.45，此后柱后尾流渦脫再次脫鎖并接近靜止渦脫頻率。

圖11 觀測(cè)點(diǎn)速度頻譜圖(B/D=3)Fig.11 Frequency spectrum of observation point’s velocity(B/D=3)

圖11為矩形柱(B/D=3)尾流速度頻譜分析圖，圖11(a)可以看出Stc=0.25時(shí)，靜止渦脫頻率在尾流流場(chǎng)中起到控制作用，雖然振動(dòng)頻率響應(yīng)部分也有明顯的峰值，但仍小于靜止渦脫頻率部分，符合超越鎖定區(qū)后的脫鎖現(xiàn)象。圖11(b)Stc=0.35時(shí)，強(qiáng)迫振動(dòng)頻率已基本控制柱后尾流渦脫，靜止渦脫頻率部分對(duì)柱后尾流渦脫的影響遠(yuǎn)小于外界激勵(lì)部分。此前研究振動(dòng)柱體繞流的文獻(xiàn)，多集中于運(yùn)動(dòng)振幅增加對(duì)方柱渦脫頻率、氣動(dòng)力系數(shù)的影響，鮮有對(duì)于邊長(zhǎng)比效應(yīng)在此方面的研究。Nakamura等[21]雖然對(duì)不同邊長(zhǎng)比柱體的相位差、升力系數(shù)頻響虛部值展開了分析，也并未明確標(biāo)示出不同邊長(zhǎng)比下柱體渦脫鎖定區(qū)情況，但其研究中升力系數(shù)頻響虛部值兩次轉(zhuǎn)正的現(xiàn)象也從側(cè)面說(shuō)明柱后尾流渦脫多級(jí)鎖定的情況。

圖12 邊長(zhǎng)比B/D=5柱體渦脫隨振動(dòng)頻率變化曲線Fig.12 Variation of frequency ratio for rectangular(B/D=5)with varying oscillating frequencies

繼續(xù)增大邊長(zhǎng)比至B/D=5，頻譜分析其速度時(shí)程，得出渦脫隨振動(dòng)頻率變化曲線，如圖12所示。圖中可以看出，柱體在強(qiáng)迫振動(dòng)頻率Stc=0.15處才開始進(jìn)入鎖定區(qū)，此時(shí)柱體渦脫頻率完全由外界振動(dòng)頻率所控制；鎖定區(qū)長(zhǎng)度有所延長(zhǎng)，Stc=0.5時(shí)柱體渦脫頻率仍處于鎖定區(qū)；強(qiáng)迫振動(dòng)頻率升高至Stc=0.6時(shí)，柱體尾部旋渦已經(jīng)脫鎖，靜止渦脫頻率重新控制尾流流場(chǎng)。

4.2 流場(chǎng)顯示

流場(chǎng)顯示對(duì)于分析不同強(qiáng)迫振動(dòng)頻率下鈍體近尾跡流場(chǎng)特征十分重要，除此之外，流場(chǎng)顯示還可用以分析鈍體在不同振動(dòng)頻率下的氣動(dòng)力系數(shù)特征，揭示鎖定區(qū)尾流渦脫作用的流動(dòng)機(jī)理。

圖13為Stc=0.1時(shí)振動(dòng)相位分別為6π、8π的方柱流場(chǎng)渦量圖，結(jié)果顯示：柱后渦脫頻率與強(qiáng)迫振動(dòng)頻率保持一致，柱后渦脫基本呈現(xiàn)周期性再現(xiàn)。

圖13 不同相位振動(dòng)方柱渦量圖(Stc=0.1)Fig.13 Vortex contour of squares under different phase(Stc=0.1)

Yang等[22]在振動(dòng)柱體研究中發(fā)現(xiàn)，振動(dòng)物體渦的放出同運(yùn)動(dòng)方向密切相關(guān)，物體的運(yùn)動(dòng)會(huì)抑制同方向渦的產(chǎn)生。圖14為本次模擬所得矩形柱不同運(yùn)動(dòng)位置處的瞬時(shí)渦量圖，當(dāng)物體向上運(yùn)動(dòng)時(shí)(圖14(b))，柱體上側(cè)壁的流動(dòng)分離則得到明顯抑制，外形呈帶狀，下側(cè)壁流動(dòng)分離尺寸增大；當(dāng)物體向下運(yùn)動(dòng)時(shí)(圖14(c))柱體下方的流動(dòng)分離得到抑制，被壓成扁平帶狀，上側(cè)壁流動(dòng)分離與零位移位置相比，尺寸增大；零位移位置處(圖14(a))柱體上下表面不再明顯出現(xiàn)上述流動(dòng)抑制現(xiàn)象。此流動(dòng)特征隨著柱體振動(dòng)位置的不同而不斷的轉(zhuǎn)換，并影響到尾流區(qū)渦脫的放出，使渦脫頻率與振動(dòng)頻率逐漸保持一致。Yang等認(rèn)為這可能就是造成柱體渦脫頻率同振動(dòng)頻率鎖定的原因。

圖14 不同運(yùn)動(dòng)相位處方柱渦量圖(B/D=1，Stc=0.2)Fig.14 Vortex contour of squares at different oscillating phase (B/D=1，Stc=0.2)

圖15～16為B/D=3，5時(shí)，不同運(yùn)動(dòng)相位處矩形柱渦量圖。與方柱所得現(xiàn)象相同，柱體表面流動(dòng)分離特征隨著柱體振動(dòng)位移的不同而不斷的轉(zhuǎn)換，并影響到尾流區(qū)渦脫的放出，使渦脫頻率與振動(dòng)頻率逐漸保持一致；但由于柱體邊長(zhǎng)比的延長(zhǎng)，上述流動(dòng)特征在運(yùn)動(dòng)方向上的影響范圍亦會(huì)增大，這也可能是渦脫鎖定區(qū)隨邊長(zhǎng)比變化的原因。

圖15 不同運(yùn)動(dòng)相位處渦量圖(B/D=3，Stc=0.2)Fig.15 Vortex contour of rectangular under different oscillating phase(B/D=3，Stc=0.2)

圖16 不同運(yùn)動(dòng)相位處渦量圖(B/D=5，Stc=0.2)Fig.16 Vortex contour of rectangular under different oscillating phase(B/D=5，Stc=0.2)

4.3 氣動(dòng)系數(shù)

圖17為強(qiáng)迫振動(dòng)方柱表面壓力分布曲線，強(qiáng)迫振動(dòng)頻率Stc=0.04時(shí)，由于外界強(qiáng)迫振動(dòng)的施加，柱后渦脫加強(qiáng)，此時(shí)壓強(qiáng)分布表現(xiàn)為較強(qiáng)的吸力，；Stc=0.1時(shí)由于柱后尾流流場(chǎng)的轉(zhuǎn)換，柱后壓強(qiáng)吸力效應(yīng)減弱；當(dāng)Stc=0.15時(shí)，強(qiáng)迫振動(dòng)頻率靠近共振鎖定點(diǎn)，柱后壓強(qiáng)發(fā)生躍變，柱后吸力效應(yīng)突然增強(qiáng)；此后隨著外界振動(dòng)頻率逐漸增強(qiáng)，柱體背部表面壓強(qiáng)發(fā)生較大的減弱，Taylor和Vazza指出這是由于高頻振動(dòng)下，側(cè)壁渦脫和柱后渦脫相互干擾，共同作用的原因。

圖17 不同振動(dòng)頻率下方柱時(shí)均表面壓力分布曲線Fig.17 Distribution of mean surface pressure Cp for squares with varying oscillating frequencies

圖18 不同振動(dòng)頻率下矩形柱時(shí)均表面壓力分布曲線(B/D=3，5)Fig.18 Distribution of mean surface pressure Cp for rectangular with varying oscillating frequencies (B/D=3，5)

如圖18所示，與方柱相同，邊長(zhǎng)比B/D=3，5時(shí)矩形柱在不同強(qiáng)迫振動(dòng)工況下，迎風(fēng)面壓力分布曲線基本保持一致；柱體側(cè)面壓力分布曲線不再保持均勻分布，出現(xiàn)曲線特征，這可能是由于邊長(zhǎng)比增大，柱體側(cè)面發(fā)生流動(dòng)分離-再附而決定的；振動(dòng)頻率越高，柱體側(cè)面前緣負(fù)壓特征越明顯，其表面壓力曲線分布特征亦同時(shí)增強(qiáng)，表明此時(shí)柱體側(cè)面渦旋尺寸發(fā)生改變，流動(dòng)分離-再附區(qū)長(zhǎng)度減小。

圖19為矩形柱時(shí)均阻力系數(shù)隨振動(dòng)頻率變化曲線。不同邊長(zhǎng)比矩形柱阻力系數(shù)均在小振動(dòng)頻率處均呈現(xiàn)先上升再降低的趨勢(shì)；在共振鎖定點(diǎn)附近處，阻力系數(shù)達(dá)到峰值，結(jié)合圖17～18可看出，柱體背后表面吸力效應(yīng)均突然加強(qiáng)，此后阻力系數(shù)迅速降低。此現(xiàn)象發(fā)生原因，可能在于共振鎖定點(diǎn)之前，柱后旋渦脫落基本呈現(xiàn)平直脫落的特征，渦脫完全再附于柱體背面(圖13)，超出此振動(dòng)頻率時(shí)，方柱近尾跡流場(chǎng)特征再次發(fā)生改變，渦脫向柱側(cè)角退縮，再附于方柱背面的現(xiàn)象消失(圖14)。

圖19 不同邊長(zhǎng)比矩形柱阻力系數(shù)隨振動(dòng)頻率變化規(guī)律Fig.19 Variation of drag coefficients for different side ratio rectangular with varying oscillating frequencies

5 結(jié) 論

本文利用IBM算法模型，在靜止網(wǎng)格體系中完成豎向強(qiáng)迫振動(dòng)柱體的繞流模擬。文章通過(guò)設(shè)置不同的柱體邊長(zhǎng)比，從旋渦脫落、流場(chǎng)特征、氣動(dòng)系數(shù)等方面入手，探究邊長(zhǎng)比效應(yīng)對(duì)振動(dòng)柱體氣動(dòng)特性的影響：

(1)柱體鎖定區(qū)隨邊長(zhǎng)比的增大而發(fā)生變化，在柱體邊長(zhǎng)比為3時(shí)，鎖定區(qū)出現(xiàn)多級(jí)鎖定現(xiàn)象，邊長(zhǎng)比升高至5時(shí)，柱體鎖定區(qū)延長(zhǎng)。

(2)柱體運(yùn)動(dòng)會(huì)抑制同方向渦的產(chǎn)生，并影響到尾流區(qū)渦脫的放出；柱體邊長(zhǎng)比增大后，柱后尾流渦脫影響范圍變大，這也可能是不同邊長(zhǎng)比柱體鎖定區(qū)發(fā)生變化的原因。

(3)邊長(zhǎng)比B/D=3，5時(shí)，矩形柱側(cè)面壓力分布曲線不再保持均勻分布，出現(xiàn)曲線特征；振動(dòng)頻率越高，柱體側(cè)面前緣負(fù)壓特征越明顯，其表面壓力曲線分布特征亦同時(shí)增強(qiáng)，表明此時(shí)柱體側(cè)面渦旋尺寸發(fā)生改變，流動(dòng)分離-再附區(qū)長(zhǎng)度減小；振動(dòng)工況下，不同邊長(zhǎng)比柱體阻力系數(shù)均在接近自然渦脫頻率處(共振鎖定點(diǎn))發(fā)生峰值效應(yīng)。

在此研究基礎(chǔ)之上，后期將開展典型橋梁斷面的振動(dòng)繞流數(shù)值模擬，利用IBM數(shù)值程序的便捷性，增加計(jì)算風(fēng)工程對(duì)實(shí)際結(jié)構(gòu)的氣動(dòng)特性研究。