變張力細長柔性立管渦激振動響應特性數值研究

高 云，潘港輝，劉 磊，柴盛林

（1.哈爾濱工業大學（威海）海洋工程學院，山東 威海 264209；2.安徽長江液化天然氣有限責任公司，安徽 蕪湖 241000；3.西南石油大學油氣藏地質與開發工程國家重點實驗室，成都 610500）

0 引 言

立管是連接海上浮式平臺與海底井口的裝置，在深海油氣工程中應用廣泛。當漩渦流經立管時，在一定的雷諾數范圍內會導致立管發生振動，這種現象被稱為渦激振動（VIV）[1]。渦激振動是一種非常典型的流固耦合問題[2]。一般認為，VIV 引起的動載荷會導致立管產生疲勞損傷甚至疲勞破壞[3-5]。由于立管處在水下，其維護費用很高，立管若發生疲勞失效，會造成巨大的經濟損失。因此，可靠地預測立管渦激振動響應是其結構設計過程中的關鍵問題之一。在過去的幾十年里，許多研究者對VIV問題進行了研究。由于立管VIV 振動具有非常復雜的流固耦合特性，早期的研究大多集中在剛性圓柱上。對剛性圓柱體的研究可以揭示振動問題的許多基本機理，如鎖定現象和滯回特性等[6-9]。

近年來，隨著深水油氣勘探的迅速發展，許多學者對大長徑比柔性結構的VIV響應特性展開了研究。根據研究方法，這些研究主要可分為物理實驗和數值模擬[10-16]。實驗方法雖然能得到豐富可靠的結果，但也存在一定的局限性，如模型規模有限、成本昂貴、復雜流剖面生成困難等。與實驗方法相比，數值方法可以克服其中的一些局限性，例如可以顯著降低研究成本，更容易實現復雜流剖面。數值計算方法基于流體力的確定方法可分為兩種。一種方法是計算流體動力學方法（CFD），它通過直接求解Navier-Stokes 方程來計算流體力；另一種方法是半經驗方法，它是根據經驗系數來確定流體力。與CFD方法相比，半經驗方法的顯著優點是計算量小。因此，為了滿足對細長柔性結構的VIV響應特性進行復雜多參數敏感性分析的要求，有必要提出一種能快速預報柔性結構VIV 響應特性的數值模型。近些年來，尾流振子法被廣泛應用于預測細長柔性結構的VIV響應特性，如表1所示[17-28]。

表1 部分基于尾流振子模型的柔性圓柱體VIV響應研究Tab.1 Selected VIV studies of a flexible cylinder using a wake oscillator model

由表1 中最后兩欄可以看出：針對柔性立管的研究，根據考慮軸向張力的特征，可分為恒張力研究和變張力研究。在變張力研究中，目前考慮的張力成分均比較單一：其中一部分研究只考慮了濕重，而另一部分研究則只考慮了彎曲振動。而在實際深海工程中，柔性立管軸向張力同時包括由濕重引起的等效張力和因結構彎曲振動而帶來的額外張力。此外，在實際深海工程中，真實海流既不是均勻流，也不是線性剪切流，而是階梯流[29-30]，這種海流的流剖面更接近指數剪切流。綜上，準確預報出真實海洋環境中變張力柔性立管的渦激振動響應特性，對海洋立管早期的合理設計和服役期間的安全工作均具有重要的理論和工程價值。

本文從實際海洋工程環境出發，建立變張力柔性立管渦激振動響應模型。考慮立管濕重和彎曲振動的作用計算軸向張力，采用二階中心差分法編程計算柔性立管的渦激振動響應。基于數值計算結果，系統地比較線性剪切、指數剪切以及真實階梯流剖面下柔性立管的VIV 響應特性，并進一步闡述三種流剖面作用下柔性立管渦激振動響應特性的特有規律。

1 模型描述

如圖1 所示，考慮一初始長度為L、外徑為D、內徑為d的細長柔性立管在來流作用下引起的橫流方向渦激振動響應問題。立管兩端采用鉸接邊界條件，取坐標原點O為立管的下端，X方向為順流方向，Z方向為鉛直方向，Y方向則是橫向振動方向，X、Y和Z三個方向形成右手直角坐標系。柔性立管上受到的張力為Θ（Z，T），柔性立管的彎曲剛度為EI，其中E和I分別是彈性模量和橫截面慣性矩。

圖1 線性剪切流作用下時變張力柔性立管模型Fig.1 Flexible riser with time-dependent varying tension subjected to a linear shear flow

