譚瀟玲，涂佳黃，雷 平，梁經(jīng)群，鄧旭輝

(1.湘潭大學(xué) 土木工程與力學(xué)學(xué)院，湖南 湘潭 411105；2.中冶長(zhǎng)天國(guó)際工程有限責(zé)任公司，長(zhǎng)沙 410205)

實(shí)際工程中，當(dāng)海洋平臺(tái)、海底管道、橋梁的斜拉索和橋塔等結(jié)構(gòu)物置于流場(chǎng)中時(shí)，周圍來流的速度往往會(huì)隨著時(shí)間和空間而發(fā)生改變，因此，僅限于對(duì)均勻來流作用下的流致振動(dòng)問題進(jìn)行研究并不足以反映工程當(dāng)中的實(shí)際情況，越來越多的學(xué)者將研究的重點(diǎn)放在了非均勻來流作用系統(tǒng)。按照來流性質(zhì)的不同，非均勻來流可分為周期性來流、Poiseuille管道流以及剪切來流。其中，對(duì)于平面剪切來流作用下的結(jié)構(gòu)系統(tǒng)，在任一高度截面圓柱體結(jié)構(gòu)兩側(cè)的來流速度由上至下呈現(xiàn)出遞減的趨勢(shì)，從而使得圓柱體結(jié)構(gòu)周圍渦量和流體速度分布呈現(xiàn)出較大的差異，速度差的存在引起的壓力差會(huì)影響結(jié)構(gòu)的表面受力以及運(yùn)動(dòng)形式。

通過對(duì)平面剪切來流作用下鈍體繞流問題進(jìn)行研究，發(fā)現(xiàn)來流剪切率的變化對(duì)結(jié)構(gòu)表面的壓力系數(shù)、升阻力系數(shù)、渦脫落頻率、邊界層分離點(diǎn)的位置以及尾流渦結(jié)構(gòu)等均有較大影響。Kang[1]采用浸泡邊界元法對(duì)低雷諾數(shù)工況下(50≤Re≤160)，單圓柱體結(jié)構(gòu)繞流特性進(jìn)行了數(shù)值模擬。研究發(fā)現(xiàn)隨著剪切率的增加(0≤K≤0.2)，結(jié)構(gòu)的渦脫落頻率和平均升力系數(shù)保持不變或者小幅下降，而平均阻力系數(shù)則明顯增大。但是，Sumner等[2]通過實(shí)驗(yàn)研究發(fā)現(xiàn)：亞臨界雷諾數(shù)(4×104≤Re≤9×104)與中低剪切率工況下(0.02≤K≤0.07)，單圓柱體結(jié)構(gòu)的渦脫落頻率和平均壓力系數(shù)隨剪切率的增加而增大，但平均阻力系數(shù)則減小。另外，結(jié)構(gòu)表面壓力的不對(duì)稱分布使得尾流渦結(jié)構(gòu)的寬度變窄。王偉等[3]對(duì)Re=100工況下，剪切率的變化(k≤0.5)對(duì)單方柱系統(tǒng)繞流特性的影響進(jìn)行了分析。研究結(jié)果表明，隨著剪切率的增加，尾流區(qū)域出現(xiàn)的卡門渦街逐漸消失，駐點(diǎn)位置會(huì)由最初向高速側(cè)分離點(diǎn)移動(dòng)后逐漸轉(zhuǎn)變?yōu)椴辉侔l(fā)生變化，而升阻力系數(shù)和斯特勞哈爾數(shù)的變化趨勢(shì)在不同的剪切率范圍內(nèi)會(huì)呈現(xiàn)出較大的差異。水慶象等[4]對(duì)不同剪切率工況下(k=0～0.4)的串列雙方柱系統(tǒng)繞流問題進(jìn)行了數(shù)值模擬，研究發(fā)現(xiàn)：剪切來流的作用會(huì)抑制結(jié)構(gòu)的渦脫落，駐點(diǎn)的位置隨剪切率的變化趨勢(shì)類似于單方柱工況。同時(shí)，剪切率的變化會(huì)改變串列雙方柱間隙流區(qū)域水平流速對(duì)結(jié)構(gòu)渦脫落的影響。

對(duì)于平面剪切來流作用下的渦激振動(dòng)問題，通過對(duì)以往的研究成果進(jìn)行總結(jié)，發(fā)現(xiàn)隨著剪切率的增大，單圓柱體結(jié)構(gòu)的渦激振動(dòng)響應(yīng)主要有以下三個(gè)方面的特征：①隨著剪切率的持續(xù)增大，結(jié)構(gòu)的順流向最大振幅值最終會(huì)遠(yuǎn)遠(yuǎn)大于橫流向最大振幅值[5]；②結(jié)構(gòu)鎖定區(qū)間的范圍會(huì)隨著剪切率的增加而擴(kuò)大[6]；③在流體從均勻來流向平面剪切來流轉(zhuǎn)變的過程中，隨結(jié)構(gòu)振動(dòng)響應(yīng)的變化，“8”字形運(yùn)動(dòng)軌跡的上半部分會(huì)逐漸消失，最終呈現(xiàn)為“雨滴”形運(yùn)動(dòng)軌跡[6-8]。而對(duì)于串列布置多柱體結(jié)構(gòu)系統(tǒng)，隨著流場(chǎng)中結(jié)構(gòu)數(shù)量的增加，振動(dòng)響應(yīng)會(huì)更加復(fù)雜，這主要是因?yàn)槌藴u激振動(dòng)之外，系統(tǒng)的振動(dòng)機(jī)制還有可能出現(xiàn)尾激振動(dòng)。在尾激振動(dòng)系統(tǒng)當(dāng)中，結(jié)構(gòu)的振動(dòng)響應(yīng)主要受到上游結(jié)構(gòu)尾流的影響，兩種振動(dòng)機(jī)理的共同作用會(huì)引起結(jié)構(gòu)更為劇烈的振動(dòng)響應(yīng)。目前，對(duì)于串列布置多柱體結(jié)構(gòu)群流致振動(dòng)響應(yīng)的研究多集中于均勻來流工況。陳威霖等[9-10]對(duì)小間距比工況下串列雙柱系統(tǒng)的渦激振動(dòng)問題進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)不同于單柱系統(tǒng)，上游圓柱體的振動(dòng)響應(yīng)僅存在上端分支。另外，捕捉到T+S和P+S兩種尾流模態(tài)?；诟褡硬柶澛椒?，Haider等[11]對(duì)Re=100時(shí)，不同間距比工況下串列雙圓柱系統(tǒng)的流致振動(dòng)響應(yīng)進(jìn)行了數(shù)值模擬，研究表明：隨著間距比的增加，上游圓柱體受到下游結(jié)構(gòu)臨近效應(yīng)的影響，振幅會(huì)呈現(xiàn)出“先大于、后小于、再接近于單柱體工況”的變化趨勢(shì)。而下游圓柱在上游結(jié)構(gòu)尾流渦干擾的影響下，振幅值在大多數(shù)間距比工況下均大于單柱體結(jié)構(gòu)。Gao等[12]也考慮了間距比變化對(duì)串列三圓柱系統(tǒng)的影響，發(fā)現(xiàn)當(dāng)L/D=2時(shí)，上游圓柱的振幅值為單柱體工況的兩倍。

