陳殿鵬，顧洪祿，李福恒，李效民，郭海燕

（中國海洋大學(xué) 工程學(xué)院，山東 青島 266100）

0 引 言

懸鏈線式立管（catenary type risers,CTRs）作為一種靈活性高、無需頂張力補償?shù)牧⒐苄问皆谏詈Ｓ蜌忾_采中得到廣泛應(yīng)用。立管在一定的流場下，管道尾部兩側(cè)會交替出現(xiàn)一對漩渦，漩渦脫落會形成渦激升力致使立管產(chǎn)生周期性振動，該現(xiàn)象稱為渦激振動（vortex-induced vibration,VIV）。當(dāng)渦脫頻率與結(jié)構(gòu)固有頻率在一定范圍內(nèi)相近時，結(jié)構(gòu)會進入鎖定（lock-in）狀態(tài)，此時結(jié)構(gòu)響應(yīng)大幅、迅速地增大，立管更易產(chǎn)生疲勞失效。立管一旦被破壞會造成巨大的經(jīng)濟損失以及嚴重的海洋污染和次生災(zāi)害。因此，對CTR的VIV響應(yīng)進行研究具有重要意義。

目前，國內(nèi)外學(xué)者對垂直立管的VIV已經(jīng)進行了廣泛的數(shù)值研究，這些研究包括了理想狀態(tài)下的均勻流[1-4]、剪切流[5-6]以及接近真實海況的振蕩流[7-8]下垂直立管VIV 響應(yīng)。但是，目前對于CTR（彎曲管）VIV 的數(shù)值模擬還較少。Srinil等[9]通過低階流固耦合模型，研究了均勻流條件下CTR 渦激振動中的多模態(tài)相互作用；饒志標等[10]基于Shear7 軟件分析了剪切流下CTR 振動位移和靜態(tài)位移的關(guān)系；Meng等[11]基于能量守恒原理，通過Van der Pol方程定性地預(yù)測了內(nèi)流速度對CTR-VIV 的影響；Tsukada 等[12]基于獨立性準則預(yù)測了均勻流下CTR 升力幅值響應(yīng)；Zhu 等[13]通過ANSYS MFX 軟件中的求解器研究了指數(shù)剪切流下CTR 的漩渦脫落與結(jié)構(gòu)運動之間的關(guān)系。綜上所述，上述研究都是CTR 在理想穩(wěn)恒流下的相關(guān)研究，對于接近真實海況的振蕩流下CTR的VIV相關(guān)數(shù)值模擬尚未見報道。

本文基于向量式有限元理論，充分考慮了結(jié)構(gòu)的幾何非線性特征以及大變形等非線性行為，應(yīng)用改進的尾流振子模型和獨立性準則，對振蕩流下頂張式立管以及均勻流下CTR 相關(guān)模型試驗進行預(yù)測，進一步對振蕩流下CTR 的VIV 特性從質(zhì)點振動時程曲線、相圖、流固能量傳遞以及響應(yīng)峰值方面進行了分析，有望為CTR工程設(shè)計提供一定的參考。

1 基于向量式有限元的模型建立和求解

初始狀態(tài)假設(shè)CTR 是一根有質(zhì)量但無重力的斜梁，將其離散成N+1個集中質(zhì)點，質(zhì)點之間通過N個無質(zhì)量梁單元連接。如圖1 所示，以立管底部質(zhì)點m1為原點建立坐標系，X軸水平向右，Z軸垂直X軸向上，Y軸垂直X-Z平面向內(nèi)。根據(jù)已知初始邊界條件，通過斜坡加載函數(shù)緩慢施加靜力荷載（重力、浮力）并移動頂端質(zhì)點mN+1至指定位置，最終得到CTR 的靜力平衡位形。初始海流U 沿+X方向，立管初始曲率平面位于X-Z平面內(nèi)。

圖1 向量式有限元中CTR立管模型Fig.1 CTR model in VFIFE

1.1 質(zhì)點控制方程

