郭 濤，張紋惠,2，王文全，羅竹梅

(1. 昆明理工大學(xué)建筑工程學(xué)院工程力學(xué)系，昆明 650500；2. 中國能源建設(shè)集團(tuán)云南省電力設(shè)計院有限公司，昆明 650000；3. 四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室，成都 610065；4. 昆明理工大學(xué)能源與動力工程系，昆明 650093)

流體-結(jié)構(gòu)互動問題是自然界中十分普遍的物理現(xiàn)象，也是一個經(jīng)典的流體力學(xué)問題。當(dāng)流體流過不具有流線型的鈍體結(jié)構(gòu)時，會在固體壁面形成黏性邊界層，在一定雷諾數(shù)范圍內(nèi)該邊界層容易與壁面發(fā)生分離，在結(jié)構(gòu)尾部形成交替式的周期性脫落漩渦，產(chǎn)生周期性的上下拖曳力，也就是升力、阻力，從而誘發(fā)結(jié)構(gòu)振動。渦激振動(Vortex-induces vibration, VIV)是促使結(jié)構(gòu)產(chǎn)生疲勞破壞的重要原因，其危害和影響十分廣泛。例如：空氣動力學(xué)中，機(jī)翼在流動空氣中的振動問題[1-3]，導(dǎo)彈無側(cè)滑大攻角分行時，背風(fēng)面分離渦的不對稱性導(dǎo)致的附加偏航力問題；水力機(jī)械中，導(dǎo)葉、葉片在水流作用下的擺動和振動問題[4-5]；橋梁工程中，震驚世界的塔科馬大橋顫振風(fēng)毀事件(1940)[6]。近期的虎門大橋渦激振動事件(2020 年5 月5 日)[7]，拉索風(fēng)雨激振問題[8-11]等；高層建筑物中，如深圳市賽格廣場大廈風(fēng)致渦激共振造成的有感振動事件(2021 年5 月18 日)[12]；以及海洋立管系統(tǒng)(如生產(chǎn)立管、輸油、輸氣立管等)和系泊系統(tǒng)[13-14]領(lǐng)域十分常見的渦激振動現(xiàn)象等。因此，結(jié)構(gòu)在流場中的渦激振動抑制以及從渦激振動中獲取能量的研究[15-22]問題成為工程界和學(xué)術(shù)界廣泛關(guān)注的問題。

1968 年Feng[23]以空氣為介質(zhì)的高質(zhì)量比渦激振動實驗，奠定了后續(xù)渦激振動的研究基礎(chǔ)。自此，各國學(xué)者對圓柱繞流渦激振動的質(zhì)量比效應(yīng)開展了廣泛的研究[24-39]。Mittal 等[24]通過有限元數(shù)值計算得出，低質(zhì)量比(m?=4.25)單圓柱會出現(xiàn)“弱鎖定”現(xiàn)象，而較高質(zhì)量比(m?=25)未發(fā)現(xiàn)“弱鎖定”現(xiàn)象。Tu 等[25]發(fā)現(xiàn)低質(zhì)量比單圓柱不僅會出現(xiàn)“弱鎖定”現(xiàn)象，同時還會伴有大幅橫流向振動。而且Tu 等[26]和Prasanth 等[27]在低質(zhì)量比串聯(lián)雙圓柱中也發(fā)現(xiàn)了“弱鎖定”現(xiàn)象。Jiang等[28]、楊驍?shù)萚29]和杜曉慶等[30]的研究也表明：質(zhì)量比對雙圓柱渦激振動有顯著影響。谷家揚等[31]和陳正壽等[32]基于CFD 法研究發(fā)現(xiàn)低質(zhì)量比單圓柱渦激振動時，鎖定區(qū)間和橫流向振幅均隨質(zhì)量比的增大而減小。卞正寧等[33]基于SST 湍流模型研究發(fā)現(xiàn)高雷諾數(shù)下，高質(zhì)量比的圓柱渦激振動會出現(xiàn)高幅分支。李小超等[34]研究發(fā)現(xiàn)質(zhì)量比控制和影響了鎖定頻率及最大能量轉(zhuǎn)換率的無量綱流速范圍。實驗方面，Williamson 等[35-36]進(jìn)行了一系列不同質(zhì)量比m?(0.5～20)下的剛性圓柱體渦激振動和尾渦脫落實驗，為后人提供了豐富的參考數(shù)據(jù)，實驗表明：在鎖定區(qū)間和振動幅值方面，低質(zhì)量比情況表現(xiàn)出來的結(jié)果與高質(zhì)量比情況有著明顯差異，而且質(zhì)量比-阻尼組合m?ζ對振動響應(yīng)也有重要影響。