本文基于四步半隱式特征線分裂算子有限元方法，對(duì)平面剪切來流作用下串列布置三圓柱體結(jié)構(gòu)單自由度流致振動(dòng)問題進(jìn)行了數(shù)值模擬，重點(diǎn)分析剪切率、雷諾數(shù)與折減速度的變化對(duì)結(jié)構(gòu)橫流向動(dòng)力響應(yīng)、流體力系數(shù)、相位階段變化、頻率特性與尾流模態(tài)等特性的影響，在實(shí)際工程中具有一定的參考價(jià)值。

1 數(shù)值方法及驗(yàn)證

1.1 控制方程

基于任意-拉格朗日歐拉(ALE)描述下的不可壓縮黏性流體的N-S方程，其無量綱形表達(dá)式為

(1)

(2)

式中：t為時(shí)間；ui和p分別為流體速度和壓力；xi和xj分別為i,j方向的空間坐標(biāo)；cj=uj-wj，其中cj為流體相對(duì)于網(wǎng)格移動(dòng)速度的對(duì)流速度，wj為網(wǎng)格移動(dòng)速度；雷諾數(shù)Re=ρU0D/v，其中U0和D分別為流場(chǎng)特征速度和特征尺寸，v為流體動(dòng)力黏性系數(shù)，ρ為流體密度。

彈性支撐圓柱結(jié)構(gòu)運(yùn)動(dòng)體系可假設(shè)為質(zhì)量-阻尼-彈簧系統(tǒng)，考慮雙自由度運(yùn)動(dòng)的結(jié)構(gòu)控制方程量綱歸一化形式為

(3)

(4)

1.2 網(wǎng)格更新

為了防止網(wǎng)格失效，求解Laplace方程得到每個(gè)動(dòng)網(wǎng)格節(jié)點(diǎn)位移，其邊值問題表達(dá)式如下

(5)

式中：τ是網(wǎng)格形變控制參數(shù)；Si是網(wǎng)格節(jié)點(diǎn)i方向位移；gi為運(yùn)動(dòng)邊界節(jié)點(diǎn)位移；Гm是網(wǎng)格動(dòng)邊界；Гf為網(wǎng)格固定邊界。

1.3 弱耦合分區(qū)算法

采用基于四步半隱式特征線分裂算子的弱耦合分區(qū)算法，運(yùn)用Fortran語言編制相應(yīng)的數(shù)值計(jì)算程序求解流固耦合問題，詳細(xì)的計(jì)算框架可參考文獻(xiàn)[13]，在此不再贅述。該算法在本課題組以前的工作中已多次驗(yàn)證其能準(zhǔn)確地求解鈍體結(jié)構(gòu)群渦激振動(dòng)問題[14]。

1.4 算例驗(yàn)證

將本文計(jì)算結(jié)果與文獻(xiàn)[15]中均勻來流作用下串列布置三圓柱體結(jié)構(gòu)的流致振動(dòng)響應(yīng)進(jìn)行對(duì)比來驗(yàn)證本文數(shù)值方法的可靠性和適用性。計(jì)算模型如圖1所示。三圓柱體均僅有橫向自由度，圓柱表面速度為u=?X/?t,v=?Y/?t。無量綱參數(shù)設(shè)置如下：Re=100，L/D=5，mr=2.0，ξ=0.0。計(jì)算域尺寸為[-30D，60D]×[-10D，10D]，上游圓柱(upstream cylinder,UC)、中游圓柱(midstream cylinder,MC)和下游圓柱(downstream cylinder,DC)的圓心位置分別為(-5D，0)，(0，0)與(5D，0)，邊界條件設(shè)置如下：入口為速度入口：ux=1，uy=0；出口為壓力出口：p=0；上下邊界為自由滑移邊界：?ux/?y=0,uy=0；結(jié)構(gòu)表面為無滑移邊界：ux=uy=0。將本文計(jì)算得出的串列三圓柱體結(jié)構(gòu)的橫流向最大振幅值與文獻(xiàn)結(jié)果進(jìn)行對(duì)比，如圖2所示，通過對(duì)比發(fā)現(xiàn)，計(jì)算結(jié)果與文獻(xiàn)結(jié)果的高吻合度驗(yàn)證了本文數(shù)值方法的可靠性和適用性。

圖1 驗(yàn)證算例計(jì)算模型示意圖Fig.1 Model schematic of validation case.

圖2 串列三圓柱體橫流向最大振幅本文結(jié)果與文獻(xiàn)結(jié)果的對(duì)比Fig.2 Comparison of the present numerical results of the maximum value of the transverse vibration amplitude(Ymax/D)with existing results

2 問題描述

2.1 計(jì)算模型

