朱 磊，李鳳明，陳永訢

（1.哈爾濱工程大學(xué)航天與建筑工程學(xué)院，哈爾濱 150001；2.海洋石油工程股份有限公司，天津 300461）

0 引 言

立管是連接海面平臺(tái)與海底井口的重要通道，其渦激振動(dòng)問題長(zhǎng)期受到工程領(lǐng)域的關(guān)注。在一定來流作用下，立管兩側(cè)產(chǎn)生周期性的旋渦瀉放，從而產(chǎn)生周期性的脈動(dòng)壓力，引發(fā)立管結(jié)構(gòu)的劇烈振動(dòng)，在復(fù)雜多變的海洋載荷作用下，立管表現(xiàn)出復(fù)雜的振動(dòng)特性，尤其是與海洋平臺(tái)相連接的立管，海洋平臺(tái)的大幅垂蕩與立管自身渦激振動(dòng)相互耦合，對(duì)立管的安全生產(chǎn)造成了嚴(yán)重影響[1-2]。

在復(fù)雜的海洋環(huán)境條件下，海洋平臺(tái)垂蕩和立管渦激振動(dòng)的耦合振動(dòng)響應(yīng)，表現(xiàn)出的振動(dòng)形式更加復(fù)雜。這一問題，尤其是立管結(jié)構(gòu)的穩(wěn)定性，引起很多國(guó)內(nèi)外學(xué)者的關(guān)注。Chatjigeorgiou 和Mavrakos[3-4]分別采用多尺度方法和數(shù)值方法分析了細(xì)長(zhǎng)柔性結(jié)構(gòu)的穩(wěn)定性問題，并討論了阻尼對(duì)穩(wěn)定性的影響，結(jié)果表明，細(xì)長(zhǎng)柔性結(jié)構(gòu)失穩(wěn)主要出現(xiàn)在一階不穩(wěn)定區(qū)，高階失穩(wěn)可能性較小，阻尼可以顯著減小不穩(wěn)定區(qū)的大?。籗inir[5]研究了微彎曲管道輸送流體時(shí)的非線性振動(dòng)，給出了不同振幅的曲率函數(shù)和流體速度的分岔圖，在微彎曲管道輸送流體的橫向振動(dòng)中，觀察到了周期運(yùn)動(dòng)和混沌運(yùn)動(dòng)現(xiàn)象；Keber 和Wiercigroch[6]研究了結(jié)構(gòu)非線性對(duì)二維輸油管道渦激振動(dòng)響應(yīng)的影響，指出內(nèi)部流體流動(dòng)可以加強(qiáng)結(jié)構(gòu)非線性對(duì)管道動(dòng)力學(xué)響應(yīng)的影響；Kuiper 等[7]采用Floquet 理論，分別討論了參激振動(dòng)穩(wěn)定性問題和劇烈海況下立管底部出現(xiàn)負(fù)張力從而產(chǎn)生局部屈曲而導(dǎo)致的失穩(wěn)問題。除穩(wěn)定性問題之外，He 等[8]采用Galerkin 離散的方法研究了平臺(tái)橫向振動(dòng)作用下立管的振動(dòng)特性；唐友剛等[9]建立了柔性立管參激—渦激耦合振動(dòng)的計(jì)算模型，計(jì)算了剪切流下立管振動(dòng)響應(yīng)特性；陳偉民等[10]考慮參激—渦激耦合效應(yīng)，采用有限元方法離散立管模型，結(jié)果表明，立管振動(dòng)響應(yīng)幅值隨振動(dòng)模態(tài)階數(shù)的降低而增大，同時(shí)出現(xiàn)模態(tài)轉(zhuǎn)換現(xiàn)象。

