水流作用下懸浮隧道耦合動力數(shù)值模擬分析

郭曉玲，金瑞佳，劉名名，耿寶磊

(1.國家知識產(chǎn)權局專利局專利審查協(xié)作天津中心，天津 300304；2.交通運輸部天津水運工程科學研究所 港口水工建筑技術國家工程實驗室 工程泥沙交通行業(yè)重點實驗室，天津 300456；3.中國海洋大學，青島 266100；4.成都理工大學，成都 610059)

水下懸浮隧道(Submerged Floating Tunnel，簡稱SFT)作為一種跨越海峽、海灣以及湖泊等的新型交通結構物，通常懸浮在水下方30 m左右。水下懸浮隧道內(nèi)部空間較大，足以滿足公路甚至鐵路等交通方式的要求。有些自然條件十分惡劣的峽灣，由于環(huán)境條件及技術的限制，傳統(tǒng)的跨越方式(如：跨海大橋、沉管隧道)并不可行，而SFT提供了跨越的可能性。耿寶磊等[1]對懸浮隧道水動力問題進行過概述，指出對于工作水深較深的懸浮隧道而言，波浪對其影響較小，而水流對其影響較大。

1 研究背景與現(xiàn)狀

針對水流與懸浮隧道的耦合動力分析，很多學者都有過相關研究。麥繼婷[2]考慮流固耦合效應，研究了表面流速、隧道斷面形式和支撐形式對懸浮隧道動力響應的影響。分析發(fā)現(xiàn)，表面流速會顯著影響懸浮隧道系統(tǒng)的響應位移，但基本不影響沿隧道軸向的應力分布形態(tài)。王廣地[3]分析了海流橫向升力作用下懸浮隧道結構所受荷載隨懸浮深度、海水深度、海流速度以及管段截面尺寸等參數(shù)的變化規(guī)律。龍旭[4]研究了不同的浮重比下水中懸浮隧道的動力響應，給出選用模型的水中懸浮隧道在水流載荷下的最優(yōu)浮重比范圍。上述研究在討論水流荷載時，均采用基于Morison公式的升力經(jīng)驗公式，并未詳細考慮懸浮隧道對流場的影響。為了準確描述水流經(jīng)過懸浮隧道后流場的變化，羅剛等[5-6]通過基于大渦模擬(Large Eddy Simulation, LES)的方法在數(shù)值上分析比較了斷面形式不同、來流速度有差異以及迎流面寬度也不盡相同的情況下固定式懸浮隧道所受的流體作用力及周圍流場分布情況，研究發(fā)現(xiàn)耳形斷面形式懸浮隧道結構物具有最好的穩(wěn)定性，并且其升力系數(shù)和阻力系數(shù)較小，是最為合理的截面形狀，然后依次為圓形、橢圓、六邊形和矩形。