本文中分別計算了三種不同類型的流剖面：線性剪切流、指數剪切流和真實階梯流剖面。這里為了研究方便，流剖面中最小流速取為最大流速的0.05 倍（即：Umin=0.05Umax）。其中線性剪切流剖面和指數剪切流剖面的計算公式可表示為

式中，U（Z）表示位于Z處流剖面的流速。對于真實階梯流剖面，通過最大流速乘以歸一化流剖面系數得到，該系數來自于Gao等[30]使用的值，如表2所示。將實際階梯流剖面沿立管軸線分為10個截面，每個截面的流速值呈線性變化特性，Nf為每段截面端點處的歸一化系數。由表2 可以看出，對于本文研究的階梯流，Umin=0.05Umax，這與前面介紹的線性剪切流和指數剪切流中最小流速與最大流速的關系相同。

表2 真實階梯流剖面歸一化系數Tab.2 Normalized coefficients for a real stepped flow

將圖1 中的柔性立管看作是細長張力梁模型，建立如下振動方程：

式中：rs和rf分別表示結構阻尼系數以及流體阻尼系數；mtotal為單位長度的振動系統質量，包括柔性立管結構質量、立管內部流體質量和立管外部流體附加質量三部分，可表示為

式中：ρs、ρf和ρw分別為柔性立管材料密度、立管內部流體密度和立管外部流體密度；CM為附加質量系數，對于圓柱體，CM=1.0。這里假設結構阻尼系數rs為0，流體阻尼系數可表示為rf=γΩfρD2，其中：Ωf為根據斯脫哈爾關系式計算得到的局部漩渦脫落頻率，Ωf=2πStU（Z）/D，St為斯脫哈爾數；γ為粘滯力系數，γ=--CD/4πSt，--CD為平均拖曳力系數，這里取為1.2。式（2）中，p（Z，T）為單位長度立管受到的升力，可表示為

式中，CL（Z，T）為升力系數，CL（Z，T）=CL0·q（Z，T）/2（CL0為靜止圓柱的升力系數，這里取為0.3；q（Z，T）表示圓柱體Z處在Y方向的尾流振子的運動）。式（2）中的Θ（Z，T）表示軸向張力，若同時考慮濕重帶來的張力變化以及結構彎曲振動帶來的額外張力變化，其表達式可寫為

式中，Θ（Z）表示僅考慮濕重時立管的軸向張力。此處引入頂張力系數K，假設立管上端所受張力為Ttop，那么Ttop可表示為Ttop=K×L×Wr，其中Wr為立管結構的濕重，可表示為

式中，g為重力加速度。由圖1 可知，坐標系原點位于結構下端，則任意坐標位置Z處的Θ（Z）可由立管上端張力和單位長度結構濕重表示為

式（5）中，ΔΘ（Z，T）為考慮彎曲變形帶來的額外張力，可表示為

式中，Ap為立管橫截面積（Ap=（D2-d2）π/4），E為楊氏彈性模量，S為結構發生彎曲變形后的瞬時長度。

在立管軸線上任取一個微元段dZ，若忽略順流X方向振動，只考慮橫流Y方向振動，則變形后的微元長度dS可表示為

針對整個立管，在某一瞬時立管的總長度S可由式（9）進行軸向積分得到，并進一步簡化為

聯立式（5）、（7）、（8）和（10），可得到同時考慮濕重以及彎曲變形情況下立管的軸向張力為

采用改進的Van der Pol方程來滿足尾流振子的非線性特征，表示如下：

式中，A和ε為經驗參數，A=12，ε=0.3。將式（2）和式（12）轉換為無量綱形式，令：

式中，Ωref為依據參考流速Uref計算得到的斯脫哈爾漩渦發放頻率，Ωref=2πStUref/D。將式（13）中三個等式分別代入式（2）和式（12）中，整理得到無量綱形式的結構振動方程以及尾流振子方程：

式中，ωf（z）表示流剖面參數，ωf（z）=Ωf/Ωref=U（z）/Uref。這里將流剖面中的最大流速Umax選作Uref。質量比μ、系統無量綱參數ML、無量綱單位長度濕重a、無量綱彎曲剛度b和無量綱變化張力c（z，t）可分別表示為

2 數值方法

