袁昱超， 薛鴻祥， 唐文勇

(上海交通大學 海洋工程國家重點實驗室； 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心， 上海 200240)

渦激振動(VIV)一直以來都是海洋細長結(jié)構(gòu)物設計過程中面臨的一個關鍵性問題，尤其對于頂張式立管的動力響應分析和疲勞壽命評估有著舉足輕重的影響.由于其極強的非線性，渦激振動在流固耦合領域仍然是一個極具挑戰(zhàn)性的研究課題.作為油氣生產(chǎn)系統(tǒng)的重要組成部分，頂張式立管是連接海上浮式平臺和海底井口的橋梁，扮演著將油氣資源由海底輸送至頂端平臺的媒介.影響立管動力響應的結(jié)構(gòu)剛度可分為由立管材料固有物理屬性決定的彎曲剛度和由立管軸向有效張力提供的附加橫向剛度2種成分，通常情況下，會在立管頂端施加一定的預緊力以提供必要的附加橫向剛度.在復雜的海洋環(huán)境中，波浪與海流的聯(lián)合作用易激發(fā)海上浮式平臺的垂蕩運動，平臺的升沉帶動張緊器的壓縮或拉伸，從而導致立管軸向張力呈現(xiàn)隨時間波動的變化特征.這不僅使得立管動力響應分析時需要引入動邊界條件，更為本質(zhì)的是結(jié)構(gòu)剛度將隨時間周期性變化.因此，計及頂端平臺垂蕩運動(即軸向張力時變效應)的頂張式立管渦激振動響應相較軸向張力恒定條件更接近于真實的海洋環(huán)境，也勢必更為復雜.

試驗研究和數(shù)值模擬是目前對渦激振動加深了解和認識的2類主要途徑.近年來，隨著相關試驗技術(shù)的不斷進步，F(xiàn)ranzini等[1]針對小尺度立管模型計及頂端垂蕩運動的渦激振動問題開展了相關試驗研究.時域數(shù)值模擬同樣是研究此類問題時一種可供選擇的有效方法，且不易受到模型尺寸的限制.時域渦激振動數(shù)值模型根據(jù)其依托的核心算法可分為計算流體力學、尾流振子模型和流體力分解模型3類.Josefsson等[2]基于二維切片理論借助計算流體力學方法研究了結(jié)構(gòu)彎曲剛度可變條件下圓柱體渦激振動響應.現(xiàn)有關于立管軸向張力時變(有時也被稱為參數(shù)激勵)條件下的渦激振動研究主要基于尾流振子模型[3-7].然而，已有的數(shù)值模擬研究鮮有與已公布的模型試驗結(jié)果進行對比，其模擬時變張力工況的有效性及預報精度仍有待進一步論證.

流體力分解模型同樣是被學術(shù)界和工程界廣泛用于預報頂張式立管渦激振動的一類較為成熟的方法，Lie[8]早在1995年就已開展了相關研究工作.在近十年間，Sidarta 等[9]和Wang等[10]相繼提出了各自的流體力分解模型用于開展立管在各種形式來流下渦激振動時域響應的相關研究.但上述流體力分解模型大多沒有計及頂端平臺垂蕩運動對立管結(jié)構(gòu)動力響應的影響，而是簡單地認為軸向張力在整個渦激振動過程中保持恒定.本文推薦了一套可供選擇的流體力分解模型，能夠用于預報張力時變條件下頂張式立管渦激振動時域響應.渦激振動流體力載荷基于受迫振蕩試驗數(shù)據(jù)，結(jié)構(gòu)剛度矩陣根據(jù)實時軸向張力沿管長方向的分布在每一分析步中進行更新.Franzini等[1]公布的模型試驗被用于驗證推薦數(shù)值模型在恒定及時變張力條件下的有效性.數(shù)值模擬與試驗觀測對比研究中發(fā)現(xiàn)，時變張力與渦激振動聯(lián)合作用下立管結(jié)構(gòu)響應呈現(xiàn)亞諧振成分、非對稱振型及Mathieu型共振等特性，本文針對以上問題作進一步展開討論，并給出了合理解釋.

1 渦激振動時域模型