可以看出，上述研究分別針對(duì)頂張力立管順流向和橫流向的非線性動(dòng)力學(xué)響應(yīng)，以及變張力作用下頂張力立管橫流向渦激振動(dòng)響應(yīng)的影響進(jìn)行研究。但是，將兩者結(jié)合研究的文獻(xiàn)很少，本文結(jié)合兩方面的工作，將目前研究變張力對(duì)頂張力立管橫流向渦激振動(dòng)響應(yīng)的影響拓展為二維問題，豐富了變張力對(duì)頂張力立管順流向渦激振動(dòng)響應(yīng)的研究?jī)?nèi)容，且以往文獻(xiàn)對(duì)于立管渦激振動(dòng)響應(yīng)的預(yù)報(bào)大多采用有限元法，本文將采用Hamilton 原理和Galerkin 方法建立結(jié)構(gòu)運(yùn)動(dòng)方程，避免了采用有限元法的繁瑣。

基于上述分析，本研究針對(duì)海洋平臺(tái)垂蕩運(yùn)動(dòng)，將其簡(jiǎn)化為軸向變化的頂張力，同時(shí)以Euler-Bernoulli 梁模型對(duì)立管結(jié)構(gòu)進(jìn)行建模，采用Hamilton 原理建立頂張力立管結(jié)構(gòu)的二維非線性運(yùn)動(dòng)方程，并對(duì)運(yùn)動(dòng)方程進(jìn)行無量綱化，采用Galerkin 方法對(duì)立管和渦尾跡的非線性方程進(jìn)行離散，并用四階Runge-Kutta方法進(jìn)行求解，對(duì)可變頂張力立管的二維非線性渦激振動(dòng)特性進(jìn)行研究。通過與已有實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證本文方法和數(shù)值計(jì)算結(jié)果的可行性和正確性。

1 建立結(jié)構(gòu)運(yùn)動(dòng)方程

1.1 結(jié)構(gòu)模型

研究中的立管模型如圖1所示。定義了一個(gè)笛卡爾坐標(biāo)系，x軸為沿立管軸向，y軸為沿來流速度方向（IL 方向），z軸為垂直于來流速度方向（CF 方向）。該系統(tǒng)由一個(gè)兩端鉸接柔性立管組成，輸送內(nèi)部流體，柔性立管系統(tǒng)置于外部流體中。

采用Hamilton原理和Galerkin方法建立立管的非線性運(yùn)動(dòng)方程。Hamilton原理表示為

式中，Ttotal和Vtotal分別表示系統(tǒng)的總動(dòng)能和總勢(shì)能，WI為外力做的功，t為時(shí)間。

系統(tǒng)的總動(dòng)能Ttotal包括空立管的動(dòng)能T1和內(nèi)部流體的動(dòng)能T2，它們分別表示為

式中：mr為單位長(zhǎng)度立管的質(zhì)量；ma為單位長(zhǎng)度附加流體質(zhì)量，ma=CaπD2ρ0/4，D為立管的外徑，ρ0為管外流體密度，Ca為附加流體質(zhì)量系數(shù)，取Ca=1；mi為單位長(zhǎng)度內(nèi)部流體的質(zhì)量；V和W分別為立管沿y軸和z軸的位移；L是立管的長(zhǎng)度；Ui為管內(nèi)流體流速。

系統(tǒng)的總勢(shì)能Vtotal由立管的彎曲勢(shì)能V1、重力勢(shì)能V2和軸向拉伸勢(shì)能V3構(gòu)成，它們分別表示為

圖1 立管示意圖Fig.1 Schematic of riser

式中，EI和EAr分別為立管的彎曲剛度和拉伸剛度，A為立管排開水的截面面積，Ar為立管管材的截面面積，g為重力加速度，Ttension為軸向拉伸引起的非線性張力，Te為平臺(tái)頂張力，εs為立管的應(yīng)變。

由于本文研究的重點(diǎn)是立管的非線性動(dòng)力學(xué)響應(yīng)，因此需要考慮立管軸向應(yīng)變引起的幾何非線性。根據(jù)Jiang等[11]的工作，由軸向拉伸引起的非線性張力表示為

立管的應(yīng)變?chǔ)舠表示為

式中，O（∈4）（∈＜＜1）為高階小量。

外力做的虛功為

式中，cs為結(jié)構(gòu)阻尼，cf為流體附加阻尼，F(xiàn)y和Fz為管外流體作用在立管上的水動(dòng)力。