由于懸浮隧道通過錨泊系統(tǒng)懸浮在水中，在流體作用力下就會發(fā)生運動，周期變化的升力使彈性支撐的懸浮隧道在垂直于水流方向發(fā)生振動，也就是“渦激振動”(Vortex-Induced Vibration，VIV)。當渦脫落的頻率接近懸浮隧道結構的自振頻率時，共振或者說鎖定(Lock-in)現(xiàn)象將會發(fā)生。鎖定的發(fā)生同時意味著結構發(fā)生了最不愿意出現(xiàn)的共振現(xiàn)象，此時運動振幅將明顯增大，更會顯著增大結構的疲勞損傷程度，從而對工程安全造成不利影響。國內(nèi)外很多學者對渦激振動問題開展了實驗研究和數(shù)值模擬。Feng[7]早在1968年就于風洞實驗室開展了彈性支撐剛性圓柱的渦激振動物理模型實驗，成為了研究這一現(xiàn)象的經(jīng)典。Williamson和Khalak[8-9]以及Govardhan等[10]于1996年開始，在水槽中進行了低質(zhì)量比彈性支撐剛性柱體的渦激振動物理模型實驗，同樣十分經(jīng)典，成為后續(xù)很多數(shù)值模擬的驗證試驗。隨著數(shù)值模擬技術的飛速發(fā)展，學者們逐漸采用數(shù)值模擬方法對渦激振動問題展開研究，Lu 和Dalton[11]對雷諾數(shù)為13 000 情況下，僅在橫流向發(fā)生運動的圓柱渦激振動問題進行了數(shù)值模擬，模型中的湍流方程采用大渦模擬的方法進行封閉。Dong 和Karniadakis[12]針對雷諾數(shù)為10 000的受迫振動圓柱受力問題，采用了直接數(shù)值模擬(Direct Numerical Simulation，DNS)的方法進行分析研究。郭曉玲等[13]研究了低雷諾數(shù)下串聯(lián)雙圓柱渦激振動機理，主要研究了雙圓柱間干涉作用對渦激振動的影響。何長江和段忠東[14]，方平治和顧明[15]同樣對渦激振動現(xiàn)象進行了數(shù)值模擬。Ge等[16]、Kang等[17]應用von der Pol方程模擬了懸浮隧道在水流作用的渦激振動問題，并研究了懸浮隧道錨鏈間距對隧道振動幅度以及錨鏈受力的影響。Su和Sun[18]應用尾流振蕩模型同樣模擬了懸浮隧道的渦激振動問題，發(fā)現(xiàn)渦激振動共振發(fā)生時，結構的運動軸向應力顯著增加，但這些研究均是采用經(jīng)驗公式研究的懸浮隧道渦激振動問題，沒有采用精確地計算流體力學(Computational Fluid Dynamic, CFD)模型對其模擬。還有其他國內(nèi)外學者都進行過相關研究，但是其研究大多針對某一特定雷諾數(shù)的渦激振動問題，而并沒有計算根據(jù)某一特定工程背景研究一系列雷諾數(shù)下的水流與結構物的耦合運動問題。

基于上述研究背景，本文通過應用有限元數(shù)值方法求解不可壓縮粘性流體雷諾平均Navier-Stokes方程，結合ALE動網(wǎng)格方法，建立了高雷諾數(shù)下水流與結構物耦合的數(shù)值模型，通過精確地CFD數(shù)值模擬，研究了雷諾數(shù)從1 000到100 000這一范圍內(nèi)的水流與懸浮隧道的耦合運動問題，并選取某一特定雷諾數(shù)，研究不同浮重比對懸浮隧道在水流作用下振動的影響，為實際工程提供參考價值。

2 數(shù)學模型的建立

對于水下懸浮隧道而言，由于結構自身尺度較大，同時隧道所處位置通常位于峽灣，水流流速較快，使懸浮隧道周圍的流動通常表現(xiàn)為湍流態(tài)。目前，對于湍流的模擬可以通過DNS或采用合適的湍流模型進行近似，本文采用湍流模型的方法進行模擬。與此同時，由于懸浮隧道結構長度遠大于斷面尺寸，因此我們可以近似為二維水流與結構物的流固耦合問題。二維模型通常會計算出較大的計算結果，但是仍然可以揭示自振頻率與振動幅值和受力特性的關系，很多學者(Lu和Dalton[11]，Dong和Karniadakis[12])同樣也采用二維模型來求解類似問題。

2.1 控制方程及湍流模型

二維不可壓縮粘性流體的雷諾平均N-S方程作為流動基本控制方程在ALE觀點下可以表示為[19-20]

(1)

(2)

(3)

式中：υt為湍流粘性系數(shù)；k為湍動能；δij為Kronecker函數(shù)。

為了封閉湍流方程，本文采用了SST(Shear-Stress Transport)k-ω湍流模型[21-22]。所有參數(shù)均與Menter[21]以及Menter 等[22]相同，其所采用的參數(shù)已被廣泛的接受和成功的應用，當模擬具有負壓力梯度的邊界層時，計算結果表現(xiàn)較好。SSTk-ω湍流模型的控制方程形式如下

(4)

(5)

其中

式中：Ω為渦量的絕對值；y*為到最近壁面的距離；F2和Dkω定義如下

其他參數(shù)通過F1求得

σk=F1σk1+(1-F1)σk2；σω=F1σω1+(1-F1)σω2；α=F1α1+(1-F1)α2；β=F1β1+(1-F1)β2