平面剪切來流作用下串列布置三圓柱體結(jié)構(gòu)流致振動(dòng)響應(yīng)的計(jì)算模型以及邊界條件，如圖3所示。三圓柱體僅具有橫流向自由度，無量綱參數(shù)設(shè)置如下：為了深入分析層流范圍內(nèi)(Re<160～210)結(jié)構(gòu)的動(dòng)力響應(yīng)以及流體與結(jié)構(gòu)之間的互擾機(jī)理，取Re=100和Re=160，由于流體狀態(tài)較為穩(wěn)定，可以采用二維模擬進(jìn)行求解。入口為速度入口：ux=Uc+ky(Uc=1.0),uy=0。由于阻塞率B=0.05，為了保證入口處來流方向一致，剪切率k=0.00，k=0.05，k=0.10。間距比選取L/D=5.5，重點(diǎn)考慮上游結(jié)構(gòu)尾流發(fā)展完全后對(duì)下游結(jié)構(gòu)振動(dòng)響應(yīng)的影響。其他參數(shù)為：折減速度Ur=3～21，mr=2.0，忽略結(jié)構(gòu)阻尼的影響，即ξ=0。UC、MC和DC三圓柱體結(jié)構(gòu)的圓心位置分別為(-5.5D，0)、(0，0)和(5.5D，0)。計(jì)算域大小與其他邊界條件的設(shè)置與1.4節(jié)中完全一致。

圖3 計(jì)算模型示意圖Fig.3 The schematic diagram of the computational model

2.2 計(jì)算參數(shù)獨(dú)立性驗(yàn)證

本文計(jì)算域的網(wǎng)格劃分采用非結(jié)構(gòu)化的三角形網(wǎng)格，為了更準(zhǔn)確地刻畫出圓柱周圍及尾流區(qū)域范圍內(nèi)渦結(jié)構(gòu)的形成和發(fā)展，對(duì)這兩個(gè)區(qū)域進(jìn)行網(wǎng)格加密處理，如圖4所示。為了確保網(wǎng)格質(zhì)量，對(duì)于單自由度串列布置三圓柱系統(tǒng)，選取Re=100，k=0.1與Ur=7工況對(duì)網(wǎng)格數(shù)量的變化對(duì)結(jié)構(gòu)振動(dòng)響應(yīng)的影響進(jìn)行分析，從而判斷哪種類型的網(wǎng)格更適用于求解該問題。由表1可知：該工況下，與粗網(wǎng)格GI相比，隨著網(wǎng)格密度的增大，GII工況下三圓柱體結(jié)構(gòu)橫流向振幅(Ymax/D)的最大變化為4.40%，而在GIII工況下，最大變化則為4.89%。由于網(wǎng)格數(shù)量從GII增大到GIII時(shí)，振幅值變化相對(duì)較小，且若采用網(wǎng)格GIII進(jìn)行計(jì)算，雖能獲得最高的計(jì)算精度，但是計(jì)算總耗時(shí)會(huì)大幅增大，計(jì)算效率降低。綜上所述，本文計(jì)算采用的網(wǎng)格密度和分布情況均與GII類似。另一方面，選取該算例，并采用GII網(wǎng)格，分析三種不同時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響，從而判斷哪種類型的時(shí)間步長(zhǎng)對(duì)求解該問題更加適用。由表2可知：與Δt=0.005工況相比，隨著時(shí)間步長(zhǎng)的減小，Δt=0.002工況下三圓柱體結(jié)構(gòu)橫流向振幅(Ymax/D)和渦脫落頻率(St)的最大變化分別為2.41%和4.65%，而在Δt=0.001工況下，最大變化則為2.61%和6.11%。因此，本文無量綱計(jì)算步長(zhǎng)選取為Δt=0.002。

圖4 網(wǎng)格劃分示意圖Fig.4 Schematic of the computational mesh

表1 網(wǎng)格獨(dú)立性驗(yàn)證：串列布置三圓柱單自由度計(jì)算結(jié)果(Re=100,k=0.10,Ur=7)Tab.1 Grid independence test:the results for the three tandem circular cylinders at Re=100,k=0.10,Ur=7

表2 時(shí)間步長(zhǎng)獨(dú)立性驗(yàn)證：串列布置三圓柱單自由度計(jì)算結(jié)果(Re=100,k=0.10,Ur=7)Tab.2 Time-step independence test:the results for the three tandem circular cylinders at Re=100,k=0.10,Ur=7

3 計(jì)算結(jié)果分析

3.1 動(dòng)力響應(yīng)

本節(jié)分析了Re=100和Re=160工況下，剪切率和折減速度的變化對(duì)串列布置三圓柱體結(jié)構(gòu)橫流向動(dòng)力響應(yīng)的影響，如圖5所示。

