懸索橋雙吊索尾流致渦激振動的大渦模擬

趙 燕， 林偉群， 杜曉慶,2， 楊 驍， 代 欽

(1. 上海大學 土木工程系，上海 200444； 2. 上海大學 風工程和氣動控制研究中心，上海 200072；3. 上海大學 上海市應用數學和力學研究所，上海 200072)

雙吊索在大跨度懸索橋上有著廣泛的應用，下游索受到上游索尾流干擾會發生尾流激振現象[1]。Fujino等[2]將雙索的尾流激振分為尾流致渦激振動、尾流馳振和尾流致顫振等振動類型。隨著懸索橋的跨度不斷增大，雙吊索發生尾流激振的可能性大幅增大，因而雙索的尾流激振研究受到眾多學者關注[3-6]。

目前，研究者對雙索的尾流馳振或尾流致顫振進行了較為深入的試驗研究。杜曉慶等研究了錯列雙圓柱的尾流馳振，并發現了兩種不同的振動形式。Fujino等研究了雙索尾流馳振的起振條件。Li等研究了風攻角和雙索間距對尾流馳振的影響等。肖春云等研究了雙圓柱處于尾流的圓柱的平均氣動力系數隨空間位置的變化規律及發生尾流弛振的不穩定區間。

相對來說，目前對雙索尾流致渦激振動的研究較少。杜曉慶等在研究了間距為4D(D為圓柱或雙索的直徑)的串列和錯列雙圓柱的尾流致渦激振動時發現，無論是錯列還是串列布置，下游圓柱的尾流致渦激共振的振幅均遠大于單圓柱渦激振動的振幅。Wong[7]研究了單圓柱及不同間距比、不同阻尼系數下各種串列圓柱組合的尾流致振。Brika等[8]設計了一套可使上下游圓柱同相位振動和反相位振動的試驗裝置，對串列和錯列的雙圓柱的渦激振動進行了研究。Kim等[9]對串列雙圓柱的風致振動進行了風洞試驗研究，并根據上下游圓柱的響應按照圓心間距劃為五個區域。

隨著流體力學理論和計算流體動力學的不斷發展，近年來研究者嘗試采用數值模擬方法來研究尾流致振問題。Carmo等[10]對串列雙圓柱渦激振動進行了數值模擬研究。結果顯示：當發生共振時，下游圓柱的振幅約為單圓柱時的2倍；下游圓柱在高折減速度下，有較大的振幅，這與同樣情況下單圓柱的較低振幅不同。Mysa等[11]對雷諾數Re=100、質量比m*=m/(0.25ρπD2)≈2.6(m為圓柱的單位長質量，ρ為流體密度)的串列雙圓柱尾流致渦激振動進行數值模擬，研究表明：流體力與結構響應之間的相位差在流體與結構之間的耦合中起著主要作用。Bao等[12]對串列雙圓柱雙自由度渦激振動進行了數值模擬，研究發現上下游圓柱同時發生共振現象在各個頻率比下都會發生，且上游圓柱的振動情況與同樣情況下的單圓柱的類似。郭曉玲等[13]對雷諾數Re=150的雙圓柱進行了數值模擬，研究了質量比、圓心間距及折減速度對下游圓柱尾流致渦激振動的影響。

與試驗相比，數值模擬方法更利于研究多圓柱風致振動的流固耦合機理。以往的數值模擬研究大多局限在低雷諾數(Re<200)的層流中進行的，而且計算大多采用二維數值模型，模型的質量比也較低(m*<20)，這與雙索的尾流致渦激振動的條件相差甚遠。圓柱型結構的氣動性能有強烈的雷諾數效應，為了揭示雙索尾流致渦激振動發生的內在機理，有必要在高雷諾數(Re>1×104)下進行尾流致渦激振動數值模擬。

本文采用大渦模擬方法，在較高的雷諾數范圍內(Re=1×104～4×104)，以圓心間距為4D的串列雙圓柱為研究對象，進行了尾流致渦激振動數值模擬，并與作者先前的試驗研究結果進行比較，研究了下游圓柱的動力響應和氣動力隨著折減風速的變化規律，探討了下游圓柱動力響應、繞流場特性以及氣動力三者之間的耦合關系，考察了上游圓柱的尾流對下游圓柱的影響方式。