湍流模型用到的常數(shù)見表1所示。

表1 SST k-ω湍流模型中的參數(shù)Tab.1 Parameters of SST k-ω turbulent model

當計算得到流場和壓力場后，通過對表面壓力和粘性剪切力積分得到結構物所受到的流體作用力，進而得到結構物受到的拖曳力和升力，再對其無量綱化得到拖曳力系數(shù)CD和升力系數(shù)CL

(6)

(7)

2.2 懸浮隧道運動方程

對于懸浮隧道在水流作用下的振動問題，由于需要保證錨索結構始終在彈性范圍內(nèi)，因此盡管懸浮隧道管體會發(fā)生運動，將整個系統(tǒng)仍簡化成質(zhì)量-阻尼-彈簧系統(tǒng)。為了研究更有針對性，本文只考慮懸浮隧道在垂直水流方向(即橫流向)的振動響應，其運動方程可表示為

(8)

式中：m，c，k分別為懸浮隧道的質(zhì)量、阻尼以及剛度；Fy為懸浮隧道受到的垂直波浪力，由上一節(jié)的流動方程確定。

根據(jù)結構動力學關系，c/m=4πξfn，k/m=(2πfn)2，并定義懸浮隧道的質(zhì)量與懸浮隧道排開水的質(zhì)量之比m*=4m/πρD2，D為懸浮隧道外徑，從而可得

(9)

式中：ζ為結構的阻尼比；fn為圓柱的固有頻率；m*為質(zhì)量比。

進一步定義如下的無因次關系式

(10)

式中：U為水流流速，同時根據(jù)升力系數(shù)的定義，F(xiàn)y=0.5ρU2DCL，將式(9)進行整理，可以得到如下無因次結構運動方程

(11)

上式右端的升力系數(shù)前面已經(jīng)給出表達式，即可通過上式計算圓截面懸浮隧道的動力響應。

2.3 計算模型及邊界條件

應用圓柱來模擬懸浮隧道在水流作用下的運動，計算模型及邊界條件如圖1所示，將坐標原點設立于圓柱初始圓心，對計算模型進行無量綱化處理，圓柱直徑D=1。在入口處指定無量綱的速度u=1，v=0；側壁采用了對稱邊界條件?u/?y=0，v=0；出口處的速度邊界條件為?ui/?t+c?ui/?xi=0，其中c為局部平均流速；圓柱表面的邊界條件為不可滑移條件u=dx/dt，v=dy/dt。計算中，在出口處指定相對壓力p=0，在其他邊界采用?p/?n=0的壓力邊界條件，n為指出流體域的外法向單位矢量。在初始時刻，流場中速度及相對壓力分布均設為零(即初始速度場滿足連續(xù)方程)。

圖1 計算域及邊界示意圖Fig.1 Sketch definition of computational domain and boundary conditions

3 數(shù)值離散及網(wǎng)格更新

本文對流擴散方程采用流線迎風有限元(SUPG-FEM)方法[23]進行求解。首先忽略壓力項，只考慮對流和擴散項的中間速度；隨后通過求解壓力方程得到下一個計算時間步的壓力；最后考慮壓力梯度項修正流場。

針對結構運動方程的Newmark-β求解方法，對于時間步進，本文采用動態(tài)時間步長進行求解，表達式如下

(12)

式中：Sc為網(wǎng)格的面積；ue為網(wǎng)格中心點的流速；min表示在計算域內(nèi)取最小值；Cs為安全系數(shù)，取Cs=0.2。而由于懸浮隧道在水流作用下會往復運動，針對該流固耦合問題采用ALE觀點下的動網(wǎng)格方法對其進行數(shù)值模擬。本文將計算域內(nèi)網(wǎng)格給予一定彈性模量，當網(wǎng)格運動時，其運動和變形可以通過求解線彈性動力學控制方程來獲得[24]。應用該網(wǎng)格更新方法，使網(wǎng)格節(jié)點的位移較為均勻，大大提高了數(shù)值計算的穩(wěn)定性；同時，還可以通過調(diào)整計算網(wǎng)格彈性模量來避免網(wǎng)格發(fā)生畸變。