由圖5可知，不同工況下，三圓柱體結(jié)構(gòu)的振幅曲線均呈現(xiàn)出“先增加后減小再趨于平穩(wěn)”的變化趨勢(shì)。根據(jù)不同的折減速度范圍內(nèi)剪切率的變化對(duì)橫流向最大振幅值影響的差異，可將三圓柱體結(jié)構(gòu)的振幅曲線劃分為Ⅰ、Ⅱ、Ⅲ和Ⅳ四個(gè)區(qū)域。對(duì)于UC，剪切率的變化對(duì)橫流向最大振幅的影響普遍較小，在區(qū)域Ⅰ(3≤Ur<5)范圍內(nèi)，振幅值隨剪切率的增加小幅下降，相反，隨折減速度的增加線性增大。在區(qū)域Ⅱ～Ⅳ(Ur≥5)范圍內(nèi)，UC橫流向最大振幅值隨剪切率的增加小幅增大。在Ur=5工況下，振幅曲線取得最大值0.59(k=0.1)，隨折減速度的進(jìn)一步增大，振幅減小并逐漸趨于穩(wěn)定。對(duì)于MC，區(qū)域Ⅰ的范圍擴(kuò)大到3≤Ur<8。特別地，Ur=6工況下，剪切率的變化對(duì)橫流向最大振幅的影響明顯增強(qiáng)，k=0.0工況下的最大振幅值達(dá)到0.79，超過k=0.1工況(Ymax/D=0.32)振幅值的兩倍。與UC不同的是，在經(jīng)歷了區(qū)域Ⅱ范圍內(nèi)(8≤Ur<19)橫流向振幅隨剪切率的增加小幅增大后，MC振幅曲線在區(qū)域Ⅳ范圍內(nèi)(Ur≥21)會(huì)再次出現(xiàn)隨剪切率的增加而減小的情況，且兩種變化模式中間存在一個(gè)過渡區(qū)域Ⅲ(19≤Ur<21)，即振幅不再隨著剪切率線性變化。另外，隨折減速度的增加，MC橫流向最大振幅整體變化趨勢(shì)與UC基本一致，分別在Ur=7(k=0.00)和Ur=8(k=0.05，k=0.10)工況下達(dá)到最大值，最大為1.04(k=0.10)。值得注意的是，Ur=4工況下，結(jié)構(gòu)發(fā)生了“弱鎖定”現(xiàn)象，對(duì)應(yīng)的橫流向振幅出現(xiàn)了極大值。對(duì)于DC，根據(jù)不同折減速度范圍內(nèi)最大振幅值隨剪切率的變化差異所劃分出的四個(gè)區(qū)域的范圍進(jìn)一步發(fā)生改變，分別為3≤Ur<10(區(qū)域Ⅰ)、10≤Ur<20(區(qū)域Ⅱ)、20≤Ur<24(區(qū)域Ⅲ)和Ur≥24(區(qū)域Ⅳ)。與中上游兩圓柱體不同的是，DC的振幅曲線中，在區(qū)域Ⅰ范圍內(nèi)部分的折減速度工況下也發(fā)現(xiàn)了振幅隨剪切率的非線性變化現(xiàn)象。這主要是因?yàn)樵谥猩嫌螆A柱尾流的影響下，DC周圍旋渦分布會(huì)更加復(fù)雜，流場(chǎng)的復(fù)雜化使得剪切率對(duì)DC振動(dòng)響應(yīng)的影響發(fā)生了改變。不同剪切率工況下，DC振幅曲線的最大值為1.32(k=0.10,Ur=11)。

圖5 當(dāng)Re=100時(shí)，串列布置三圓柱體結(jié)構(gòu)橫流向最大振幅隨剪切率和折減速度的變化情況Fig.5 The variation of the maximum crossflow vibrating amplitudes with shear ratio and reduced velocity at Re=100

當(dāng)Re=160時(shí)，隨著雷諾數(shù)的增大，三圓柱體結(jié)構(gòu)周圍流場(chǎng)會(huì)逐漸出現(xiàn)三維效應(yīng)，不穩(wěn)定性明顯增強(qiáng)。流場(chǎng)的復(fù)雜化使得中下游兩圓柱體的橫流向最大振幅值隨剪切率的變化情況在不同的折減速度工況下具有較大的差異。與Re=100工況不同的是，三圓柱體結(jié)構(gòu)的最大振幅曲線大致可以劃分為隨剪切率減小(區(qū)域Ⅰ)、增大(區(qū)域Ⅱ)以及非線性變化(區(qū)域Ⅲ)三個(gè)不同的區(qū)域，且區(qū)域Ⅰ和Ⅱ范圍內(nèi)最大振幅值隨剪切率的非線性變化現(xiàn)象出現(xiàn)的頻率明顯增大，如圖6所示。另外，中上游兩圓柱體振幅曲線隨折減速度的整體變化趨勢(shì)與振幅峰值均接近于Re=100工況。不同的是，當(dāng)Re=160時(shí)，UC達(dá)到峰值以及開始穩(wěn)定所對(duì)應(yīng)的折減速度減小，意味著流場(chǎng)的變化使得UC發(fā)生共振所對(duì)應(yīng)的頻率值會(huì)增大。與僅受到來流影響的UC不同，中下游兩圓柱體周圍流場(chǎng)在旋渦脫落影響下的復(fù)雜程度遠(yuǎn)遠(yuǎn)大于來流，使得中下游兩圓柱體振幅難以較快地穩(wěn)定下來，尤其是k=0.10工況。當(dāng)中下游兩圓柱體的振幅曲線達(dá)到峰值后，k=0.10工況所對(duì)應(yīng)結(jié)構(gòu)的振幅值明顯大于其他兩個(gè)剪切率工況，說明雷諾數(shù)的增加會(huì)對(duì)大剪切率工況下結(jié)構(gòu)的動(dòng)力響應(yīng)有較大的影響。當(dāng)k=0.10時(shí)，DC分別達(dá)到峰值1.54(Ur=12)與1.17(Ur=18)，兩個(gè)區(qū)域中間過渡所對(duì)應(yīng)的極小值點(diǎn)的位置為Ur=17。

圖6 當(dāng)Re=160時(shí)，串列布置三圓柱體結(jié)構(gòu)橫流向最大振幅隨剪切率和折減速度的變化情況Fig.6 The variation of the maximum crossflow vibrating amplitudes with shear ratio and reduced velocity at Re=160

3.2 柱體表面流體力特性

圖7～圖10分別給出了Re=100和Re=160工況下，串列布置三圓柱體結(jié)構(gòu)表面阻力系數(shù)平均值和升力系數(shù)均方根值隨剪切率與折減速度的變化情況。

圖10 當(dāng)Re=160時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)升力系數(shù)均方根值隨折減速度的變化情況Fig.10 The variation of the rms value of the lift coefficient with shear ratio and reduced velocity when Re=160