1 數值方法

1.1 湍流模型的選擇

目前湍流數值模擬方法主要分為三種：直接數值模擬(Direct Numerical Simulation，DNS)、雷諾平均法(Reynolds Average Navier-Stokes，RANS)和大渦模擬(Large Eddy Simulation，LES)。DNS方法可以獲得湍流場的精確信息，但現有的計算資源往往難以滿足對高雷諾數流動模擬的需要。RANS方法可以計算高雷諾數的復雜流動，但只能給出平均運動結果。LES方法采用瞬時的N-S方程直接模擬湍流中的大尺度旋渦，而小尺度旋渦可以通過近似模型來考慮。與DNS法相比，LES法可以節約計算機資源。與RANS法相比，LES法可以更準確地模擬圓柱的邊界層分離、渦的形成及其發展過程，因而本文采用LES方法。

繞流問題中黏性流體不可壓連續方程和瞬時N-S方程濾波整理后，為[14-16]

(1)

(2)

(3)

構造亞格子應力τij的Smagorinsky-Lilly封閉模式為

(4)

亞格子Smagorinsky-Lilly模型中的亞格子湍流黏性系數為

(5)

Ls=min(κd,CsV1/3)

(6)

式中：κ為Von Karman常數；d為到最近壁面的距離；Cs為Smagorinsky常數； 本文計算中Cs取0.1；V為計算單元的面積。

1.2 結構控制方程

考慮圖1所示上游圓柱固定、下游圓柱可作順風向和橫風向的運動的雙圓柱系統，并將下游圓柱等效為質量-彈簧-阻尼系統，其兩自由度振動方程為

(7)

(8)

圖1 雙圓柱計算模型示意圖Fig.1 Computational model of two circular cylinders

本文的數值模擬圓柱與流場之間的流固耦合作用，是通過動網格技術來實現的，其步驟如下：

步驟1運用大渦模擬方法求解流體控制方程，獲得流場速度場、壓力場和下游圓柱表面的流體力；

步驟2將流體力作用于下游圓柱，以四階Runge Kutta法求解下游圓柱運動控制方程式(7)和式(8)，得到下游圓柱橫風向和順風向的動力響應；

步驟3通過動網格技術，將下游圓柱的速度傳遞于網格系統，更新網格位置；

步驟4返回步驟1開始計算下一個時間步的響應，如此循環可得各時間步的下游圓柱動力響應，實現上述流固耦合算法。

1.3 網格劃分及計算參數

計算域與網格劃分如圖2和圖3所示。將計算域分塊進行網格劃分，在雙圓柱附近采用圓形網格劃分，并在靠近圓柱表面處進行加密。所有工況的圓柱近壁面的0.35D的厚度內采用結構化網格，如圖3所示。下游圓柱在圖2所示的設置為界面(interface)邊界的圓形流域內運動，其他區域為靜止區域。

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

圖3 計算網格Fig.3 Computational grid

針對圓心間距為4D的串列圓柱尾流致渦激振動，采用杜曉慶等的風洞試驗相關的參數：圓柱直徑D=0.18 m，雷諾數為1×104～4×104，折減風速Ur=3～10，在靜止空氣中的自振頻率fn=1.73 Hz。由于大渦模擬的計算量較大，為了減小起振時間，降低計算量，阻尼比取0，質量比m*取40。

壓力和速度耦合采用SIMPLEC法求解，動量方程和湍動能耗散率方程采用二階精度的離散格式。無量綱時間步長Δt*為0.005(Δt*=ΔtU/D， 其中Δt為非定常計算時間步長)。

邊界條件設定：入流面采用速度入口(velocity-inlet)邊界條件；出流面采用自由出流(outflow)邊界條件；頂面和底面采用對稱(symmetry)邊界條件；兩側面(圓柱展向兩端的平面)采用周期性(periodic)邊界條件；串列圓柱表面采用無滑移壁面(wall)邊界條件；圓柱體的運動區域與靜止區域之間的交界面采用界面(interface)邊界條件。