假設柔性立管無量綱總長度L/D可劃分為M段，計算無量綱總時間ttotal可劃分為N段，那么，計算空間步長為Δz=L/（D×M），計算時間步長Δt=ttotal/N。被離散后的M+1 個空間點可表示為：z=zi（i=0，1，2，…，M）；被離散后的N+1 時間點可表示為：t=tj（j=0，1，2，…，N）。假設tn時刻zm位置處對應的無量綱參數y和q可表示為和，那么式中各偏導數項的二階精度差分格式可表示為

將式（17）中的差分格式代入式（14）和式（15）可得

式（18）中c（zm，tn）表示如下：

由式（16）可以看出，表達式Θ的第二項ΔΘ（z，t）含有積分項，這里為了簡化計算，將式（17）中的差分格式代入積分項，并對其進行累加求和處理，可得Θ（zm，tn）的表達式為

將式（18）～（19）中的迭代表達式形式作進一步簡化，轉化為顯性表達式形式表示如下：

假設y的初始條件（t=t0）位移和速度均為0，即y=?y/?t=0；q的初始條件設為一個波數的微幅擾動，且?q/?t=0，結合式（17）可得到：

將式（24）依次代入式（22）和式（23）中得到t=t1時刻y和q的值，表示如下：

至此求解得到t=t0以及t=t1時刻y和q的值。在t=t2時刻以后，需要使用到邊界條件。邊界條件設為兩端鉸接，即y在z=0時刻和L/D處的位移為0，彎矩為0，可表示為

當m=0和m=M時，需要用到位移為0的邊界條件，即

當計算m=1和m=M-1時，需要用到彎矩為0的邊界條件，由彎矩為0，聯立式（17）可得到

對于式（22），先賦予初始條件式（25）和邊界條件式（27），即對于t0時刻zm處的y，t1時刻zm處的y，tn時刻z0處的y，以及tn時刻zM處的y在求解開始前都為已知。當n≥1且0≤m≤M時，通過式（22）循環迭代求解y，需要注意當求解m=1 和m=M-1 這兩個位置的y時，需要用到式（28）彎矩為0 的邊界條件。當tn+1時刻zm處的y已知，便可依據式（23）求得tn+1時刻zm處的q，依此類推對式（22）和式（23）進行迭代求解。本文的數值方法已經在過去的研究中得到了充分的驗證［24-26］。因此，在這里并沒有再次對數值模型進行相關驗證。

3 結果與討論

在以往的研究中，我們主要研究了具有恒張力模型的柔性圓柱體VIV響應特性。然而，這里我們著重研究柔性立管在時變軸向張力下的VIV 振動響應特性。因此，這里使用的數值模型中軸向張力部分比以往多考慮了濕重和彎曲振動帶來的影響[24-26]。研究過程中使用的主要參數與Xu等（2017）[22]使用的參數值相同，如表3所示。

表3 立管模型參數Tab.3 Parameters of the riser model

3.1 渦激振動位移研究

圖2給出了三種不同線性剪切流剖面下細長柔性圓柱體的無量綱振動位移響應。三種線性剪切流剖面的最大流速Umax依次為1.5 m/s、2.0m/s 和2.5 m/s，最小流速均為最大流速的0.05 倍，即Umin=0.05Umax。圖2中，第一列給出了三種不同最大流速對應的線性流剖面，第二列給出了結構振動位移響應RMS值，第三列給出了結構振動位移響應最大值，第四列給出了結構振動位移響應包絡線。由圖中的第二列可看出：當最大流速為1.5 m/s、2.0 m/s 和2.5 m/s 時，結構振動主導模態依次為10 階、13 階和16 階，這說明結構振動主導模態階數隨最大流速的增加呈上升趨勢。由結構振動位移RMS 值可看出：在結構兩端，節點處與腹點處的位移值差別較大，說明此時結構振動位移響應由駐波占主導；與結構兩端區域相比，中間區域附近節點處與腹點處的位移值差別相對較小，這說明此時結構振動位移由行波占主導。由圖2中第四列的振動位移包絡線同樣可看出：所有節點處的振動位移均明顯不為零，說明沿整個軸線方向上，結構振動位移響應均呈現出明顯的行波特性。

圖2 線性剪切流時三種不同最大流速下細長柔性圓柱的均方根、最大位移值以及位移包絡線的變化規律Fig.2 RMS,maximum values of the VIV displacements and VIV displacement envelopes of a long flexible cylinder under three different Umax for a linear shear flow