在向量式有限元中，每個質(zhì)點在三維空間下有三個平移自由度以及三個轉(zhuǎn)動自由度。質(zhì)點在一個時段內(nèi)，位置發(fā)生改變，這個時段叫做途徑單元。任意質(zhì)點J在途徑單元內(nèi)的控制方程如下：

1.2 獨立性原則

King[14]和Ramberg[15]試驗研究發(fā)現(xiàn)，海流作用下傾斜立管的渦激升力以及旋渦脫落頻率與海流垂直分量作用在與其等效的直管時相同，這個假設(shè)被稱為獨立性原則（independence principle）。

該原則首先將海流矢量分解為平行和垂直于立管平面的分量，如圖2 所示。其中垂直于立管平面的海流分量Unh=Uel·sin(θ)，平行于立管平面的海流分量為Uel·cos(θ)。平面內(nèi)分量可進一步分解為垂直于立管單元的分量，大小為Unv=Uel·cos(θ)·sin(α)，以及沿立管單元軸向的分量Uel·cos(θ)·cos(α)，如圖3 所示。對于立管VIV 來說，其升力取決于垂直于該立管單元的有效海流矢量，并且其方向垂直于有效海流分量與該單元軸向組成的平面。垂直于立管的有效海流速度矢量由兩個分量組成：垂直于立管平面的Unh和立管平面內(nèi)的Unv。垂直于立管單元的有效海流速度Un=((Uel·sin(θ))2+(Uel·cos(θ)·sin(α))2)0.5，渦激升力方向如圖4所示。

圖2 平面外海流分解Fig.2 Decomposition of out-plane current

圖3 平面內(nèi)海流分解Fig.3 Decomposition of in-plane current

圖4 升力方向Fig.4 Lift force direction

1.3 尾流振子方程

對于均勻流下CTR的VIV模擬，本文采用Fachinetti等[16]改進的尾流振子方程：

式中：q為尾流振子變量；εy為非線性項中的小參數(shù)，一般取0.3；Ay為流體動力參數(shù)，Ay= 12；ωf=2πStUn/D為漩渦脫落頻率，St為Strouhal數(shù)，D為立管外徑；Y為立管橫向振動位移。

但上述尾流振子方程主要適用于穩(wěn)恒流下立管VIV 的模擬，對于時變振蕩流，許京[17]進一步對其進行了改進，引入無量綱時間tref=tωref和無量綱位移y=Y/D，其中ωref= 2πStUref/D為振蕩流下漩渦脫落參考頻率，Uref表示流速在一個變化周期內(nèi)的平均值，也就是參考值。然后引入振蕩海流參數(shù)得到如下方程：

1.4 控制方程求解

為了克服隱式時間積分帶來的反復(fù)迭代求解以及收斂問題，本文通過中央差分法對n+1 時刻的質(zhì)點位置向量、轉(zhuǎn)角向量以及尾流振子變量q進行求解：

2 模型驗證

2.1 振蕩流下水平直管VIV驗證

本文選取文獻[19]的試驗結(jié)果，將本文模型得到的模擬結(jié)果與其對比，結(jié)果如圖5所示，左列為試驗結(jié)果，右列為模擬結(jié)果。在該試驗中，流速波動方程為

圖5 工況1、2、3和4下z=2 m處VIV響應(yīng)Fig.5 VIV response for Cases 1,2,3 and 4 at z=2 m

式中：Amax、Umax分別為海流最大振蕩幅值和瞬時最大流速，T為振蕩周期。對于振蕩流來說，KC數(shù)是一個重要的描述參數(shù)，一般用KC數(shù)表示結(jié)構(gòu)在水中運動的相對幅度。

立管具體參數(shù)見表1，根據(jù)文獻[8]，本算例Strouhal數(shù)取0.2，同樣選取兩個典型工況組別：KC=178（大KC數(shù)工況）和KC=31（小KC數(shù)工況）進行驗證分析，詳細工況見表2。

表1 振蕩流下柔性立管模型主要參數(shù)Tab.1 Main parameters of flexible riser model under oscillatory flow

表2 對比試驗工況列表Tab.2 List of comparative test conditions