典型的平臺-頂張式立管-海床系統(tǒng)可簡化為如圖1所示布置形式.圖中：Cartesian坐標體系用以描述推薦渦激振動時域模型；x軸順應來流方向；y軸垂直于流速方向；z軸為立管軸向即豎直向上方向.

圖1 平臺-頂張式立管-海床系統(tǒng)示意圖Fig.1 Sketch of platform-TTR-seabed system

1.1 渦激振動流體力載荷模型

基于渦激振動流體力分解模型，流體力載荷

Fhydro,y=Fv+Fd+Fm=

(1)

式中：Fv為激勵力；Fd為阻尼力；Fm為慣性力；Cv為激勵力系數(shù)；A為響應位移幅值；D為立管直徑；f為響應頻率；V為遭遇流速；t為時間；ρf為流體密度；cf為水動力阻尼系數(shù)；Ca為附加質(zhì)量系數(shù)，本文假定Ca為經(jīng)驗值 1.0.

由式(1)可知，渦激振動流體力載荷系數(shù)Cv和cf的確定是流體力分解模型的核心內(nèi)容，其在每一分析步內(nèi)均同時受控于響應幅值和頻率.圖2所示為本文選取的激勵力系數(shù)云圖[11].圖中：粗實線標注出對應Cv=0的能量輸入及輸出邊界；當Cv<0時，激勵力將遵循下文規(guī)定法則轉(zhuǎn)化為水動力阻尼；當超出圖2所示的試驗數(shù)據(jù)范圍時，同樣認為激勵力將依據(jù)下文提及的經(jīng)驗阻尼模型轉(zhuǎn)化為水動力阻尼.上述受迫振蕩試驗工況測得的斯坦頓數(shù)(St)約為 0.193.

圖2 渦激振動激勵力系數(shù)庫[11]Fig.2 VIV excitation force coefficient database [11]

計算流體力載荷時，立管的振動幅值和頻率是最為關鍵的決定因素，因此，在每一分析步中，需提取各節(jié)點位移(y)、速度(v)和對應時刻(t)求解振動幅值和頻率分別為

A=|yb-ya|/2,f=1/[2(tb-ta)]

(2)

式中：(ya,yb)和(ta,tb)分別為相鄰兩個瞬時速度為零(即v=0)點a和b對應的位移和時刻.

1.2 鎖定判定準則

鎖定(Lock-in)是渦激振動問題中尤其需要引起重視的響應特性，本文參照文獻[12]中雙自由度圓柱體受迫振蕩試驗研究成果，選取無因次頻率帶寬[0.125，0.25]作為鎖定區(qū)間，上述鎖定區(qū)間需要考慮真實環(huán)境及試驗條件下St的差異進行修正.修正方法為

(fr/St)test=(fr/St)actual

(3)

式中：fr為無因次頻率且有fr=fD/V.

若無因次頻率落在渦激振動激發(fā)帶寬內(nèi)，鎖定將發(fā)生.本文假定立管單元的無因次頻率將被鎖定到渦激振動激勵中心對應的 0.17 (見圖2)以確保激勵力系數(shù)取到峰值，此處的 0.17 同樣需要根據(jù)式(3)進行修正，此時對應的頻率則作為渦激振動主導頻率用于根據(jù)圖2計算激勵力系數(shù).

水動力阻尼同樣是渦激振動數(shù)值模擬領域較為重要的一個課題，當激勵力系數(shù)為正時，激勵力與立管速度同向，整個系統(tǒng)內(nèi)能量由流體傳入立管，而當激勵力系數(shù)為負時，水動力阻尼力作用于立管，能量也將由立管傳回流體.參考渦激振動頻域方法中的處理手段，本文假定流體傳入立管的等效能量在一個振動周期內(nèi)傳回流體，根據(jù)能量守恒原理，水動力阻尼系數(shù)計算公式可表示為

(4)

若無因次頻率位于激發(fā)帶寬外，則可根據(jù)圖2求得Cv<0的激勵力系數(shù)，進而基于式(4)得到對應的水動力阻尼系數(shù).當立管響應幅值和頻率超出圖2所示試驗數(shù)據(jù)范圍，采用Venugopal提出的渦激振動經(jīng)驗阻尼模型[13]模擬水動力阻尼力，其中的經(jīng)驗系數(shù)Chf、Csw和Clf推薦選取 0.18、0.25 和 0.2.