將式（2）～（8）代入到式（1）中，利用變分原理，立管的IL 方向和CF 方向耦合的非線性偏微分運(yùn)動(dòng)方程可表示為[12]

式中，m為單位長(zhǎng)度系統(tǒng)的總質(zhì)量，m=mr+mi+ma。結(jié)構(gòu)阻尼cs和流體附加阻尼cf可表示為[12]

式中：ζ為結(jié)構(gòu)阻尼比，取ζ=0.05；ωr為立管振動(dòng)的角頻率；?為水動(dòng)力阻尼系數(shù)，根據(jù)Skop 等[13-14]的研究，取?=0.8；ωs為旋渦脫落的角頻率，可表示如下：

式中：Ue為來流速度；St為斯特勞哈爾數(shù)（Strouhal number），取St=0.2。

頂端海洋平臺(tái)垂蕩運(yùn)動(dòng)作用下立管頂張力可表示為[9]

式中，T0為預(yù)張力，T為可變張力幅值，?為可變張力的頻率(即海洋平臺(tái)垂蕩的頻率)。

本模型的簡(jiǎn)支邊界條件可表示為[15]

1.2 水動(dòng)力和尾流振子模型

式（9）和式（10）中的Fy和Fz分別為IL 方向和CF 方向作用在立管上的水動(dòng)力。對(duì)于固定圓柱結(jié)構(gòu)，F(xiàn)y和Fz分別與拖拽力fD和渦激升力fL相關(guān)，fD和fL的作用方向與IL和CF方向一致。然而，當(dāng)立管在水動(dòng)力作用下發(fā)生振蕩時(shí)，拖拽力fD和渦激升力fL不再完全作用在y軸和z軸方向上。Fy和Fz可以表示為[16]

式中，fD可以看作是尾渦脫落的脈動(dòng)拖拽力，fD0是作用在立管上的平均拖拽力，fL是作用在立管上的渦激升力。

此外，fD、fD0、fL可以分別進(jìn)一步表示為[16]

式中，CDi是尾渦脫落的拖拽力系數(shù)，CD0是靜止圓柱體的平均拖拽力系數(shù)，CL是升力系數(shù)。

立管漩渦脫落區(qū)尾流特性本文采用的是基于van der Pol方程的尾流振子模擬，該方法是一種半經(jīng)驗(yàn)?zāi)Ｊ椒椒?，不考慮具體的流場(chǎng)結(jié)構(gòu)，而將流體和其中的振蕩物體視為一個(gè)整體系統(tǒng)，把尾流看成是一個(gè)非線性振子，尾流振子的振動(dòng)引起結(jié)構(gòu)的振動(dòng)，反過來，結(jié)構(gòu)的振動(dòng)又對(duì)尾流有一個(gè)反饋的作用。

根據(jù)Facchinetti等[16]和Srinil等[17-18]的研究，van der Pol尾流振子模型方程如下：

式中，εy、εz、Ay和Az是由實(shí)驗(yàn)得到的無因次參數(shù)，p(x,t)和q(x,t)分別表示IL和CF方向的尾流變量。

通過公式p(x,t)= 2CDi/CDi0和q(x,t)= 2CL/CL0，結(jié)構(gòu)模型與尾流振子模型建立了聯(lián)系，且Fy和Fz整理后可以進(jìn)一步表示為

式中，根據(jù)文獻(xiàn)[17]，取CDi0=0.1，CD0=2.0，CL0=0.3。

Wang等[12]的研究中，討論了三維柔性管振動(dòng)中參數(shù)的選取，當(dāng)四個(gè)無量綱參數(shù)設(shè)置為εy=0.02，εz=0.04，Ay=96和Az=12時(shí)，理論結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合較好。需要注意，van der Pol尾流振子模型是半經(jīng)驗(yàn)半解析的，該經(jīng)驗(yàn)?zāi)Ｐ秃雎粤藴u型和渦強(qiáng)等細(xì)節(jié)。然而，用van der Pol 模型預(yù)測(cè)的立管渦激振動(dòng)響應(yīng)與直接數(shù)值模擬預(yù)測(cè)的結(jié)果吻合較好。