對于時變的振蕩流下的結(jié)構(gòu)響應(yīng)來說，其響應(yīng)模式極為復(fù)雜。但是，如圖5所示，通過對比可以發(fā)現(xiàn)本文模型可以對上述特征現(xiàn)象進行較好的模擬，并且本文模型模擬的振動幅值和周期與試驗結(jié)果吻合度較高。因此，一定程度上說明了本文模型在振蕩流下預(yù)測結(jié)構(gòu)VIV響應(yīng)方面的有效性和準確性。

2.2 均勻流下CTR-VIV驗證

為了進一步驗證本文模型對于預(yù)測具備幾何非線性結(jié)構(gòu)形式（CTR）的VIV 響應(yīng)有效性，選取文獻[20]的試驗結(jié)果與本文結(jié)果進行對比驗證，其中立管模型參數(shù)見表3，本例質(zhì)點總數(shù)取34。CTR 初始位形如圖6 所示，其中海流流速U=0.26 m/s。由圖7 可知，本文模型預(yù)測的模態(tài)數(shù)目和整體振動包絡(luò)峰值與試驗較為一致，并且與文獻[20]的數(shù)值模擬結(jié)果相比，本文模型結(jié)果與試驗結(jié)果更為吻合。經(jīng)過上述驗證，一定程度上證明了本文模型在計算CTR-VIV響應(yīng)方面的有效性。

表3 CTR模型主要參數(shù)Tab.3 Main parameters of CTR model

圖6 均勻流下CTR初始位形Fig.6 Initial configuration of CTR under uniform flow

圖7 文獻[20]與本文得到的VIV包絡(luò)圖對比Fig.7 Comparison between results in Ref.[20]and the VIV envelope obtained in this paper

3 振蕩流下CTR渦激振動分析

本章繼續(xù)選取文獻[20]中CTR 模型進行振蕩流下的VIV 動力響應(yīng)分析，分別取KC=20、40、60 和178 四種KC 數(shù)作為典型的振蕩流工況，其中KC=178 代表的是較為典型的大KC數(shù)工況，KC=20、40 和60 代表小KC數(shù)工況。本章將從質(zhì)點振動時程、相位圖、流固能量傳遞以及結(jié)構(gòu)響應(yīng)峰值四個方面對振蕩流下CTR-VIV進行分析。

3.1 質(zhì)點振動時程響應(yīng)分析

懸垂點是CTR 設(shè)計中關(guān)鍵的控制點，因此本文對CTR 懸垂點（Pt-30）在穩(wěn)定時間30-40 s 的VIV響應(yīng)進行分析，如圖8所示，其中一、三、五行為流速U的時程變化，二、四、六行為懸垂點無量綱位移的時程變化。從圖中可以看出，小KC數(shù)范圍內(nèi)，在相同Umax下，KC數(shù)越大，振蕩海流激勵下質(zhì)點產(chǎn)生的響應(yīng)幅值越小。這是因為隨著KC數(shù)的增大，海流振蕩周期隨之增加，在一個較長的振蕩周期內(nèi)，海流更接近于均勻流的性質(zhì)，這削弱了海流振蕩幅值對立管的激勵作用。同時，該現(xiàn)象也說明海流振蕩周期和幅值同時控制著立管的VIV響應(yīng)。

圖8 不同KC數(shù)和不同最大流速下CTR懸垂點（Pt-30）的VIV響應(yīng)Fig.8 VIV response at the CTR’s overhang point(Pt-30)at different KC numbers and Umax