截至目前，對較高質(zhì)量比m?情況下的渦激振動研究較少，1992 年Anagnostopoulos 等[37]做的比較經(jīng)典的彈性支撐剛性圓柱體渦激振動實驗，質(zhì)量比非常高，達(dá)到m?=148。趙劉群等[38]和鄧躍等[39]采用ALE 和CFD 法對該實驗參數(shù)下的低雷諾數(shù)圓柱體VIV 進(jìn)行了數(shù)值計算，得到了非鎖定區(qū)間、鎖定區(qū)間和“拍頻”現(xiàn)象。除此以外，很多數(shù)值研究工作涉及到的質(zhì)量比m?并不高，從0.5～25 不等。2021 年5 月18 日，深圳市賽格廣場大廈在低速風(fēng)場(風(fēng)速發(fā)生時風(fēng)速約為10 m/s)作用下產(chǎn)生了強烈的有感振動問題。經(jīng)多方論證認(rèn)為是由于樓頂桅桿的第四階非對稱風(fēng)致渦激振動(2.12 Hz 頻率)和大廈及桅桿動力特性改變的耦合，誘發(fā)結(jié)構(gòu)的高階彎扭組合振動[40]。以此問題為出發(fā)點，為進(jìn)一步澄清渦激振動的高質(zhì)量比效應(yīng)，本文以彈性支撐形式下的“質(zhì)量-彈簧-阻尼系統(tǒng)”單圓柱體結(jié)構(gòu)(便于和前人研究作對比)為研究對象，基于銳利界面浸入邊界法，采用C++編程計算了低雷諾數(shù)(Re=80～150)流場中，不同高質(zhì)量比(m?=14.8 ～280 )和阻尼比( ζ=0.0012 ～0.036)作用下的結(jié)構(gòu)振動響應(yīng)、結(jié)構(gòu)受力、“鎖定區(qū)間”、相位突變以及尾渦脫落問題。研究了質(zhì)量比對結(jié)構(gòu)振動響應(yīng)特性、升阻力系數(shù)、渦脫頻率等水力不穩(wěn)定現(xiàn)象的影響規(guī)律，分析了尾渦脫落模態(tài)的演變過程，探討了導(dǎo)致結(jié)構(gòu)振動的流固耦合機(jī)制。

1 計算方法

為避免復(fù)雜貼體網(wǎng)格的更新和畸形對動邊界流場計算效率、精度的影響，以充分掌握結(jié)構(gòu)場的受力特性，本文采用浸入邊界法進(jìn)行流固耦合計算。浸入邊界法(Immersed boundary method，IBM)最初由Peskin[41]提出，并成功應(yīng)用到模擬人類心臟血液流動中的瓣膜運動中。早期的IBM 法，是將固體邊界離散為一組由彈簧連接的單元，適用于柔性體[42]，但變形不易通過流場物理量求得，不具有一般性。Goldstein 等[43]發(fā)展了虛擬邊界法來求解剛性邊界問題，主要使用反饋力法，配合譜方法來漸進(jìn)地施加所需的邊界條件，也稱反饋力法。后來Mohd-Yusof[44]提出了直接構(gòu)造受力區(qū)的速度場來確保滿足邊界條件的直接力法。Fadlum等[45]、Le 等[46]、Wu 等[47]、Uhlmann[48]、王文全等[49]和郭濤等[50]將其應(yīng)用于流固耦合研究，獲得了與貼體網(wǎng)格下計算結(jié)果一致性較好的結(jié)果。以上浸入邊界法為擴(kuò)散界面法(Diffused-interface method)，均通過計算力源項來施加浸入邊界條件。與擴(kuò)散界面法相對的另一類方法為銳利界面法(Sharpinterface method)，該方法將力源項的影響通過界面單元處的速度插值直接施加到計算域，即用插值方程替代浸入界面附近處的流體控制方程，實施較簡單。本文采用一種改進(jìn)的銳利界面浸入邊界法[51]模擬復(fù)雜邊界的流激振動問題。該方法將整個物理區(qū)域劃分成純流體區(qū)域以及包含固體的次流體區(qū)域。在流體次區(qū)域，構(gòu)造“虛擬點—受力點—垂足點”的計算結(jié)構(gòu)。利用雙線性插值法對虛擬點最近的4 個歐拉網(wǎng)格點進(jìn)行插值計算，得到虛擬點的速度，并通過強制滿足固體邊界的無滑移條件計算出受力點的速度，將其作為流動邊界條件，實現(xiàn)復(fù)雜動邊界的流動數(shù)值模擬。