由圖7可知，不同剪切率工況下，三圓柱體結(jié)構(gòu)的阻力系數(shù)平均值隨折減速度的變化趨勢(shì)與橫流向最大振幅曲線基本一致，均為先增加到最大值后逐漸減小并趨于穩(wěn)定，且取得峰值所對(duì)應(yīng)的折減速度工況也相當(dāng)接近。值得注意的是，當(dāng)結(jié)構(gòu)置于流場(chǎng)中時(shí)，會(huì)在來流的作用下產(chǎn)生順流向的推力，而對(duì)于串列布置多柱體結(jié)構(gòu)，受到來流影響越大的圓柱所產(chǎn)生的推力就會(huì)越大。從圖中可以看出，自由來流對(duì)三圓柱體結(jié)構(gòu)的作用大小是依次遞減的，當(dāng)Re=100時(shí)，中上游兩圓柱體的最小阻力系數(shù)平均值分別接近于1.00(UC)和0.40(MC)，會(huì)大于DC。當(dāng)Ur=5～6時(shí)，部分剪切率工況下，下游圓柱體的阻力系數(shù)平均值接近于零，這主要是因?yàn)橄掠螆A柱完全處于中上游兩圓柱體的尾流中而被屏蔽。另外，隨著剪切率的增加，下游圓柱的阻力系數(shù)平均值曲線所對(duì)應(yīng)的最大值分別為1.61(k=0.00)，1.21(k=0.05)和0.91(k=0.10)，說明剪切來流作用下流體速度的不均勻分布減小了結(jié)構(gòu)順流方向所受的流體力。由圖8可知，不同剪切率工況下，三圓柱體結(jié)構(gòu)的升力系數(shù)均方根值隨折減速度的變化曲線呈現(xiàn)出“峰+平臺(tái)”的變化趨勢(shì)，且中間存在一個(gè)過渡區(qū)域Ur=5～8。剪切率的變化對(duì)中上游兩圓柱體升力系數(shù)均方根值的最大值影響較小。而對(duì)于下游圓柱，Ur=4工況下，最大值隨著剪切率的增加明顯減小。在過渡區(qū)域內(nèi)，上游圓柱的升力系數(shù)均方根值接近于零，而中下游兩圓柱體正處于共振區(qū)域內(nèi)，升力系數(shù)的波動(dòng)性較大，結(jié)構(gòu)的振動(dòng)較不穩(wěn)定。剪切率的變化對(duì)“平臺(tái)”范圍內(nèi)上游圓柱升力系數(shù)均方根值的影響較小，而中下游兩圓柱體的升力系數(shù)均方根值則隨著剪切率的增加而減小，下游圓柱更為明顯，這與其振幅曲線的變化情況基本一致。這可能是因?yàn)檎蹨p速度較大時(shí)，結(jié)構(gòu)剛度較小，流體與結(jié)構(gòu)之間的耦合作用減弱，自由來流對(duì)下游圓柱體的影響增強(qiáng)，再加上圓柱表面渦量分布的差異，使得剪切來流作用下下游圓柱的升力系數(shù)均方根值會(huì)接近于上游圓柱穩(wěn)定后的數(shù)值。

圖7 當(dāng)Re=100時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)阻力系數(shù)平均值隨折減速度的變化情況Fig.7 The variation of the mean value of the drag coefficient with shear ratio and reducedvelocity when Re=100

圖8 當(dāng)Re=100時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)升力系數(shù)均方根值隨折減速度的變化情況Fig.8 The variation of the rms value of the lift coefficient with shear ratio and reduced velocity when Re=100

當(dāng)Re=160時(shí)，雷諾數(shù)的增大對(duì)上游圓柱體升阻力系數(shù)的影響較小，而對(duì)于中下游圓柱的影響較大，如圖9和10所示。由于來流不穩(wěn)定性增強(qiáng)，中游圓柱最小阻力系數(shù)平均值減小50%。而對(duì)于下游圓柱，當(dāng)Ur=4時(shí)，與Re=100工況不同，阻力系數(shù)平均值僅在k=0.10工況下出現(xiàn)極大值0.84，而在其他剪切率工況下則接近于零，這與振幅曲線的變化基本一致。另外，不同剪切率工況下，下游圓柱阻力系數(shù)平均值的最大值有所減小，特別是k=0.05工況，大約減小了20%。當(dāng)k=0.10,Ur>15時(shí)，下游圓柱的阻力系數(shù)平均值略微大于其他剪切率工況，這與結(jié)構(gòu)振動(dòng)的第二個(gè)大振幅區(qū)域?qū)?yīng)。另一方面，k=0.00和k=0.05工況下，下游圓柱的升力系數(shù)均方根值峰值的位置轉(zhuǎn)移到了過渡區(qū)域范圍內(nèi)且數(shù)值明顯減小，而k=0.10工況下則恰好相反，與低雷諾數(shù)工況相比，最大值增加了37%。另外，不同剪切率工況下，由于下游圓柱產(chǎn)生劇烈振動(dòng)所對(duì)應(yīng)的區(qū)域范圍擴(kuò)大，升力系數(shù)均方根值曲線過渡區(qū)域的范圍擴(kuò)大到Ur=5～15范圍內(nèi)。值得注意的是，平穩(wěn)狀態(tài)的出現(xiàn)均會(huì)遲于Re=100工況，且剪切率的變化對(duì)升力系數(shù)均方根值的顯著影響會(huì)發(fā)生在大折減速度工況。

圖9 當(dāng)Re=160時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)阻力系數(shù)平均值隨折減速度的變化情況Fig.9 The variation of the mean value of the drag coefficient with shear ratio and reduced velocity when Re=160

3.3 相位階段變化

圖11與圖12分別給出了Re=100和Re=160工況下，串列布置三圓柱體結(jié)構(gòu)體系橫流向位移和升力系數(shù)之間的相位差隨剪切率和折減速度的變化情況。柱體結(jié)構(gòu)的力和位移之間相位變化會(huì)影響柱體與流場(chǎng)之間能量的傳遞方式，從而影響結(jié)構(gòu)的振動(dòng)響應(yīng)[16]。