另一方面，隨著KC數(shù)增大，質(zhì)點振動逐漸呈現(xiàn)典型的振幅調(diào)制現(xiàn)象。值得注意的是，這種振動幅值波動現(xiàn)象與均勻流下立管因模態(tài)轉(zhuǎn)換而發(fā)生的拍頻現(xiàn)象有著本質(zhì)的區(qū)別。在KC=178（大KC數(shù)）工況下，這種振幅調(diào)制現(xiàn)象更為劇烈，質(zhì)點振動產(chǎn)生如文獻[19]所述的建立、鎖定、衰減三個特征階段。海流速度接近為0的時候，質(zhì)點幾乎不振動，出現(xiàn)“間歇性VIV”現(xiàn)象，這可能是由于立管響應(yīng)與時變海流之間存在遲滯現(xiàn)象，文獻[8]同樣指出這種現(xiàn)象在大KC數(shù)下更為明顯。在KC=20 下，Umax的增加對立管VIV 響應(yīng)幅值影響較大，但隨著KC數(shù)增大，Umax對于立管VIV 響應(yīng)無明顯影響，這與均勻流下流速作為關(guān)鍵參數(shù)控制結(jié)構(gòu)響應(yīng)有著顯著的不同。本文也研究了其他位置點，結(jié)果都與懸垂點響應(yīng)特征一致。并且，上述相關(guān)特征與振蕩流下直管的VIV響應(yīng)特征較為相近，這說明CTR的曲率特征對振蕩流下結(jié)構(gòu)響應(yīng)特性無明顯影響。

3.2 相圖分析

圖9是CTR懸垂點Pt-30在KC=20、40、60和178下，流速分別為Umax=0.18 m/s、0.26 m/s和0.34 m/s，時間間隔為30～40 s 下的質(zhì)點振動相圖。其中橫坐標是無量綱橫向振動位移y，縱坐標是無量綱振動速度?y/?t。從圖中可以看出，不同KC數(shù)下質(zhì)點的相圖主形狀有著顯著差異，并且這種差異與Umax關(guān)聯(lián)性較小。KC=20、40、60 和178 下質(zhì)點相圖形狀分別呈“蘋果形”、“心形”、“雙圓環(huán)形”和“玫瑰形”。其中，由于海流的時變性導(dǎo)致立管響應(yīng)呈現(xiàn)一定的隨機性，進而導(dǎo)致相圖呈混沌狀態(tài)，KC數(shù)越大，混沌狀態(tài)越明顯。

圖9 不同KC數(shù)和不同最大流速下CTR懸垂點（Pt-30）的振動相圖Fig.9 Phase portrait of the Pt-30 at different KC numbers and Umax

3.3 流固能量傳遞特性分析

流體與固體之間的能量傳遞不僅可以從VIV 機理角度較好地解釋立管振動特性，而且可以為實際立管設(shè)計提供理論指導(dǎo)。瞬時流固能量傳遞公式為[21]

式中，WL為瞬時能量傳遞量，v是質(zhì)點橫向振動速度，CL為升力系數(shù)。圖10 給出了KC=20、40、60 和178，Umax=0.26 m/s 下CTR 的能量傳遞時空分布，在小KC數(shù)范圍內(nèi)（KC=20、40 和60），能量傳遞峰值分別為0.16、0.11 和0.05，隨著KC數(shù)增大能量傳遞峰值呈遞減趨勢。在KC=178 下能量傳遞峰值并沒有進一步增大，這是因為結(jié)構(gòu)在振蕩流下產(chǎn)生的間歇性VIV 會導(dǎo)致能量在流體阻尼以及自身阻尼的影響下產(chǎn)生折損。并且在KC=60 時，能量傳遞隨時間變化無周期性可言，處于不穩(wěn)定的傳遞狀態(tài)，這可能是由于此時海流激勵能量與結(jié)構(gòu)阻尼處于不平衡轉(zhuǎn)換狀態(tài)導(dǎo)致。當(dāng)KC數(shù)進一步增大時，能量傳遞再次呈現(xiàn)一定的周期性。

圖10 不同KC數(shù)下流固能量傳遞時空分布Fig.10 Space-time distribution of energy transfer at different KC numbers

在四種KC數(shù)下，能量傳遞峰值都出現(xiàn)在立管懸垂段。能量傳遞沿立管全長呈現(xiàn)非對稱性，逐漸向支座兩邊遞減。這是由于模型中考慮了CTR幾何非線性特征，根據(jù)獨立性原則，立管懸垂段處的海流有效分量最大，因此懸垂段受海流激勵產(chǎn)生的能量更多。并且值得注意的是，不同KC數(shù)下，立管都產(chǎn)生了負能量傳遞現(xiàn)象，即能量從立管向海流傳遞，這與我們的直觀判斷相悖。根據(jù)公式（10）可知，這是由于升力系數(shù)與速度產(chǎn)生相位差導(dǎo)致其乘積為負值，這種現(xiàn)象在文獻[22]中剪切流下直管VIV模擬中也出現(xiàn)過，并且這種現(xiàn)象同樣發(fā)生在等效低流速區(qū)域（立管下部），繼而形成結(jié)構(gòu)阻尼區(qū)。這主要是與有效流速剖面的不均勻分布有關(guān)，才導(dǎo)致了立管向流體傳遞能量，即負能量傳遞現(xiàn)象。