1.3 Galerkin離散化

借助于下列無因次量：

式（9）、式（10）、式（19）、式（20）的耦合方程可以寫成以下無量綱形式：

為了離散立管和尾流的無量綱運(yùn)動(dòng)方程，采用以下Galerkin離散化方案[19]：

式中，φi(ξ)= 2 sin( )iπξ是簡(jiǎn)支梁的本征函數(shù)(τ)、(τ)、(τ)、(τ)是相應(yīng)的廣義坐標(biāo)。

為了保證結(jié)果的收斂性，在接下來的分析中，截取前六階模態(tài)N=6。將式（27）～（28）代入式（23）～（26），等式兩邊乘以φi(ξ)，從0到1積分，然后利用模態(tài)的正交性，可以將偏微分方程轉(zhuǎn)化為一組常微分方程：

利用四階Runge-Kutta法可以進(jìn)一步數(shù)值求解離散化的常微分方程，得到結(jié)構(gòu)的位移響應(yīng)等物理量。

2 模型驗(yàn)證

?？松梨诠綥ehn[20]在挪威海洋工程研究中心測(cè)試了垂直立管模型在均勻流和剪切流中的渦激振動(dòng)響應(yīng)。Wang 和Xiao[21]利用CFD 方法模擬了立管在均勻流作用下的渦激振動(dòng)響應(yīng)。立管的主要參數(shù)：立管長(zhǎng)度為9.63 m，外徑為20 mm，壁厚為0.45 mm，質(zhì)量比為2.23，長(zhǎng)細(xì)比為481.5，頂張力為817 N，彈性模量E=1.025×1011N/m2，抗彎剛度EI=135 N·m2。不考慮結(jié)構(gòu)阻尼，立管垂直安裝，兩端鉸接。

針對(duì)不同的來流速度，IL 和CF 方向的渦激振動(dòng)運(yùn)動(dòng)軌跡沿立管長(zhǎng)度的分布分別如圖2～3 所示。當(dāng)Ue=0.2 m/s 時(shí)，計(jì)算結(jié)果如圖2 所示：對(duì)于IL 方向，用本理論模型預(yù)測(cè)的運(yùn)動(dòng)軌跡表明，振動(dòng)以第二階模態(tài)為主，這與實(shí)驗(yàn)數(shù)據(jù)和CFD 仿真結(jié)果一致，且IL 方向振幅的數(shù)值計(jì)算結(jié)果也是可以接受的；對(duì)于CF方向，用本理論模型預(yù)測(cè)的運(yùn)動(dòng)軌跡表明，振動(dòng)以第一階模態(tài)為主，這與實(shí)驗(yàn)數(shù)據(jù)和CFD仿真結(jié)果一致，且CF 方向振幅的數(shù)值計(jì)算結(jié)果也是可以接受的。當(dāng)Ue=0.42 m/s時(shí)，計(jì)算結(jié)果如圖3所示，對(duì)于IL 方向，用本理論模型預(yù)測(cè)的運(yùn)動(dòng)軌跡表明，振動(dòng)以第三階模態(tài)為主，這與實(shí)驗(yàn)數(shù)據(jù)和CFD 仿真結(jié)果一致；對(duì)于CF 方向，用本理論模型預(yù)測(cè)的運(yùn)動(dòng)軌跡表明，振動(dòng)以第二階模態(tài)為主，這與實(shí)驗(yàn)數(shù)據(jù)和CFD 仿真結(jié)果一致。由于本文考慮了立管非線性效應(yīng)的影響，致使采用本文方法計(jì)算的立管響應(yīng)包絡(luò)圖幅值小于文獻(xiàn)[20-21]中的實(shí)驗(yàn)和CFD仿真計(jì)算結(jié)果。

圖2 立管運(yùn)動(dòng)軌跡的比較（Ue=0.2 m/s）Fig.2 Comparison of the riser’s response envelopes for Ue=0.2 m/s

圖3 立管運(yùn)動(dòng)軌跡的比較（Ue=0.42 m/s）Fig.3 Comparison of the riser’s response envelopes for Ue=0.42 m/s