1.1 控制方程

將整個物理區(qū)域(包括流體、固體)看作不可壓縮牛頓流體的粘性流動。在整個計算域內(nèi)，流體采用Euler 描述，固體采用Lagrange 描述。其連續(xù)性方程和動量方程可表示為：

1.2 時間推進(jìn)

采用二階時間分布投影法對時間進(jìn)行離散，來求解控制式(1)和式(2)，具體步驟如下：

首先，在次區(qū)域內(nèi)通過算法判斷和標(biāo)識流體節(jié)點、固體節(jié)點，確定受力點。在任意的流體網(wǎng)格點周圍可以找到相鄰的4 個網(wǎng)格點，若相鄰的這4 個網(wǎng)格點中只要有一個為固體點，那么認(rèn)為該點為考慮力源項影響的流體網(wǎng)格點即受力點，如圖1(a)所示；然后，構(gòu)造如圖1(b)所示的“虛擬點—受力點—垂足點”計算結(jié)構(gòu)：過受力點F 向固體邊界做垂線交固體邊界于垂足點B，反向延長垂線到虛擬點I，使點I～點F的距離等于受力點F到垂足點B的距離。利用線性插值來確定受力點F的速度uf。其中，垂足點B在固體邊界上，當(dāng)圓柱為固定邊界時，其速度恒等于零。本文算例中，固體為具有單自由度彈性支承的圓柱體，發(fā)生渦激振動后，固體邊界存在被動振動的橫流向位移。所以，B點的初始速度為零，之后每一時間步下的速度由式(17)求得。

圖1 固體邊界處理Fig. 1 Solid boundary

虛擬點I的速度利用其周圍4 個Euler 網(wǎng)格點速度進(jìn)行雙線性插值得到：

式中：un為距離虛擬點最近的4 個Euler 網(wǎng)格上的點(三點不可在同一直線上)； φi為常系數(shù)：

式中：

式中，下標(biāo)1～4 分別指距離虛擬點I最近的4 個Euler 網(wǎng)格點。

3)考慮新的壓力項，更新流場，求得下一步的速度u，即：

式中，pn+1為通過解耦壓力Poisson 方程式(11)得到的下一步壓力值。

1.3 空間離散

對于對流項，設(shè) φ為流場內(nèi)的任意變量，對于流體區(qū)域內(nèi)部點，采用如下格式進(jìn)行離散：

在邊界位置，對流項采用一階迎風(fēng)格式進(jìn)行離散。

對于擴(kuò)散項，中間節(jié)點處采用中心差分格式，即：

在邊界上，采用一階迎風(fēng)格式。

在對壓力的離散過程中，采用五點差分離散壓力泊松方程：

對于Neumann 邊界條件，在邊界點處采用增設(shè)虛點的方法[52]。對式(14)等號右邊項，內(nèi)部節(jié)點采用中心差分格式：

邊界節(jié)點采用一階迎風(fēng)格式，如在左邊界上(i=0)，有：

2 渦激振動算例

2.1 計算模型