3.4 CTR最大振動幅值分析

圖11 是KC=20、40、60 和178 下CTR 整體響應(yīng)最大值ymax隨最大流速Umax的變化趨勢圖。由圖11可知，在小KC數(shù)范圍內(nèi)，在相同Umax下，KC數(shù)越大，立管VIV 響應(yīng)峰值越小。這是因為在較高KC數(shù)下，水動力阻尼對于結(jié)構(gòu)響應(yīng)有明顯抑制作用，隨著KC數(shù)減小，水動力阻尼對結(jié)構(gòu)影響逐漸減弱。而大KC數(shù)下，立管響應(yīng)峰值再次回落，這可能與其渦脫落對數(shù)量以及特征模式有關(guān)。在振蕩流下，由于立管在時變激勵下會產(chǎn)生較為不穩(wěn)定的模態(tài)切換，因此其振動響應(yīng)峰值會隨著流速變化出現(xiàn)“多峰”現(xiàn)象，這種現(xiàn)象在文獻[23]中同樣出現(xiàn)，在KC=20 時，“多峰”現(xiàn)象更為顯著。

圖11 不同KC數(shù)下VIV響應(yīng)峰值隨最大流速變化Fig.11 Peak of VIV response versus Umax at different KC numbers

4 結(jié) 論

本文基于向量式有限元理論對振蕩流下CTR-VIV 響應(yīng)進行了研究，將本文結(jié)果與相應(yīng)的試驗結(jié)果對比，在一定程度上驗證了本模型的有效性。進一步對KC=20、40、60 和178 的振蕩流下的CTRVIV響應(yīng)特性進行了分析，由分析結(jié)果可得到如下結(jié)論：

（1）海流振蕩周期和幅值同時制約著立管振動響應(yīng)幅值；在小KC數(shù)下，立管在振蕩流激勵下產(chǎn)生振幅調(diào)制現(xiàn)象，在大KC數(shù)下質(zhì)點振動響應(yīng)進一步呈現(xiàn)明顯的間歇性渦激振動現(xiàn)象。值得說明的是，CTR的曲率特征對振蕩流下結(jié)構(gòu)響應(yīng)特性無明顯影響。

（2）KC=20、40、60 和178 下質(zhì)點振動相圖主形狀分別呈“蘋果形”、“心形”、“雙圓環(huán)形”和“玫瑰形”，其中Umax對質(zhì)點振動相圖形狀影響不大；振蕩流下質(zhì)點振動具有一定的隨機性，隨著KC數(shù)的增大，相圖混沌狀態(tài)的趨勢更為顯著。

（3）在小KC數(shù)范圍內(nèi)（KC=20、40、60），隨著KC數(shù)增大，流固能量傳遞峰值逐漸減小，這是由于KC數(shù)較大時，結(jié)構(gòu)的不穩(wěn)定振動狀態(tài)導(dǎo)致能量傳遞損耗；由于CTR 的曲率特征導(dǎo)致的不均勻有效海流剖面，使立管在底部低等效流速區(qū)域產(chǎn)生局部阻尼區(qū)，因此產(chǎn)生負能量傳遞現(xiàn)象。

（4）在小KC數(shù)范圍內(nèi)，KC數(shù)越大，由于水動力阻尼的影響，在相同Umax下立管被激勵產(chǎn)生的響應(yīng)峰值ymax越小，隨著KC數(shù)減小，水動力阻尼對其影響衰減；由于大KC數(shù)下結(jié)構(gòu)具備特征的渦脫模式，因此導(dǎo)致其響應(yīng)與小KC數(shù)下有著本質(zhì)的不同。