3 頂張力對(duì)立管渦激振動(dòng)的影響

本章計(jì)算采用的立管模型及工況主要參數(shù)如表1所示。

表1 立管模型的參數(shù)Tab.1 Parameters of the riser model

3.1 頂張力對(duì)立管固有頻率的影響

在研究有支撐的管道輸送流體的動(dòng)力學(xué)過程中，可以很容易地用線性理論預(yù)測(cè)最低的幾階固有頻率。根據(jù)線性理論，將式（22）和式（23）中的非線性項(xiàng)、外流激勵(lì)項(xiàng)和重力項(xiàng)去掉，整理后得到

圖4 為Ui=0，T0=0～460.42 kN，T=0 時(shí)，頂張力在該區(qū)間變化時(shí)對(duì)立管固有頻率的影響。選取立管前四階固有頻率進(jìn)行研究，結(jié)果表明：相同階次下的立管固有頻率隨頂張力增加逐漸增大，并且隨著立管固有頻率階次增加，立管的固有頻率隨頂張力的變化越快。立管固有頻率隨頂張力的變化規(guī)律可以很好地解釋下文的模態(tài)轉(zhuǎn)換現(xiàn)象。

3.2 頂張力作用下的立管渦激振動(dòng)模態(tài)轉(zhuǎn)換

為了考察平臺(tái)運(yùn)動(dòng)對(duì)立管渦激振動(dòng)的影響，取兩種計(jì)算狀態(tài)：（1）立管上只作用渦升力，海洋平臺(tái)固定不動(dòng)；（2）立管上作用渦升力，同時(shí)海洋平臺(tái)沿軸向（x軸）做簡(jiǎn)諧運(yùn)動(dòng)。

立管在水中的濕重為（mr+mi-ρ0A）gL=160.42 kN，文中Te（t）為立管頂部的張力，為了保證立管底部不出現(xiàn)負(fù)張力，必須滿足Te（t）=T0+Tsin（2πft）≥160.42 kN。第一種計(jì)算狀態(tài)：不考慮平臺(tái)運(yùn)動(dòng)，即T=0，為了滿足Te（t）=T0+Tsin（2πft）≥160.42 kN，取T0=160.42 kN。第二種計(jì)算狀態(tài)：可變張力部分取T=150 kN，為了滿足Te（t）=T0+Tsin（2πft）≥160.42 kN，取T0=310.42 kN。

圖5 給出了當(dāng)預(yù)張力T0=160.42 kN，變張力幅值T=0，外部流體流速Ue=0.2 m/s 時(shí)，立管的位移響應(yīng)。從圖5（c）中可以看出，立管IL 方向的位移響應(yīng)對(duì)應(yīng)第三階模態(tài)，從圖5（d）中可以看出立管CF 方向的位移響應(yīng)對(duì)應(yīng)第二階模態(tài)。

圖4 頂張力對(duì)立管固有頻率的影響Fig.4 Effect of top tension on the natural frequency of the riser

圖5 立管的位移響應(yīng)（T0=160.42 kN）Fig.5 Displacement response of the riser(T0=160.42 kN)

圖6 給出了當(dāng)預(yù)張力T0=460.42 kN，變張力幅值T=0，外部流體流速Ue=0.2 m/s 時(shí)，立管的位移響應(yīng)。從圖6（c）中可以看出立管IL方向的位移響應(yīng)對(duì)應(yīng)第二階模態(tài)，從圖6（d）可以看出立管CF方向的位移響應(yīng)對(duì)應(yīng)第一階模態(tài)。通過對(duì)比可以發(fā)現(xiàn)，不同頂張力作用下立管的渦激振動(dòng)模態(tài)發(fā)生了變化。

圖6 立管的位移響應(yīng)（T0=460.42 kN）Fig.6 Displacement response of the riser(T0=460.42 kN)