1.3 結(jié)構(gòu)動力響應求解

頂張式立管渦激振動響應通常遵循小位移假定，立管模型可認為是滿足Euler-Bernoulli梁理論的抗彎彈性結(jié)構(gòu)物，其在橫流向的運動控制微分方程

(5)

式中：m為立管單位長度質(zhì)量；cs為結(jié)構(gòu)阻尼系數(shù)且有cs=2mωξ，ξ為結(jié)構(gòu)阻尼比，ω為圓頻率；E為彈性模量；I為慣性矩；T0為恒定軸向有效張力；ΔT、PT和φ分別為時變張力的幅值、周期和初始相位.

式(5)左側(cè)的第3和4項分別對應由立管固有屬性EI決定的彎曲剛度和由軸向有效張力提供的附加橫向剛度.為了更好地說明本文的推薦方法，圖3給出渦激振動時域分析中每一分析步的完整流程圖.

圖3 渦激振動時域分析流程圖Fig.3 Flow-chart of VIV time domain analysis

本文推薦的渦激振動時域數(shù)值模型基于直接積分動力分析，采用Hilber-Hughes-Taylor (HHT)方法[14]求解式(5)中的結(jié)構(gòu)運動控制方程.HHT法作為一種隱式數(shù)值求解方法已較為成功地應用于經(jīng)典動力學問題的仿真模擬，由于其具有更好的收斂性、計算精度和效率，也被認為是一種改進的Newmark-β法.

2 數(shù)值模型試驗驗證

本文選取由文獻[1]中公布的一組小尺度立管模型(未拉伸長度L0為2.552 m)計及頂端垂蕩運動的渦激振動試驗結(jié)果用以驗證推薦數(shù)值模型預報恒定及時變張力工況的有效性.

2.1 試驗描述

試驗立管模型采用豎直放置形式，部分浸沒于水中，模型頂部由伺服電機裝置施加垂蕩幅值At/L0=1%的簡諧運動.沿立管跨長分布43個測量點，測量點位移通過水下相機借助光學追蹤技術(shù)直接獲取.頂端布置測壓元件用以測量軸向張力.圖4給出試驗裝置布置示意圖.表征立管模型物理特性的相關參數(shù)如表1所示.

結(jié)合上述信息，可估算出該試驗模擬的軸向張力時變工況下ΔT/T0約為 0.51.除一個恒定張力工況外，文獻[1]中僅公布了2個垂蕩頻率(即ωT/ωCT=2,3)下試驗觀測結(jié)果.其中：ωT為軸向張力變化圓頻率；ωCT為恒定張力條件下由渦激振動誘發(fā)的立管響應圓頻率.3個工況對應的雷諾數(shù)(Re)介于[700, 6 100]，本文選取實際St近似為 0.21.

圖4 模型試驗裝置布置圖[1]Fig.4 Experimental facilities arrangement of the model test[1]

模型參數(shù)數(shù)值未拉伸長度/m2.552拉伸長度/m2.602浸沒長度/m2.257外徑/m0.022彎曲剛度/(N·m2)0.056質(zhì)量比3.48浸沒質(zhì)量/kg2.05定常預張力/N40

2.2 固有頻率對比

自由衰減試驗表明立管模型n階固有頻率fN,n為 0.84nHz[1].本文基于FEM方法進行模態(tài)分析得到模型固有頻率預報值，與實測值對比情況如圖5所示.計算所得第1、2、3階固有頻率分別為 0.85、1.70 和 2.57 Hz，與試驗數(shù)據(jù)吻合較好.而從第4階開始，計算值略大于實測值.事實上，對于一端受拉的簡支梁，在立管固有彎曲剛度和軸向張力提供的附加橫向剛度共同作用下，結(jié)構(gòu)各階固有頻率理論上應呈現(xiàn)預報結(jié)果所示的增幅逐階遞增規(guī)律，而自由衰減試驗揭示的嚴格線性關系更像是僅考慮附加橫向剛度所得結(jié)果.由于本文涉及的3個工況下，主激發(fā)模態(tài)最高為3階，因此固有頻率對比所示誤差對下文結(jié)構(gòu)動力響應研究影響甚微.