由圖可知，不同雷諾數(shù)和剪切率工況下，三圓柱體結(jié)構(gòu)橫流向位移與升力系數(shù)之間相位差的變化趨勢(shì)基本一致：當(dāng)折減速度較小時(shí)，相位差接近于0°，結(jié)構(gòu)會(huì)從流體當(dāng)中獲得能量，且隨著系統(tǒng)剛度逐漸減小，結(jié)構(gòu)受周圍流體的影響增強(qiáng)，從而引起共振現(xiàn)象的發(fā)生，激發(fā)出較大的振幅。隨著折減速度的增加，結(jié)構(gòu)的振動(dòng)模式發(fā)生轉(zhuǎn)變的過程中，相位差會(huì)逐漸轉(zhuǎn)變?yōu)?80°，發(fā)生轉(zhuǎn)變所對(duì)應(yīng)的工況稱為“相位開關(guān)”。當(dāng)Re=100時(shí)，不同剪切率工況下，三圓柱體結(jié)構(gòu)“相位開關(guān)”所對(duì)應(yīng)的位置均為Ur=6工況，流場(chǎng)的不穩(wěn)定性以及剪切率的變化會(huì)對(duì)相位的轉(zhuǎn)變過程產(chǎn)生較大的影響，中游圓柱體力和位移之間的相位差從0°完全轉(zhuǎn)變?yōu)?80°所對(duì)應(yīng)折減速度的范圍明顯擴(kuò)大，如圖11所示。而對(duì)于下游圓柱，影響則更為顯著，相位差隨折減速度的變化曲線波動(dòng)性較大，力和位移的相位關(guān)系在結(jié)構(gòu)不規(guī)則運(yùn)動(dòng)的過程中會(huì)發(fā)生顯著變化，使得大折減速度工況下，難以快速轉(zhuǎn)變?yōu)橥耆聪嗟臓顟B(tài)，且這種現(xiàn)象隨著剪切率的增大愈加明顯。當(dāng)Re=160時(shí)，上游圓柱以及均勻來流作用下下游圓柱的“相位開關(guān)”提前至Ur=5工況，如圖12所示。剪切率的變化對(duì)中游圓柱相位差的影響增強(qiáng)，Ur≥8范圍內(nèi)，相位差隨折減速度的變化曲線會(huì)呈現(xiàn)出下凹趨勢(shì)，且下凹的幅度以及范圍會(huì)隨剪切率的增加而增大。相反，與Re=100工況相比，下游圓柱的相位關(guān)系受剪切率的影響減小。

圖11 當(dāng)Re=100時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)橫流向位移與升力系數(shù)相位差隨折減速度的變化情況Fig.11 The variation of the phase difference between the crossflow displacement and the lift coefficient with shear ratio and reduced velocity when Re=100

圖12 當(dāng)Re=160時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)橫流向位移與升力系數(shù)相位差隨折減速度的變化情況Fig.12 The variation of the phase difference between the crossflow displacement and the lift coefficient with shear ratio and reduced velocity when Re=160

3.4 頻譜特性

圖13與圖14選取了Re=100和Re=160時(shí)，部分工況下中下游兩圓柱體橫流向位移能量譜密度(power spectral density,PSD)曲線的變化情況。在結(jié)構(gòu)橫流向位移的PSD曲線當(dāng)中，振動(dòng)頻率的階數(shù)越高，結(jié)構(gòu)的振動(dòng)響應(yīng)就會(huì)越復(fù)雜，主頻所對(duì)應(yīng)的能量越大，結(jié)構(gòu)的振動(dòng)越劇烈。

由圖13(a)可知，當(dāng)Re=100時(shí)，均勻來流作用下，中下游兩圓柱體振動(dòng)位移所對(duì)應(yīng)的PSD曲線隨折減速度的變化趨勢(shì)和3.1小節(jié)結(jié)構(gòu)的振動(dòng)響應(yīng)一一對(duì)應(yīng)。隨折減速度的增加，中游圓柱體的PSD曲線中主頻的大小基本一致，且振動(dòng)主頻所對(duì)應(yīng)的能量在共振區(qū)間內(nèi)明顯增大。而下游圓柱的主頻則隨Ur增加逐漸減小，且共振區(qū)間內(nèi)次頻的數(shù)量明顯增大。當(dāng)Ur=3時(shí)，結(jié)構(gòu)的動(dòng)力響應(yīng)較弱，PSD曲線中僅出現(xiàn)一個(gè)振動(dòng)頻率，且對(duì)應(yīng)的能量較小。隨Ur增加，次頻數(shù)量的增多以及主頻能量的增大使得結(jié)構(gòu)的振動(dòng)出現(xiàn)“調(diào)制周期”、不規(guī)則振動(dòng)等不同的模式，動(dòng)力響應(yīng)更加劇烈。當(dāng)中下游兩圓柱體的振動(dòng)響應(yīng)隨折減速度的變化曲線處于下降段時(shí)，PSD曲線中次頻的數(shù)量以及主頻所對(duì)應(yīng)的能量值也隨之減小。大振幅區(qū)域內(nèi)，下游圓柱主頻所對(duì)應(yīng)的能量值明顯大于中游圓柱，與兩圓柱體最大振幅值對(duì)應(yīng)，如圖5所示。剪切來流作用時(shí)，部分折減速度工況下，兩圓柱體均會(huì)出現(xiàn)一個(gè)數(shù)值接近于零的振動(dòng)頻率，且當(dāng)Ur=20和Ur=40時(shí)，該頻率所對(duì)應(yīng)的能量值相當(dāng)大，但由于頻率值較小，對(duì)于結(jié)構(gòu)振動(dòng)響應(yīng)的影響可以忽略不計(jì)，如圖13(b)所示。

圖13 當(dāng)Re=100時(shí)，不同剪切率工況下中下游兩圓柱體結(jié)構(gòu)橫流向位移能量譜密度隨折減速度的變化情況Fig.13 The variation of the power spectral density of the crossflow displacement with shear ratio and reduced velocity when Re=100

當(dāng)Re=160時(shí)，均勻來流作用下，除了Ur=12工況外，其他折減速度工況下中游圓柱所對(duì)應(yīng)的振動(dòng)主頻均略微大于Re=100工況，如圖14(a)所示。而當(dāng)Ur=12時(shí)，會(huì)出現(xiàn)兩個(gè)能量接近的頻率值。對(duì)于下游圓柱，Ur=12工況下結(jié)構(gòu)振動(dòng)主頻所對(duì)應(yīng)的能量值由0.8(Re=100)增加到1.0(Re=160)。剪切來流作用下，中下游兩圓柱體橫流向位移PSD曲線對(duì)應(yīng)的主頻均減小，如圖14(b)所示。另外，大振幅區(qū)間內(nèi)下游圓柱體PSD曲線中出現(xiàn)復(fù)雜的“多頻”現(xiàn)象，會(huì)影響能量在流體與結(jié)構(gòu)之間的傳遞，從而改變結(jié)構(gòu)的振動(dòng)特性。