2 計算結果及分析

2.1 計算結果驗證

為了保證本文所采用的計算方法和網格模型的可靠性，以固定單圓柱為研究對象進行結果驗證。在雷諾數為4×104時，分析了周向網格數量和展向網格數量對計算結果的影響，并將計算結果與文獻內的風洞試驗和數值模擬結果進行了比較。

表1列出了本文計算得到的平均阻力系數、脈動升力系數和Strouhal數(St)等結果與文獻結果的比較，而圖4和圖5給出了本文計算得到的圓柱表面平均風壓系數和脈動風壓系數與文獻值的比較。

從表1、圖4和圖5可知：隨著周向網格從120增加至256，本文的三個工況A1，A2和A3的風壓系數、氣動力系數和St數較為接近，體現了本文計算結果的網格獨立性；而這三個工況的數值與文獻結果相比發現，St數較為吻合，但平均阻力、脈動升力和表面風壓系數有一定的偏差。為了驗證展向網格對結果的影響，A4工況在A2工況的基礎上，將展向的網格數量從10個增加至15個，由表1、圖4和圖5可知，A4工況的氣動力系數和表面風壓系數更為接近于文獻值。這說明增加展向的網格數量可以提高計算的精度。但限于計算資源，下文仍參考A2工況的計算參數和網格來建立尾流致渦激振動的計算模型。

表1 靜止圓柱的確定與驗證

圖6給出了串列雙圓柱的尾流模態圖，在振幅較小、折減風速較低的Ur=3.0情況下，下游圓柱的尾流呈現出“2S”的模態(“2S”的模態，即：每半個振動周期從下游圓柱脫落單個旋渦)；在振幅較大的Ur=5.8時，下游圓柱的尾流呈現出“2P”的模態(“2P”的模態，即：每半個振動周期從下游圓柱脫落一對旋渦)。這與文獻[22]的試驗結果相吻合，再次驗證了計算結果的正確性。

圖4 平均風壓系數Fig.4 Mean pressure coefficient

圖5 脈動風壓系數Fig.5 RMS pressure coefficient

圖6 Ur=3和Ur=5.8的串列雙圓柱的尾流模態Fig.6 Wake mode of tandem cylinders under Ur=3 and Ur=5.8

2.2 動力響應特性

圖7為下游圓柱的振幅隨折減風速的變化曲線，圖7也給出了杜曉慶等研究中的風洞試驗結果。本文中所提到的振幅指的是無量綱化的振動幅度的最大值。其中，xmax/D,ymax/D分別為順風向和橫風向的無量綱化的最大振幅。

由圖7可知，下游圓柱的振動主要發生在橫風向，下游圓柱的順風向的振幅較小，振幅隨折減風速的變化趨勢與風洞試驗結果相似，這說明數值模擬結果與試驗結果吻合良好。

圖7 振幅隨折減風速的變化情況Fig.7 Variation of the amplitude with the reduced velocity

本文的計算工況在折減風速Ur<5時，下游圓柱的振動幅值較小，還未進入渦激共振區；隨著折減風速的進一步著增大，振幅開始有明顯增大，在Ur=5.0時，振幅達到0.35D； 當折減風速Ur=5.8時，下游圓柱的橫風向的振幅達到最大值，約為0.47D；隨著折減風速的進一步增大，圓柱的橫風向的振幅開始減小。其中，共振發生在Ur=5.0～7.0。

而文獻[1]的最大振幅發生在Ur=6.0時，最大振幅達到0.35D，共振發生在Ur=5.5～7.0。相比而言，本文的計算工況在橫風向發生共振的折減風速范圍和最大振幅均較大，最大振幅發生的折減風速較小。這可能是由于本文所采用的質量比和阻尼較文獻[1]小的緣故。