圖5 立管模型固有頻率Fig.5 Natural frequencies of the riser model

2.3 響應位移時歷與幅值譜對比

3個計算工況下監(jiān)測點(z/L0=0.43)處立管結(jié)構(gòu)響應的預報結(jié)果與試驗觀測對比情況如圖6所示.圖6(a)表征恒定張力工況，圖6(b)和(c)分別表征ωT/ωCT=2,3的時變張力工況.其中：VR,1為約化速度，可通過VR,1=V/(fN,1D)計算得到.圖6(a)表明，均勻流條件下恒定張力工況中立管渦激振動響應展現(xiàn)出明顯的單一模態(tài)激發(fā)特征，位移時歷嚴格服從正弦變化規(guī)律，預報和實測響應位移幅值均大致穩(wěn)定在 0.6D，預報值略大于實測值.預報所得主導頻率同樣與試驗數(shù)據(jù)吻合，均接近于結(jié)構(gòu)一階固有頻率，由此可知，約化速度VR,1=5.63 時，渦激振動穩(wěn)定地激發(fā)該立管模型一階固有模態(tài).

圖6 監(jiān)測點處響應位移時歷曲線及對應幅值頻響譜Fig.6 Time histories of response displacement and corresponding amplitude spectra at monitor point

由圖6(b)可知，變張力使得預報與實測結(jié)果中的響應位移幅值放大約1倍.試驗觀測中捕捉到位移振幅劇烈調(diào)制的現(xiàn)象且調(diào)制過程無明顯的規(guī)律性或周期性，而對應的預報結(jié)果中雖然也存在一定的振幅調(diào)制現(xiàn)象，但就整個時間歷程而言位移響應表現(xiàn)得更為穩(wěn)定.盡管如此，預報所得位移幅值頻響譜仍與實測結(jié)果保持較高的一致性，能量響應集中在ω/ωCT=1處并伴隨著少量成分出現(xiàn)在ω/ωCT=2, 3, 4, 5等處，兩者之間細微的區(qū)別在于實測結(jié)果中其他頻率成分弱于預報結(jié)果所示，且如圖中虛圓圈標注所示，預報結(jié)果中ω/ωCT=3,5的能量響應較之ω/ωCT=2,4更強.

由圖6(c)可知，預報和實測位移幅值均在0.6D附近輕微波動，而對應的幅值頻響譜中揭示出有別于前2個工況的新現(xiàn)象.預報和實測所得頻響譜在ω/ωCT=1，3處均出現(xiàn)突出的能量響應峰值，且ω/ωCT=3處的能量峰值大于ω/ωCT=1處對應值，即立管動力響應不再主要展現(xiàn)一階固有模態(tài)而變?yōu)橛山Y(jié)構(gòu)3階固有模態(tài)主導.一些其他的頻率成分同樣出現(xiàn)在預報所得頻響譜中，而實測頻響譜似乎并未發(fā)現(xiàn)類似的亞諧振響應.Franzini等[15]指出，根據(jù)模型試驗研究成果，在ωT/ωCT=3的純頂端垂蕩運動條件(即流速為0，無渦激振動響應)下發(fā)現(xiàn)了穩(wěn)定的橫向振動響應，但僅限于第3階模態(tài)而對于第1、2階模態(tài)并不適用.由此認為，在VR,1=5.63 這一約化速度下第3階模態(tài)響應伴隨著第1階渦激振動同時發(fā)生，對于造成這一特殊現(xiàn)象的深層次原因，文獻[15]中并未提及.上述事實也在一定程度上解釋了為何ωT/ωCT=3的時變張力工況中體現(xiàn)的頻響特征與前2個工況有所區(qū)別.

圖7 響應位移幅值沿管長包絡線Fig.7 Amplitude envelops of response displacement along the riser

2.4 響應幅值包絡線對比