圖14 當(dāng)Re=160時(shí)，不同剪切率工況下中下游兩圓柱體結(jié)構(gòu)橫流向位移能量譜密度隨折減速度的變化情況Fig.14 The variation of the power spectral density of the crossflow displacement with shear ratio and reduced velocity when Re=160

3.5 尾流模態(tài)

圖15和圖16分別給出了Re=100和Re=160時(shí)，部分工況下三圓柱體結(jié)構(gòu)群尾流場(chǎng)瞬時(shí)渦量圖，所對(duì)應(yīng)的時(shí)刻均為下游圓柱最大正位移處。由圖可知，剪切率的變化對(duì)三圓柱體結(jié)構(gòu)上下兩側(cè)尾流的分布影響較大。當(dāng)Re=100時(shí)，均勻來流作用下，周圍流場(chǎng)強(qiáng)弱分布的差異較小，尾流場(chǎng)的旋渦分布相對(duì)來說較為規(guī)則。而剪切來流作用下，柱體表面自由來流的速度由上至下逐漸減小，使得三圓柱體結(jié)構(gòu)上表面的旋渦明顯強(qiáng)于下表面，下側(cè)旋渦在向尾流區(qū)域發(fā)展的過程中受到上側(cè)渦流的擠壓、結(jié)構(gòu)的運(yùn)動(dòng)以及流體之間慣性力的影響較大，在尾流場(chǎng)的分布會(huì)呈現(xiàn)出較大的不規(guī)則性，如圖15所示。尾流模態(tài)隨折減速度的具體變化過程如下：均勻來流作用下，當(dāng)Ur=3時(shí)，結(jié)構(gòu)振動(dòng)較為穩(wěn)定，尾流模態(tài)呈現(xiàn)為標(biāo)準(zhǔn)的雙行2S模式(S：?jiǎn)螠u)。當(dāng)Ur=5時(shí)，上游圓柱的振幅增大，2S模態(tài)的橫流向間距也隨之增大，中上游兩圓柱體尾流對(duì)下游圓柱的“屏蔽”作用進(jìn)一步增強(qiáng)，周圍流場(chǎng)速度的減小使得結(jié)構(gòu)的橫向流體力減小，從而減弱振動(dòng)幅度Ymax/D=0，如圖5(c)所示。隨折減速度進(jìn)一步增大，脫落旋渦與結(jié)構(gòu)之間的耦合作用增強(qiáng)，中下游兩圓柱逐漸激發(fā)出大幅振動(dòng)。特別是當(dāng)Ur=7時(shí)，中游圓柱的振幅達(dá)到峰值，其上下表面脫落的旋渦受到結(jié)構(gòu)強(qiáng)烈的擺動(dòng)作用而明顯加長(zhǎng)，從而正面撞擊下游圓柱并在其上下兩側(cè)分裂成一大一小兩個(gè)子旋渦，在下游圓柱表面形成一對(duì)相反的作用力，使得結(jié)構(gòu)表面升力系數(shù)接近于零，振動(dòng)幅度減弱，如圖8(c)所示。當(dāng)Ur≥9時(shí)，結(jié)構(gòu)尾流場(chǎng)呈現(xiàn)為規(guī)則的2S模態(tài)，且隨著結(jié)構(gòu)振動(dòng)響應(yīng)逐漸減小而趨于穩(wěn)定，渦街的橫向間距以及旋渦數(shù)量明顯減小。另外，下游圓柱周圍渦街的分布均呈現(xiàn)為中上游圓柱下表面脫落的旋渦恰好或即將到達(dá)下游圓柱下表面，而上側(cè)旋渦即將脫落的形式，但是四種折減速度工況下下游圓柱的振動(dòng)幅度卻有較大的差異，這主要是旋渦在下游圓柱表面所產(chǎn)生的壓力區(qū)的強(qiáng)弱差異造成的。剪切來流作用下，大多數(shù)折減速度工況下，由于下側(cè)旋渦強(qiáng)度較弱且距離較遠(yuǎn)，下游圓柱的運(yùn)動(dòng)主要受到上側(cè)渦流的影響。當(dāng)Ur=3和Ur=5時(shí)，三圓柱體結(jié)構(gòu)尾流模態(tài)呈現(xiàn)為上下兩側(cè)渦量差異較大的2S模態(tài)，由于旋渦的強(qiáng)度以及脫渦頻率較大，上側(cè)渦流在向遠(yuǎn)尾流區(qū)域發(fā)展的過程中隨著能量的耗散會(huì)出現(xiàn)渦量堆積的現(xiàn)象。另外，中上游兩圓柱體對(duì)下游圓柱的“屏蔽”作用會(huì)持續(xù)至Ur=6工況。結(jié)構(gòu)產(chǎn)生大振幅響應(yīng)所對(duì)應(yīng)的折減速度區(qū)域內(nèi)，尾流場(chǎng)旋渦強(qiáng)度較大且分布的不規(guī)則性進(jìn)一步增強(qiáng)。值得注意的是，當(dāng)Ur=11時(shí)，中下游兩圓柱體周圍渦流的數(shù)量明顯增大，在MC尾流區(qū)呈現(xiàn)出P+S模態(tài)(P：對(duì)渦)。同時(shí)，耦合作用的變化使得柱體上下兩側(cè)渦量差減小，引起DC發(fā)生劇烈的振動(dòng)響應(yīng)Ymax/D=1.32，如圖5(c)所示。隨著折減速度的進(jìn)一步增大，到Ur=40時(shí)，中游圓柱與周圍旋渦的耦合作用與均勻來流工況類似，然而下游圓柱卻有很大的差異，尾流模態(tài)類似于Ur=5工況，使得下游圓柱的振幅逐漸趨近于零，遠(yuǎn)小于k=0.00工況。