圖8為Ur=5時下游圓柱的橫風向的位移y/D的時程曲線。其中，y/D為無量綱化的橫風向位移。圓柱的橫向振動過程可以觀察到“拍”的現象，而對于Ur<5及Ur≥7的工況，則沒有“拍”的現象。“拍”是渦激振動的一個重要的特征，是由于圓柱的固有頻率與渦激振動頻率相近而產生的共振不穩定導致的。

圖8 “拍”現象Fig.8 The “beat” phenomenon

2.3 氣動力特性

圖9給出了下游圓柱的橫風向的振幅ymax/D、平均阻力系數CD2,mean和脈動升力系數CL2,rms隨折減風速的變化曲線。圖中也給出了本文對雷諾數為40 000的固定的圓心間距為4D的串列雙圓柱的大渦模擬的計算結果，即下游圓柱的平均阻力系數CD2，mean-sta和脈動升力系數CL2，rms-sta。

從圖9可知，橫風向振幅較小時(在Ur<5和Ur≥7時)，帶兩自由度的下游圓柱的平均阻力系數基本與固定圓柱的結果較為接近。而在發生大幅渦激振動時(5

當Ur=3和Ur=10時，隨著時間的推移，下游圓柱的升力系數的幅值變化與固定雙圓柱的基本類似，且升力系數的幅值大于1；而當Ur=5.8時，隨著時間的推移，升力系數的幅值時而增大至與圖10其他工況變化類似，時而減小至幅值小于1(圖10(b)中t=13.5～14.5 s附近)，其波動與圖10中的其他的工況相比，較為明顯。Ur=5.8工況的這種隨時間變化，升力系數時而增大時而減小，可能是造成Ur=5.8的脈動升力系數并沒有顯著增大的原因。另外，從圖10(b)可知，在Ur=5.8時，隨著下游圓柱瞬時的橫風向的位移的增大，脈動升力系數反而減小，這是由于氣動力與位移之間的相位差造成的。

圖9 氣動力和最大振幅隨折減風速的變化情況Fig.9 Variation of the maximum amplitude and aerodynamic forces with the reduced velocity

為了進一步解釋脈動升力系數的特點，圖10給出了固定的串列雙圓柱的下游圓柱的升力系數的時程曲線及三個典型計算工況(Ur=3，5.8，10)的下游圓柱的橫風向位移及升力系數的時程曲線。

2.4 氣動力與振動響應關系

圖11給出了在Ur=5.8時，下游圓柱的氣動力系數、振幅、瞬時的能量輸入(其中，正值代表氣動力對結構做功，負值代表氣動力消耗結構的能量，FL(t)代表瞬時的氣動升力，y(t)代表瞬時的橫風向位移)、采用希爾伯特變換[23]得到的橫風向的位移與氣動升力的相位差的時程曲線(其中，正值代表橫風向位移的瞬時相位滯后于氣動升力，負值代表橫風向位移的瞬時相位領先于氣動升力)。

圖10 氣動力和位移時程曲線Fig.10 Time history of the displacement and the lift

由圖11的瞬態能量輸入圖可以看出：對于高雷諾數的串列雙圓柱的渦激振動，下游圓柱的氣動力對下游圓柱輸入能量或者消耗下游圓柱的能量，這會引起下游圓柱的氣動力與振動之間產生不穩定的相位差，同時，相位差對下游圓柱的振幅的變化也有重要影響。

當尾流致渦激振動發生時，在下游圓柱的橫風向的振幅逐漸增大過程中，位移的瞬時相位領先于升力，在一個振動周期內升力對下游圓柱做正功，而位移與升力之間的相位差則逐漸增大；反之，下游圓柱的橫風向的振幅逐漸減小時，位移的相位滯后于氣動升力，在一個振動周期內氣動升力對結構系統做負功，下游圓柱的氣動阻力系數出現負值，氣動升力與位移間的相位差也隨之減小。

圖11 Ur=5.8時的位移、氣動力系數的時程曲線及能量輸入模式及相位差Fig.11 Time history of displacement and aerodynamic force coefficients, energy input mode and the phase difference between the displacement and the lift under Ur=5.8

2.5 尾流干擾模態

在工況Ur=5.8的兩個連續的振動周期內選取如圖12所示九個不同典型時刻可知，兩個振動周期包含了振幅達到最大及開始減小的過程。圖13為圖12中T1～T9的典型時刻的瞬態渦量圖。其中，圖13中的“1”、“2”、“3”、“4”代表上游圓柱依次脫落的旋渦的編號。

圖12 Ur=5.8時的位移、氣動力系數的時程曲線Fig.12 Time history of displacement and aerodynamic force coefficients under Ur=5.8

從圖13(a)～圖13(f)可知：“1”號旋渦在T2時刻自上游圓柱脫落后，在T4時刻，由于下游圓柱偏離平衡位置相對較遠，旋渦從下游圓柱的外側流經，并與下游圓柱的剪切層發生相互作用。在T1～T5時刻，氣動升力與位移之間基本是屬于同相位，同時振幅在T2時刻達到最大。

從圖13(a)～圖13(i)可知：“2”、“3”和“4”號旋渦自上游圓柱脫落后，分別在T6、T7和T9時刻時，直接撞擊到下游圓柱的表面，這是由于下游圓柱在平衡位置附近。不同時刻撞擊的位置也不同，并與下游圓柱不同位置的剪切層(取決于上游圓柱的旋渦撞擊到下游圓柱表面的位置)發生作用。從T5～T9時刻的推移，氣動升力與位移之間的相位差有明顯增加，振幅逐漸減小。

綜上，在來流的作用下，上游圓柱的上下兩側會形成交替脫落的旋渦，它們對下游圓柱的影響有上述兩種不同的方式。其中，上游圓柱的旋渦撞擊到下游圓柱的表面的這種影響方式對下游圓柱的氣動力和位移之間的相位差影響較大。

圖13 Ur=5.8時不同時刻的渦量Fig.13 Vorticity at different moments under Ur=5.8

3 結 論

本文采用CFD(Coputation Fluid Dynamics)大渦模擬方法及動網格技術，對串列雙圓柱尾流致渦激振動的流固耦合機理進行了研究。主要結論如下：

(1) 大渦模擬得到的振動響應與文獻中的趨勢類似，在尾流致渦激振動的鎖定區域，橫風向的振動響應能很好地捕捉到“拍”的現象，尾流的模態也與文獻結果一致。表明數值模擬結果與試驗結果吻合良好，在高雷諾數下采用大渦模擬分析雙圓柱尾流致渦激振動問題可獲得可靠的結果。

(2) 在不同的折減風速下，帶兩自由度的下游圓柱的平均阻力系數與其振幅的變化趨勢相同，即下游圓柱的平均阻力系數隨著最大振幅的增加而增大，但下游圓柱的脈動升力系數沒有這樣的變化趨勢，這是由于氣動力與位移之間的相位差造成的。

(3) 當尾流致渦激振動發生時，在下游圓柱的振幅增大的過程中，位移的瞬時相位領先于氣動升力，在一個振動周期內氣動力對下游圓柱做正功，下游圓柱的位移與氣動升力的相位差逐漸增大；反之，在下游圓柱的振幅逐漸減小的過程中，位移的相位相對于氣動升力的相位滯后時，氣動力對下游圓柱做負功，下游圓柱的位移與氣動升力的相位差逐漸減小。

(4) 上游圓柱的尾流對下游圓柱的影響有兩種不同的方式，一種是當下游圓柱偏離平衡位置較遠時，從上游圓柱脫落的上游圓柱脫落的旋渦從下游圓柱的外側流經，與下游圓柱的剪切層發生相互作用；另外一種是當下游圓柱在平衡位置附近時，上游圓柱的漩渦則直接撞擊到下游圓柱的迎風面，不同時刻撞擊的位置也不同。

需要指出的是，本文僅考慮了兩個圓柱和來流處在同一直線上的串列布置雙圓柱的尾流致渦激振動，實際吊索會受到來自不同風向的來流作用，還需考慮風向角對尾流激振的影響。