圖3給出了三種不同線性剪切流剖面下某個穩態時間段內無量綱振動位移隨時間和空間的變化云圖。圖2 中第二列和第四列中所體現出的結構振動的行波響應特性在圖3 的云圖中得到了更為明顯的體現。云圖不僅能反應出結構振動的行波響應特性，更能非常方便地判斷出行波的傳播方向和傳播速度。如圖中白色箭頭所示，行波的傳播方向均為從流速較大的上端區域傳播到流速較小的下端區域。通過對比三種不同流速下的位移變化云圖可看出：在選取的某穩定無量綱時間段（無量綱時間段為20）內，隨著最大流速的增加，結構行波在結構軸線位置上傳播的距離逐漸減小，這說明隨著最大流速的增加，結構振動響應的行波傳播速度呈下降趨勢。

圖3 線性剪切流時三種不同最大流速下細長柔性圓柱的位移隨時間（t）和位置（z）的變化規律Fig.3 Displacement evolutions versus time and span location of a long flexible cylinder under three different Umax for a linear shear flow

圖4給出了三種不同指數剪切流剖面下細長柔性圓柱體的無量綱振動位移響應，圖5給出了三種不同指數剪切流剖面下某個穩態時間段內無量綱振動位移隨時間和空間的變化云圖。由圖4（a）的第二列、第四列和圖5看出：當最大流速為1.5 m/s時，在結構上端點附近，振動位移響應由駐波和行波共同主導；而在結構的中部以及下部區域，振動位移響應完全由駐波占主導，但同時具備輕微的行波特性。由圖4（b）和圖4（c）可看出：隨著流速的增加，當最大流速為2.0 m/s和2.5 m/s時，結構振動位移響應在整個軸線上均為行波占主導。圖6給出了三種不同真實階梯流剖面下細長柔性圓柱體的無量綱振動位移響應，圖7 給出了三種不同真實階梯流剖面下某個穩態時間段內無量綱振動位移隨時間和空間的變化云圖。結合圖6和圖7可明顯看出：對于最大流速較小的兩種工況（即：Umax=1.5 m/s和Umax=2.0 m/s 時），在結構上端區域（流速較大區域），結構振動位移響應由駐波和行波共同主導；在結構下端區域（流速較小區域），結構振動位移響應則由駐波占主導。而對于最大流速較大的工況（即Umax=2.5 m/s時），結構振動位移響應在沿整個軸線方向上均呈現出行波占主導的特性。

圖4 指數剪切流時三種不同最大流速下細長柔性圓柱的均方根、最大位移值以及位移包絡線的變化規律Fig.4 RMS,maximum values of the VIV displacements and VIV displacement envelopes of a long flexible cylinder under three different Umax for an exponential shear flow

圖5 指數剪切流時三種不同最大流速下細長柔性圓柱的位移隨時間（t）和位置（z）的變化規律Fig.5 Displacement evolutions versus time and span location of a long flexible cylinder under three different Umax for an exponential shear flow

圖6 真實階梯流時三種不同最大流速下細長柔性圓柱的均方根、最大位移值以及位移包絡線的變化規律Fig.6 RMS,maximum values of the VIV displacements and VIV displacement envelopes of a long flexible cylinder under three different Umax for a real stepped flow

圖7 真實階梯流時三種不同最大流速下細長柔性圓柱的位移隨時間（t）和位置（z）的變化規律Fig.7 Displacement evolutions versus time and span location of a long flexible cylinder under three different Umax for a real stepped flow

綜合比較圖2、圖4 和圖6 的第二列可看出：對于某個選定的最大流速，線性剪切流剖面下的結構振動位移響應RMS 值與指數剪切流剖面和真實階梯流剖面下的結構振動位移響應RMS 值均存在明顯的差異。對于同一個最大流速，指數剪切流剖面和真實剪切流剖面下的結構振動位移響應RMS 值非常接近，均比線性剪切流剖面下的結構振動位移響應RMS 值要小很多。由圖2、圖4 和圖6 的第二列同樣可以看出：對于某個選定的最大流速，三種不同流剖面下結構振動位移的主導模態階數非常接近，只有輕微的區別。以Umax=2.0 m/s 為例，當流剖面分別為線性剪切、指數剪切和真實階梯流剖面時，結構振動主導模態階數依次為13階、14階和14階。

3.2 結構振動響應的頻譜分析