圖15 當(dāng)Re=100時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)瞬時(shí)渦量圖隨折減速度的變化情況Fig.15 The variation of the instantaneous vorticity contourswith shear ratio and reduced velocity when Re=100

當(dāng)Re=160時(shí)，流場(chǎng)的復(fù)雜化使得尾流分布以及向遠(yuǎn)尾流區(qū)域發(fā)展過程當(dāng)中的不規(guī)則性、渦流的數(shù)量和強(qiáng)度明顯增強(qiáng)，如圖16所示。最顯著的變化就是受到結(jié)構(gòu)運(yùn)動(dòng)與慣性力的影響，流場(chǎng)當(dāng)中大部分旋渦的長(zhǎng)度增大，發(fā)生分裂之后會(huì)導(dǎo)致近尾流區(qū)域分布較多的強(qiáng)度較弱的子旋渦。均勻來流作用下，當(dāng)Ur=3和Ur=5時(shí)，中下游圓柱的尾流模態(tài)均呈現(xiàn)出2S模式，但渦街向下游發(fā)展的過程中有發(fā)生偏斜且轉(zhuǎn)變?yōu)榻惶娣植?、間隔較遠(yuǎn)的雙行渦街的趨勢(shì)。當(dāng)Ur=6時(shí)，中游圓柱的振動(dòng)幅度明顯大于Re=100工況，尾流區(qū)域渦街的橫向間距增大。隨折減速度的進(jìn)一步增大，雷諾數(shù)的增大改變了三圓柱體結(jié)構(gòu)運(yùn)動(dòng)與旋渦脫落之間的關(guān)系，從而改變結(jié)構(gòu)周圍旋渦和壓力區(qū)的分布。當(dāng)Ur=9和Ur=18時(shí)，中游圓柱尾流區(qū)均出現(xiàn)了P+S旋渦模態(tài)，但對(duì)下游圓柱的影響有所不同，其橫向振幅相差近兩倍，即Ymax/D=1.17(Ur=9)與Ymax/D=0.43(Ur=18)。剪切來流作用下，當(dāng)Ur=3時(shí)，下游圓柱體尾流區(qū)域的旋渦在較大的慣性力的影響下難以維持2S模態(tài)，在近尾流區(qū)域呈現(xiàn)為不規(guī)則的分布。當(dāng)Ur=5時(shí)，上下兩側(cè)渦流的橫向間距進(jìn)一步增大，這主要是因?yàn)橄聜?cè)強(qiáng)度較弱的渦流發(fā)生偏斜，與均勻來流工況下僅受到結(jié)構(gòu)振幅的影響不同。值得注意的是，當(dāng)Ur=18時(shí)，中游圓柱的尾流會(huì)從下游圓柱上下兩側(cè)通過，有利于其表面旋渦的發(fā)展，從而加劇結(jié)構(gòu)的振動(dòng)，形成了第二個(gè)峰值，如圖6(c)所示。其余折減速度工況下，尾流區(qū)域旋渦的分布則更加的不規(guī)則。

圖16 當(dāng)Re=160時(shí)，不同剪切率工況下串列三圓柱體結(jié)構(gòu)瞬時(shí)渦量圖隨折減速度的變化情況Fig.16 The variation of the instantaneous vorticity contourswith shear ratio and reduced velocity when Re=160

4 結(jié) 論

基于四步半隱式特征線分裂算子有限元方法，本文對(duì)平面剪切來流作用下串列布置三圓柱體結(jié)構(gòu)的單自由度流致振動(dòng)響應(yīng)進(jìn)行了數(shù)值模擬，深入探討了雷諾數(shù)、剪切率與折減速度的變化對(duì)結(jié)構(gòu)群橫流向動(dòng)力響應(yīng)、流體力系數(shù)、相位特征、頻譜特性與尾流模態(tài)的影響，闡釋了流體與結(jié)構(gòu)群之間的耦合效應(yīng)，得到的主要結(jié)論如下：

(1)雷諾數(shù)的增加對(duì)大剪切率工況下結(jié)構(gòu)的振幅響應(yīng)有明顯的促進(jìn)作用。不同的雷諾數(shù)工況下，三圓柱體結(jié)構(gòu)的振幅曲線隨剪切率呈階段性變化，最大振幅值會(huì)隨剪切率的增加而小幅增大。

(2)雷諾數(shù)和剪切率的變化對(duì)下游圓柱的影響大于中上游兩圓柱體。隨雷諾數(shù)的增加，下游圓柱體在部分剪切率工況下升阻力系數(shù)的最大值及其對(duì)應(yīng)的位置會(huì)發(fā)生明顯的改變。剪切率的變化對(duì)下游圓柱體共振區(qū)間內(nèi)的升阻力系數(shù)以及大折減速度范圍內(nèi)的升力系數(shù)均方根隨剪切率變化的波動(dòng)性較大。

(3)隨著雷諾數(shù)和剪切率的增加，中下游兩圓柱體發(fā)生相位轉(zhuǎn)變所對(duì)應(yīng)的區(qū)域范圍明顯擴(kuò)大。值得注意的是，由于中下游兩圓柱體在整個(gè)振動(dòng)周期內(nèi)相位變化的差異較大，使得大多數(shù)工況下兩圓柱體的相位差不會(huì)呈現(xiàn)出完全同相或反相的狀態(tài)。

(4)雷諾數(shù)和剪切率的增大會(huì)顯著增大橫流向位移能量譜密度曲線中頻率的階數(shù)，從而影響能量在流體與結(jié)構(gòu)之間的傳遞，改變結(jié)構(gòu)群的動(dòng)力響應(yīng)。

(5)三圓柱體結(jié)構(gòu)尾流場(chǎng)的旋渦分布隨折減速度變化規(guī)律性較強(qiáng)，主要會(huì)出現(xiàn)2S與P+S等尾流模態(tài)。但隨雷諾數(shù)和剪切率的增加，尾流場(chǎng)的旋渦分布在較大的慣性力以及來流速度差的影響下變得更加復(fù)雜。