如圖2 所示，計算域為一矩形區(qū)域，長×寬=15D×10D，流體次區(qū)域為邊長3D的正方形，次區(qū)域左側(cè)距流場入口的距離為3.5D。次區(qū)域中央浸沒的剛性圓柱設(shè)置為“質(zhì)量-彈簧-阻尼結(jié)構(gòu)”的彈性支撐形式，如圖3 所示。圓柱直徑為D，圓心坐標(biāo)為(0, 0)。計算域左側(cè)為均勻來流入口邊界，采用Dirichlet 邊界條件，即u=U，v=0；上下兩側(cè)均為無穿透邊界；右側(cè)為自由出流邊界，采用Neumenn 邊界，即 ?u/?x=0 ， ?v/?x=0。整個流場采用一套間距為 Δx=Δy=0.025的均勻四邊形網(wǎng)格，固體邊界采取與流體網(wǎng)格尺度相等的等間距離散，時間步長為 Δt=0.002 ，時間歷程t=160 s，采用本文算法計算機(jī)時為29 h (32 G 內(nèi)存的普通PC 機(jī))。

圖2 整體計算域及邊界條件Fig. 2 Computational domain and boundary conditions

圖3 渦激振動計算模型Fig. 3 Computational model of vortex-induced vibration

在流體作用下，圓柱發(fā)生被動形式的剛性縱向振動，采用一般的運動方程進(jìn)行描述：

式中：y為圓柱縱向位移；m為結(jié)構(gòu)的質(zhì)量；c為結(jié)構(gòu)的阻尼；k為結(jié)構(gòu)的剛度；Fy為縱向的流體力(升力)。 計算過程中所涉及到的參數(shù)均做了無量綱化處理，見表1。

表1 各參數(shù)的無量綱化處理Table 1 The dimensionless parameter

將表1 中的無量綱化參數(shù)代入式(17)中得到無量綱化的運動方程，即：

2.2 不同來流情況對流場的影響

圖4 和圖5 為固定參數(shù)質(zhì)量比m?=148，阻尼比 ζ=0.0012，約化速度U?=1 時，渦脫頻率fs、振幅比Y隨雷諾數(shù)(Re=80～150)的變化規(guī)律。

圖4 渦脫頻率 fs/fn 隨Re 數(shù)的變化規(guī)律Fig. 4 Variation of vortex shedding frequency ( fs/fn)with Reynolds number

圖5 振幅比Y 隨Re 數(shù)的變化規(guī)律Fig. 5 Variation of amplitude (Y) ratio with Reynolds number

從圖4 可以看出，渦脫頻率fs隨著雷諾數(shù)Re的變化趨勢可以分成三個階段。初期：fs隨著Re數(shù)的增大而增大，fs/fn<1這一區(qū)域是“弱鎖定”區(qū)間。第二階段：當(dāng)Re數(shù)進(jìn)入某一區(qū)域后，fs保持不變開始穩(wěn)定下來，“渦脫頻率fs與被圓柱的固有頻率fn所捕獲，fs與fn基本相同”，這一區(qū)域便是渦激振動中的“鎖定區(qū)間”，振動頻率單一，而振幅迅速增加。后期：當(dāng)Re數(shù)超過這一區(qū)間后，渦脫頻率將“擺脫”圓柱的固有頻率，fs再次隨著Re數(shù)的增大而增大，“鎖定現(xiàn)象”消失。實驗[37]所得鎖定區(qū)間為Re=105～125；ALE 方法[38]約為Re=98～105，偏離較遠(yuǎn)；CFD+動網(wǎng)格方法[39]得到的鎖定區(qū)間為Re=90～125，區(qū)間較廣，上限值與實驗相同；而本文方法所得區(qū)間為Re=100～130，區(qū)間更加接近實驗值。從圖5 中也可看出，在鎖定區(qū)間內(nèi)，渦激振動強度將達(dá)到最大，最大振動幅值均在0.5D左右。由于計算方法的差異，達(dá)到最大振幅所對應(yīng)的Re有所不同，但都在鎖定區(qū)間內(nèi)靠近Re數(shù)較小的一側(cè)。本文方法和實驗所得最大振幅對應(yīng)的Re數(shù)都為110。