圖7給出上述3個工況下預報和實測所得響應位移幅值沿管長包絡線對比情況.圖中：L為立管拉伸長度.對于恒定張力工況，預報所得響應位移幅值包絡線與試驗觀測較為吻合，均為規(guī)則的1階振型，僅在立管上半部分略大于實測值.兩者振幅峰值均出現(xiàn)在z/L=0.42 附近，偏離立管中點.相比之下，在ωT/ωCT=2的時變張力工況中，幅值包絡線被顯著放大(1倍甚至更多).Franzini等[1]推斷這一放大現(xiàn)象是由于結(jié)構(gòu)發(fā)生了由Mathieu不穩(wěn)定引起的參數(shù)共振.實測結(jié)果不再是標準的1階模態(tài)振型，尤其是在z/L位于[0.25, 0.6]范圍內(nèi)，一些額外模態(tài)明顯參與立管渦激振動.預報結(jié)果同樣顯示出高階模態(tài)被激發(fā)的響應特征，在z/L=0.4 附近(即1階振型極值出現(xiàn)位置)最大響應幅值略大于實測值，呈現(xiàn)1階模態(tài)主導3階模態(tài)參與的整體振型，與圖6(b)中幅值頻響譜所示頻響特性對應.而對于ωT/ωCT=3的時變張力工況，預報和實測所得響應幅值包絡線吻合較好，均展現(xiàn)出明顯的3階模態(tài)振型，預報值僅在z/L=[0, 0.3] 范圍內(nèi)略大于試驗數(shù)據(jù).預報與實測振幅峰值出現(xiàn)位置同樣偏離標準3階模態(tài)振型對應位置，且兩者振幅包絡線均呈現(xiàn)一定程度的非對稱性.對比恒定張力工況，ωT/ωCT=3的時變張力對結(jié)構(gòu)響應同樣具有放大效應，只是不如ωT/ωCT=2條件下劇烈.

3 聯(lián)合作用下結(jié)構(gòu)響應特性研究

本節(jié)將就前文對比研究中發(fā)現(xiàn)的時變張力與渦激振動聯(lián)合作用下結(jié)構(gòu)動力響應展現(xiàn)的亞諧振響應、非對稱振型及Mathieu型共振等特性展開討論，探究其深層次原因.

3.1 亞諧振響應

由圖6(b)中幅值頻響譜所示，除渦激振動鎖定頻率外，ω/ωCT=3，5處同樣存在著明顯的能量響應，而這兩處恰好對應ωCT+ωT和ωCT+2ωT.事實上，時變張力與渦激振動聯(lián)合作用下ωCT±kωT的亞諧振頻響成分的確可能被激發(fā)，現(xiàn)基于解析推導給出相應的理論解釋.式(5)中的水動力阻尼主要起到限制立管振動幅值的作用，對于結(jié)構(gòu)共振條件的推導影響較小，現(xiàn)忽略水動力阻尼項cf?y/?t并假定y(z,t)=φ(z)Y(t)對動力學方程進行解耦.其中：φ(z)為模態(tài)振型函數(shù)；Y(t)為幅值函數(shù).對于兩端簡支梁可表示為sin(iπz/L)，式(5)可整理為

(6)

將Y寫為級數(shù)展開形式

Y=Y0+εY1+ε2Y2…

(7)

將式(7)代入式(6)，可得到關于ε的各階方程.出于簡潔考慮，僅以第0階和第1階表達式為例進行推導：

通過求解式(8)和式(9)可得各自穩(wěn)態(tài)解

Y0=A0cosωCTt+B0sinωCTt

(10)

(11)

式中：

關于ε的一階方程的穩(wěn)態(tài)解表明，額外的頻率成分ωCT±ωT滿足諧振條件.若求解關于ε的k階方程，參照上述理論推導不難得出在渦激振動與時變張力聯(lián)合作用下亞諧振頻率成分ωCT±kωT的確可能被激發(fā)的結(jié)論.

3.2 非對稱振型

