計及平臺垂蕩的立管渦激振動模擬與試驗驗證

袁昱超， 薛鴻祥， 唐文勇

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

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

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

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

1 渦激振動時域模型

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

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

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

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

Fhydro,y=Fv+Fd+Fm=

(1)

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

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

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

計算流體力載荷時，立管的振動幅值和頻率是最為關鍵的決定因素，因此，在每一分析步中，需提取各節點位移(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]作為鎖定區間，上述鎖定區間需要考慮真實環境及試驗條件下St的差異進行修正.修正方法為

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

(3)

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

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

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

(4)

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

1.3 結構動力響應求解

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

(5)

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

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

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

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

2 數值模型試驗驗證

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

2.1 試驗描述

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

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

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

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

2.2 固有頻率對比

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

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

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

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

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

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

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

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

2.4 響應幅值包絡線對比

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

3 聯合作用下結構響應特性研究

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

3.1 亞諧振響應

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

(6)

將Y寫為級數展開形式

Y=Y0+εY1+ε2Y2…

(7)

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

通過求解式(8)和式(9)可得各自穩態解

Y0=A0cosωCTt+B0sinωCTt

(10)

(11)

式中：

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

3.2 非對稱振型

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

3.3 Mathieu型共振

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

4 結論

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

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

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

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

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