從圖6 中也可看出升力隨著雷諾數(shù)Re的增大而增大。初期Re數(shù)較小時，流場對圓柱所產(chǎn)生的升力系數(shù)也比較小，而影響結(jié)構(gòu)渦激振動縱向位移的主要作用力就是升力，因此振幅Y也很小，基本趨近于0。進(jìn)入鎖定區(qū)間后，升力增加幅度有限，而振幅卻急劇增加(最大達(dá)十幾倍)，此時結(jié)構(gòu)振動不在依賴于升力，結(jié)構(gòu)發(fā)生了共振，Re=110時，橫流向振幅達(dá)到最大；當(dāng)Re=130 時，結(jié)構(gòu)開始擺脫鎖定狀態(tài)，振動開始變?nèi)?，而且此時橫流向振幅與升力系數(shù)發(fā)生了明顯的“相位突變”現(xiàn)象；當(dāng)Re=150 時，結(jié)構(gòu)已經(jīng)完全擺脫鎖定狀態(tài)，此時不僅渦激振動變?nèi)?、消失，振幅與升力系數(shù)仍然處于反相位狀態(tài)，而且升力系數(shù)隨時間的推進(jìn)逐漸出現(xiàn)了“拍頻”現(xiàn)象。

圖6 不同雷諾數(shù)下升力系數(shù)CL 及振幅比Y 時程曲線Fig. 6 The time-history curve of lift coefficients (CL) and amplitude ratio (Y) with different Reynolds number

圖7 為鎖定區(qū)間內(nèi)，當(dāng)振幅達(dá)最大時，一個周期內(nèi)尾跡渦的演化過程。從圖7 中可知。圓柱后面正、負(fù)漩渦交替脫落，向下游發(fā)展，呈現(xiàn)出典型的卡門渦街現(xiàn)象，其渦脫落模態(tài)為2S 型，圓柱表面的邊界層分布、分離剪切層的卷起和圓柱后方的分離區(qū)域清晰可見。與固定繞流不同的是，隨著圓柱的上下振動，漩渦中心與圓柱中心不再總是處于一條水平線上。

圖7 一個振動周期T 內(nèi)尾跡渦隨時間的演化(Re=110，范圍： -8≤w≤8，間隔1，圖中實線為正值，虛線為負(fù)值)Fig. 7 The evolution of vorticity Isolines with time in a vibration period (Re=110, Range:-8≤w≤8, Spacing 1, Solid lines are positive isolines,Dashed lines are negative isolines)

2.3 質(zhì)量比 m? 、阻尼比 ζ 和質(zhì)量-阻尼比組合m?ζ對渦激振動的影響

圖8 給出了不同質(zhì)量比結(jié)構(gòu)下，振動振幅隨雷諾數(shù)Re的變化規(guī)律。

由圖8 可知，不同高質(zhì)量比m?下，結(jié)構(gòu)振幅曲線趨勢基本相同。當(dāng)m?參數(shù)發(fā)生成倍變化時，曲線的變化并不明顯。結(jié)構(gòu)的鎖定區(qū)間都位于Re=100～130，而且發(fā)生共振時對應(yīng)的雷諾數(shù)均為110，位于鎖定區(qū)間靠近Re較小的一側(cè)，最大振幅也都在0.5D左右，無太大變化，即m?對振幅的影響并不大?，F(xiàn)如今，新能源已被廣泛利用，從渦激振動中獲取海洋能的多少與結(jié)構(gòu)受流場力的作用而產(chǎn)生的振幅大小息息相關(guān)。因此，若只考慮減小結(jié)構(gòu)質(zhì)量比m?的方法，對提高結(jié)構(gòu)振動幅度的效果有限。

圖8 不同質(zhì)量比 m?對比振幅比Y 的影響Fig. 8 Variation of amplitude ratio (Y) with mass ratio ( m?)and Reynolds number

圖9 為相同質(zhì)量比m?=14.8 、不同阻尼比ζ(考慮了振動系統(tǒng)阻尼、彈簧剛度的綜合系數(shù))對渦激振動幅值的影響。

圖9 不同阻尼比 ζ對振幅比Y 的影響Fig. 9 Variation of amplitude ratio (Y) with damping ratio ( ζ)and Reynolds number