由圖7可知，3個工況(尤其是時變張力工況)下立管整體振型均呈現(xiàn)一定的非對稱性，主要表現(xiàn)為振幅峰值出現(xiàn)位置偏離理論值且立管底端與頂端振型明顯不對稱.本文認為上述現(xiàn)象主要是由于該試驗立管模型采用豎直放置形式，立管結(jié)構(gòu)自重使得軸向張力定常成分沿管長(即垂向)分布不均勻，自上而下線性減小.因此，各分段對整體附加橫向剛度矩陣的貢獻也不相同.事實上，豎直放置立管的固有模態(tài)滿足具有非正交解的Bessel類方程，模態(tài)振型峰值出現(xiàn)位置的確會向底端偏移.另外，T0(z) 沿管長不等，ΔT/T0自上而下逐漸變大，因此時變張力對立管底端的影響效應也會比頂端更為明顯.Chen等[7]研究表明，隨著ΔT/T0的增加，時變張力對結(jié)構(gòu)響應的放大效應也越發(fā)顯著.這也能夠解釋為何底端響應幅值大于頂端.文獻[4]中采用尾流振子模型同樣得出與本文預報結(jié)果趨勢相同的非對稱整體振型，一定程度上也佐證了本文數(shù)值模型的正確性.立管模型部分浸沒水中導致的附加質(zhì)量沿管長分布不均也會對結(jié)構(gòu)整體振型產(chǎn)生一定影響，但由于該試驗中浸沒比達到87%，本文認為該影響效應與上述兩點相比是次要的.

3.3 Mathieu型共振

圖7(b)表明，ωT/ωCT=2的時變張力工況下的位移振幅包絡線相較恒定張力工況放大1倍甚至更高，有理由認為此時的時變張力與已有的渦激振動協(xié)同作用下，立管結(jié)構(gòu)發(fā)生強烈共振，出于結(jié)構(gòu)安全性的考慮，共振狀態(tài)是立管在工程設計階段希望極力避免的，本文對于此類共振現(xiàn)象查閱相關文獻給出內(nèi)在機理解釋.事實上，ωT=2ωCT的時變張力正是觸發(fā)結(jié)構(gòu)Mathieu失穩(wěn)的最危險條件，文獻[15]中試驗研究表明僅ωT=2ωCT的時變張力即可導致結(jié)構(gòu)發(fā)生參數(shù)共振.這種由于Mathieu不穩(wěn)定造成的結(jié)構(gòu)共振伴隨著渦激振動鎖定現(xiàn)象的發(fā)生被進一步放大，最終導致結(jié)構(gòu)動力響應顯著加劇.因此，本文認為此類時變張力與渦激振動聯(lián)合作用下發(fā)生的共振現(xiàn)象可定義為Mathieu型共振.張力變化頻率在ωT=2ωCT附近均可能誘發(fā)Mathieu型共振，故ωT/ωCT=3的時變張力同樣對立管振動幅值具有一定的放大效應.

4 結(jié)論

計及頂端平臺垂蕩運動的頂張式立管渦激振動問題仍然是現(xiàn)階段值得深入研究的前沿課題.本文基于渦激振動流體力分解模型提出一套可用于預報頂張式立管在海流及時變軸向張力聯(lián)合作用下結(jié)構(gòu)動力響應的時域數(shù)值方法.在對其恒定及時變張力條件下的有效性進行試驗驗證的基礎上，詳細討論了對比研究中發(fā)現(xiàn)的若干響應特性，并得到以下結(jié)論：

(1) 基于本文數(shù)值方法的預報結(jié)果與試驗觀測吻合較好，推薦數(shù)值模型可有效預報軸向張力恒定及時變條件下頂張式立管渦激振動時域響應.

(2) 通過理論推導與數(shù)值模擬均可發(fā)現(xiàn)，當計及軸向張力時變效應時，除渦激振動鎖定頻率外，對應ωCT±kωT的亞諧振頻率成分也可能被激發(fā).

(3) 立管自重導致軸向張力定常成分沿管長分布不均勻，加之張力時變效應的影響，立管結(jié)構(gòu)響應呈現(xiàn)非對稱振型.

(4) 時變軸向張力與渦激振動聯(lián)合作用下，結(jié)構(gòu)會在張力變化頻率ωT=2ωCT附近發(fā)生Mathieu型共振，響應位移幅值將明顯放大.