圖8～10 依次給出了線性剪切、指數剪切和真實階梯流剖面下結構的振動位移時間歷程曲線、響應幅值譜和相圖。值得注意的是，對于每種不同類型的流剖面，分別選取了3 種最大流速（即Umax=1.5 m/s，2.0 m/s，2.5 m/s）加以研究。為了研究方便，這里僅對振動位移RMS 值的最大值加以研究。研究過程中選取的時間段為2000～2500 的振動穩態時間段（如圖中的第一列所示），對穩態段振動位移時間歷程曲線做快速傅里葉變換便可得到響應幅值譜（如圖中第二列所示）；為了更透徹地研究振動位移響應特性，這里進一步對振動位移相圖進行了分析（如圖中第三列所示）。由圖8可看出：當結構處在線性剪切流剖面中時，對于三種不同的流速工況，結構振動均呈現出明顯的概周期振動特性；振動響應呈現出明顯的寬帶分布特性，且存在多個明顯的峰值頻率；相圖軌跡由多個不可重復的橢圓組成，且各個橢圓的中心位置較為接近。

圖8 線性剪切流時三種不同最大流速下最大RMS值位置處VIV位移響應的時間歷程、振幅譜和相圖Fig.8 Time histories,amplitude spectra and phase portraits of the VIV displacement response at the maximum RMS value under three different Umax for a linear shear flow

由圖9 可以看出：對于指數剪切流剖面，當Umax=1.5 m/s 時，結構振動呈現出規律的周期性振動特性，結構振動僅存在一個峰值頻率，相圖軌跡為一系列接近重合的橢圓組成；隨著流速的增加，當Umax=2.0 m/s，2.5 m/s時，結構振動呈現出不規則的混亂特性，且振動響應頻率分布帶非常寬，分布在整個0～1 區間內。與圖9 中的線性剪切流剖面相比，類似的地方是相圖軌跡均由多個不可重復的橢圓組成，但不同之處在于線性剪切流剖面下相圖軌跡中不同橢圓的中心位置較為接近，而指數剪切流剖面下相圖軌跡中不同橢圓的中心位置差別很大。由圖10 可以看出：當結構處在真實階梯流剖面時，對于三種不同的流速工況，結構振動均呈現出明顯的概周期振動特性。

圖9 指數剪切流時三種不同最大流速下最大RMS值位置處VIV位移響應的時間歷程、振幅譜和相圖Fig.9 Time histories,amplitude spectra and phase portraits of the VIV displacement response at the maximum RMS value under three different Umax for an exponential shear flow

圖10 真實階梯流時三種不同最大流速下最大RMS值位置處VIV位移響應的時間歷程、振幅譜和相圖Fig.10 Time histories,amplitude spectra and phase portraits of the VIV displacement response at the maximum RMS value under three different Umax for a real stepped flow

4 結 論

本文對三種不同流剖面（線性剪切、指數剪切和真實階梯流剖面）、三種不同的流速工況（即Umax=1.5 m/s，2.0m/s，2.5 m/s）下的細長柔性結構渦激振動響應特性進行了數值研究。該數值研究可為柔性立管早期的合理設計和服役期間的安全工作提供理論支持和技術保障。基于上述數值研究結果，可以得到如下結論：

（1）對于線性剪切流剖面，在結構兩端結構振動位移由駐波主導，而在結構中間區域，結構振動位移則由行波主導；對于指數剪切流剖面和真實階梯流剖面，隨著最大流速的增加，振動位移響應沿整個軸線方向上由駐波主導，但具備輕微行波特性逐漸向完全由行波主導轉變。因此，立管設計過程中（尤其在高流速下），要重點關注立管軸線方向上的行波傳播對立管渦激振動響應特性的影響。

（2）當使用線性剪切流剖面去預測真實階梯流剖面時，得到的振動位移響應預測值與真實值存在明顯差別，而得到的主導模態階數預測值與真實值差別不大；當使用指數剪切流剖面去預測真實階梯流剖面時，得到的振動位移響應預測值以及主導模態階數預測值均與對應的真實值非常接近。

（3）當使用線性剪切流剖面去預測真實階梯流剖面時，得到的振動響應幅值以及振動頻率分布的預測值均與真實值存在很大差別；當使用指數剪切流剖面去預測真實階梯流剖面時，得到的振動響應幅值與真實值非常接近，但振動頻率分布的預測值與真實值仍存在較大差別。因此，在立管早期設計以及后期分析中，為了保證數值計算結果的可靠性，應盡量使用具有分段表達式的真實階梯流剖面來進行分析，而不是使用具有單一表達式的簡單流剖面來進行分析。