由圖9 可知，在質(zhì)量比相同的情況下，不同阻尼比結(jié)構(gòu)系統(tǒng)的振動響應(yīng)隨雷諾數(shù)的變化趨勢大致相同。 ζ的變化并沒有對渦激振動的鎖定區(qū)間造成影響；雖然鎖定區(qū)間和變化趨勢相同，但當(dāng)ζ發(fā)生變化時，在鎖定區(qū)間內(nèi)，圓柱發(fā)生渦激振動的最大響應(yīng)幅值是有差別的，隨著 ζ的增加，振幅比Y逐漸減小，振動受到了抑制。所以，質(zhì)量比m?相同的情況下，減小阻尼比ζ，無法使圓柱的鎖定區(qū)間增大，但是對提高結(jié)構(gòu)振動幅度的效果明顯。

圖10 討論了高質(zhì)量比m?=14.8和較高質(zhì)量比m?=148結(jié)構(gòu)系統(tǒng)，在質(zhì)量阻尼比相同情況下(m?ζ=0.1776)，渦激振動振幅隨雷諾數(shù)變化的關(guān)系曲線。

圖10 質(zhì)量阻尼比對縱向振幅的影響Fig. 10 Variation of amplitude ratio (Y) with mass damping ratio ( m?ζ) and Reynolds number

從圖10 可以得到，當(dāng)m?ζ保持不變時，圓柱的渦激振動幅值基本相同，最大值還是在0.5D左右；但是，鎖定區(qū)間卻發(fā)生了較大變化。相對低的質(zhì)量比結(jié)構(gòu)系統(tǒng)，由非鎖定狀態(tài)進(jìn)入鎖定區(qū)間的雷諾數(shù)更低，鎖定狀態(tài)持續(xù)時間長，擺脫鎖定狀態(tài)也晚。因此，質(zhì)量-阻尼比相同時，質(zhì)量比低的結(jié)構(gòu)系統(tǒng)發(fā)生渦激振動時的鎖定區(qū)間更廣，雷諾數(shù)大致為Re=90～140，比m?=148 的較高質(zhì)量比結(jié)構(gòu)系統(tǒng)提高了1.67 倍，能更早的進(jìn)入鎖定區(qū)間。而且，質(zhì)量比低時，結(jié)構(gòu)系統(tǒng)發(fā)生最大振動響應(yīng)時對應(yīng)的雷諾數(shù)也小。

3 結(jié)論

本文利用一種銳利界面浸入邊界法對二維彈性支承圓柱流固耦合問題進(jìn)行了數(shù)值計算，并討論了不同參數(shù)(雷諾數(shù)Re、阻尼比 ζ 、質(zhì)量比m?)對圓柱縱向渦激振動的影響，通過對比分析，得到以下結(jié)論：

(1)通過與其它文獻(xiàn)和實驗結(jié)果的對比，驗證了本方法的準(zhǔn)確性和有效性。

(2)對于圓柱繞流問題，當(dāng)圓柱發(fā)生渦激振動時，在鎖定區(qū)間內(nèi)，振動強度將達(dá)到最大。本文計算結(jié)果與實驗結(jié)果[37]均表明最大振動幅值均在0.5D左右，此時所對應(yīng)的雷諾數(shù)都為Re=110，都在鎖定區(qū)間內(nèi)靠近雷諾數(shù)較小的一側(cè)。實驗所得鎖定區(qū)間為Re=105～125，本文方法所得區(qū)間為Re=100～130。

(3)質(zhì)量比m?對鎖定區(qū)間、最大振動幅值以及發(fā)生共振時對應(yīng)的雷諾數(shù)都沒有太大影響；同樣，阻尼比 ζ對鎖定區(qū)間并無太大影響。但是，對鎖定區(qū)間內(nèi)的最大振動幅值是有影響的，隨著 ζ的增加，振動受到了抑制。因此，相反若是對于海洋能發(fā)電等新能源領(lǐng)域，應(yīng)該考慮從減小振動系統(tǒng)阻尼、彈簧剛度等方面著手，以便獲取更多海洋能。

(4)質(zhì)量阻尼比m?ζ相同情況下，雖然振動振幅不會有明顯變化，但是質(zhì)量比低的結(jié)構(gòu)系統(tǒng)鎖定區(qū)間更廣。發(fā)生最大振動響應(yīng)時對應(yīng)的雷諾數(shù)也小。