時(shí)變頂張力作用下立管的耦合位移響應(yīng)見圖7，此時(shí)，外部流體流速Ue=0.2 m/s，預(yù)張力T0=310.42 kN,變張力幅值T=150 kN，平臺(tái)垂蕩頻率f=0.002 Hz。從圖7（a）可見，t=4850 s附近立管IL方向的位移響應(yīng)V=0.3368～0.4122 m，此刻IL 方向的耦合位移響應(yīng)與圖5 相比基本相同；t=4850～5122 s，立管IL 方向的位移響應(yīng)先減小后增大（V=Vmax-（Vmax+Vmin）/2）；t=5122 s 附近立管IL 方向的耦合位移響應(yīng)V=0.0873～0.1283 m，此刻IL 方向的耦合位移響應(yīng)與圖6 相比略有減小。從圖7（b）可以看出，t=4850 s 附近立管CF 方向的耦合位移響應(yīng)W=-0.172～0.1722 m，此刻CF 方向的耦合位移響應(yīng)與圖5 相比基本相同；t=4850～5122 s，立管CF方向的耦合位移響應(yīng)先減小后增大；t=5122 s附近立管CF方向的耦合位移響應(yīng)W=-0.0893～0.0893 m，此刻CF方向的耦合位移響應(yīng)與圖6相比基本相同。經(jīng)上述對(duì)比發(fā)現(xiàn)，立管IL方向的耦合位移響應(yīng)在第二模態(tài)和第三模態(tài)之間周期性交替出現(xiàn)，立管CF方向的耦合位移響應(yīng)在第一模態(tài)和第二模態(tài)之間周期性交替出現(xiàn)，且隨著張力的變化，立管的耦合位移響應(yīng)幅值周期性地先減小后增大。這種現(xiàn)象的發(fā)生，主要是由于立管軸向張力的變化，使得立管固有頻率隨之變化，從而導(dǎo)致立管發(fā)生耦合位移響應(yīng)的主要模態(tài)交替出現(xiàn)。

圖7 立管的耦合位移響應(yīng)（f=0.002 Hz）Fig.7 Coupled-displacement responses of the riser(f=0.002 Hz)

與圖7相同工況下，僅將平臺(tái)垂蕩頻率變?yōu)閒=0.008 Hz，立管的耦合位移響應(yīng)見圖8。經(jīng)對(duì)比可以發(fā)現(xiàn)，平臺(tái)垂蕩頻率對(duì)立管IL 方向的耦合位移響應(yīng)影響較小，但對(duì)立管CF 方向的耦合位移響應(yīng)影響較大。這可能是由于高頻激勵(lì)作用下張力變化較快，所產(chǎn)生的影響尚未充分發(fā)展，張力就已經(jīng)改變，因此沒有產(chǎn)生較大的耦合振動(dòng)響應(yīng)。

圖8 立管的耦合位移響應(yīng)（f=0.008 Hz）Fig.8 Coupled-displacement responses of the riser(f=0.008 Hz)

4 結(jié) 論

本研究考慮了內(nèi)部流體作用下平臺(tái)頂張力對(duì)立管非線性動(dòng)力學(xué)特性的影響，研究了平臺(tái)垂蕩和立管渦激振動(dòng)的耦合振動(dòng)響應(yīng)，并與立管單獨(dú)渦激振動(dòng)結(jié)果進(jìn)行了比較。通過與已有實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，驗(yàn)證了本文分析方法的合理性。本研究得到了如下結(jié)論:

（1）在考慮海洋平臺(tái)垂蕩和立管自身渦激振動(dòng)的耦合振動(dòng)響應(yīng)時(shí)，由于海洋平臺(tái)的劇烈運(yùn)動(dòng)導(dǎo)致頂張力大幅度周期性變化，從而改變了立管的固有頻率，在一定來流速度下，引發(fā)立管產(chǎn)生不同模態(tài)的渦激振動(dòng)。

（2）隨著頂張力的變化，立管的耦合響應(yīng)在相鄰模態(tài)間轉(zhuǎn)化，立管的耦合位移響應(yīng)幅值周期性地先減小后增大，且隨著平臺(tái)垂蕩頻率的增加，立管CF方向的耦合位移響應(yīng)反而減小。