4 數(shù)值模型驗證

表2 網(wǎng)格收斂性驗證結果Tab.2 Grid convergence verification results with different mesh solutions

應用計算收斂的網(wǎng)格，計算不同固有頻率下的彈性支撐圓柱在水流作用下的振動，并同Khalak和Williamson[8]做過經(jīng)典試驗結果進行對比，計算參數(shù)與Khalak和Williamson試驗中完全一致，雷諾數(shù)12 000，質(zhì)量比m*=2.4，質(zhì)量阻尼比m*ξ=0.013，計算結果如圖2所示，圖中橫坐標表示不同的自振頻率，縱坐標為垂直水流方向的振動幅值。從對比結果看，本文的模型所得到的幅值以及“鎖定”區(qū)間和Khalak和Williamson的實驗結果吻合良好。說明本文所建立的模型可以用于高雷諾數(shù)流固耦合問題的研究。

圖2 彈性支撐圓柱計算結果同物理模型試驗結果對比圖Fig.2 Calculation result comparisons of cylinder with spring and damping between numerical model and experimental test

5 算例分析

5.1 不同雷諾數(shù)下懸浮隧道動力分析

基于上述數(shù)值模型，本文計算了懸浮隧道在不同約束剛度，以及不同雷諾數(shù)下的運動情況。雷諾數(shù)從1 000計算到100 000，根據(jù)之前學者的描述，在這一區(qū)間內(nèi)流體流動呈現(xiàn)湍流性質(zhì)，尾流發(fā)生分離。在實際自然環(huán)境中，由于懸浮隧道尺度較大，在自然環(huán)境的水流作用下符合上述雷諾數(shù)。首先介紹雷諾數(shù)50 000，懸浮隧道質(zhì)量比m*=2.5，系統(tǒng)的阻尼ξ=0.007、自振頻率分別為 0.500 Hz、 0.200 Hz和0.083 Hz情況下的垂直水流方向的運動時間歷程曲線，如圖3所示。

3-c 自振頻率fn=0.083 Hz(Re=50 000，m*=2.5，ζ=0.007)圖3 不同自振頻率情況下懸浮隧道運動時間歷程曲線Fig.3 Time histories of SFT motion under different natural vibration frequencies

通過圖3，可以發(fā)現(xiàn)當結構自振頻率為0.5 Hz和0.083 Hz時，懸浮隧道垂直水流方向振動幅值較小，而當結構自振頻率為0.2 Hz時，結構振動幅值較大。隨后取振動幅度較大和較小的工況，分別對其流場進行分析，如圖4和圖5所示。

通過上述對比可以看出，當懸浮隧道系統(tǒng)自振頻率為0.500 Hz和0.083 Hz時，此時渦脫落頻率遠離結構自振頻率，懸浮隧道的渦激振動現(xiàn)象不明顯，僅發(fā)生小幅值振動，其尾渦模式如圖5所示，在同一振動周期內(nèi)尾渦形態(tài)也很規(guī)律；但是當通過調(diào)整彈簧剛度使結構自振頻率為0.200 Hz時，渦脫落頻率與該頻率較為接近，懸浮隧道在水流升力作用下發(fā)生了大幅度的VIV現(xiàn)象，懸浮隧道振動幅值達到了0.7 倍的懸浮隧道外徑，其尾渦模式如圖4所示，可以看出尾渦在剝離開后，有較長的“尾巴”，尾渦形態(tài)也沒有規(guī)律。接下來我們比較該雷諾數(shù)下結構不同自振頻率情況下水流與懸浮隧道的耦合作用，統(tǒng)計得到懸浮隧道在不同頻率下的振動幅值，如圖6所示。

從計算結果可知，當懸浮隧道結構頻率從0.1～0.25 Hz時，在水流作用下結構發(fā)生了“鎖定”現(xiàn)象，而在其他頻率時，結構運動幅值較小，隨后比較相同質(zhì)量比和阻尼情況下，不同雷諾數(shù)下的懸浮隧道渦激振動現(xiàn)象，計算結果如圖7所示。從不同雷諾數(shù)下懸浮隧道的渦激振動結果來看，雷諾數(shù)對振動幅值影響較大。而對發(fā)生“鎖定”的區(qū)間影響不大，無論雷諾數(shù)多大時，發(fā)生“鎖定”的結構自振頻率為0.12～0.25 Hz，因此懸浮隧道的自振頻率要盡量避免在此頻率范圍內(nèi)。

5.2 不同浮重比情況下懸浮隧道動力分析

浮重比對懸浮隧道影響很大，假設懸浮隧道水中尺寸不變，則其受到浮力不變，懸浮隧道浮重比越大，表明自身質(zhì)量越輕，所需建造成本也就越小，但是帶來的弊端就是錨索系統(tǒng)的預張力也就越大，因此為懸浮隧道找到一個合適的浮重比十分重要，接下來分析懸浮隧道浮重比(BWR)1.1、1.3、1.5、1.7和1.9五種情況，對應的質(zhì)量比(m*)分別為2.1、2.3、2.5、2.7和2.9，現(xiàn)計算質(zhì)量比為上述各值時其在水流下的動力響應和受力特征。計算中，隨著懸浮隧道浮重比的變化，其受到的約束條件會相應改變，從而保證結構自振頻率研究區(qū)間不變。再研究雷諾數(shù)Re=100 000，阻尼ζ=0.007的情況，計算結果如圖8和圖9。

圖8 不同浮重比情況下懸浮隧道在不同自振頻率下的振動幅值(Re=100 000，ζ=0.007)Fig.8 Vibration amplitude of SFT versus natural frequency at different BWRs (Re=100 000, ζ=0.007)

圖8為不同浮重比情況下懸浮隧道在不同自振頻率下的振動幅值分析，從結果可以看出，浮重比對振動幅值影響不大，但是對引起“鎖定”的頻率區(qū)間有較大影響，當浮重比較小時，其發(fā)生“鎖定”的區(qū)間較大，而當浮重比大于1.3后，發(fā)生“鎖定”的結構自振頻率區(qū)間區(qū)域一致，因此建議懸浮隧道的浮重比不要選擇過小，這樣可以更好地避免懸浮隧道發(fā)生渦激振動“鎖定”現(xiàn)象。

圖9為不同浮重比情況下懸浮隧道在不同自振頻率下的受力特性分析，從結果可以看出，浮重比對平均拖曳力系數(shù)和升力系數(shù)均方根影響不大，只是對個別頻率下的升力系數(shù)有較大影響，建議在懸浮隧道設計時結構遠離該自振頻率。

6 結論

本文應用有限元數(shù)值方法求解不可壓縮粘性流體雷諾平均Navier-Stokes方程，結合任意拉格朗日-歐拉(ALE)動網(wǎng)格方法，通過精確的CFD數(shù)值模擬，對不同雷諾數(shù)下懸浮隧道的渦激振動問題進行了研究，主要結論如下：

(1)通過網(wǎng)格敏感性分析證明數(shù)值模型的穩(wěn)定性并找到了兼顧計算效率和精度的網(wǎng)格，隨后通過模擬圓柱在高雷諾數(shù)下的渦激振動問題，與其他學者的實驗結果對比，證明本文所建立的模型可以用于高雷諾數(shù)流固耦合問題的研究。

(2)通過懸浮隧道在水流作用下的渦激振動研究，分析了在不同結構自振頻率下懸浮隧道的受力和運動情況，并研究了不同浮重比情況下懸浮隧道的受力和運動的結果。研究表明，雷諾數(shù)對發(fā)生“鎖定”的區(qū)間影響不大；對于不同浮重比的懸浮隧道，浮重比對振動幅值和受力特性影響不大，僅對個別頻率下的升力系數(shù)有較大影響，建議在懸浮隧道設計時結構遠離該自振頻率，但是對引起“鎖定”的區(qū)間有一定影響，當浮重比較小時，其發(fā)生“鎖定”的區(qū)間較大，因此建議設計時在保證錨索張力安全的情況下，盡量避免選取較小浮重比情況，既可以節(jié)省成本，又可以減小“鎖定”區(qū)間。