頂張力立管三維渦激振動(dòng)研究?

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

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

當(dāng)海流流經(jīng)立管時(shí)，在一定條件下會(huì)發(fā)生漩渦脫落的物理現(xiàn)象，漩渦脫落會(huì)對(duì)立管產(chǎn)生周期性的可變力，導(dǎo)致立管發(fā)生垂直于來(lái)流方向的多模態(tài)振動(dòng)，該現(xiàn)象稱為渦激振動(dòng)。渦激振動(dòng)是立管疲勞破壞的主要因素之一，當(dāng)漩渦脫落頻率與立管結(jié)構(gòu)自振頻率接近時(shí)，漩渦的泄放過(guò)程將被結(jié)構(gòu)的振動(dòng)頻率所控制，發(fā)生頻率“鎖定”現(xiàn)象，此時(shí)立管振動(dòng)強(qiáng)烈，振幅增大，加劇疲勞。立管一旦出現(xiàn)安全問(wèn)題，將會(huì)造成巨大的經(jīng)濟(jì)損失和嚴(yán)重的海洋污染和次生災(zāi)害，因此必須加深對(duì)立管渦激振動(dòng)動(dòng)力特性的研究，確保其在服役期間的具有足夠的安全性和可靠性。

針對(duì)立管渦激振動(dòng)的數(shù)值模擬，國(guó)內(nèi)外許多學(xué)者進(jìn)行了大量研究。Park等[1]采用有限元方法研究了細(xì)長(zhǎng)立管的渦激振動(dòng)動(dòng)力響應(yīng)，研究中考慮了浮體運(yùn)動(dòng)對(duì)立管動(dòng)力響應(yīng)的影響，并進(jìn)行了波流激勵(lì)分析和參數(shù)激勵(lì)分析。Holmes等[2]對(duì)網(wǎng)格劃分技術(shù)進(jìn)行改進(jìn)，并提出了一種能夠在網(wǎng)格劃分稀疏的條件下分析立管的方法，該方法能夠模擬任何攻角的位置、斜立立管和其他的幾何形狀，它將三維計(jì)算流體力學(xué)和頂張力立管結(jié)構(gòu)模型相結(jié)合來(lái)預(yù)測(cè)渦激振動(dòng)，克服了二維分析方法的部分缺點(diǎn)。Huang等[3]使用張拉梁模型建立了頂張力立管控制方程，采用了隱式的中央差分法對(duì)立管的運(yùn)動(dòng)控制方程進(jìn)行離散，結(jié)合有限元方法對(duì)長(zhǎng)細(xì)比在482～3 350范圍內(nèi)的3根立管進(jìn)行了數(shù)值模擬，得到的結(jié)果與已發(fā)表的試驗(yàn)數(shù)據(jù)和商用軟件吻合良好，該方法在空間上具有二階精度，在時(shí)間上具有一階精度。

國(guó)內(nèi),徐萬(wàn)海等[4]在研究影響張力腿因素的研究中，基于Euler-Bernouli梁模型分別建立了線性模型和非線性模型，采用有限差分法對(duì)方程進(jìn)行求解，分析了兩種模型在預(yù)測(cè)結(jié)構(gòu)動(dòng)力響應(yīng)方面的差異性，非線性模型的計(jì)算結(jié)果相比線性模型偏小，當(dāng)平臺(tái)運(yùn)動(dòng)幅值增大時(shí)，2種模型計(jì)算結(jié)果差異明顯。唐國(guó)強(qiáng)等[5]采用梁結(jié)構(gòu)模型建立一種經(jīng)驗(yàn)?zāi)Ｐ停瑢?duì)柔性立管渦激振動(dòng)進(jìn)行了分析。唐友剛等[6]考慮張緊器張力變化，結(jié)合尾流振子模型，建立了剪切流場(chǎng)的頂張力立管模型，對(duì)影響立管渦激振動(dòng)特性的各種因素進(jìn)行了參激分析。萬(wàn)德成等[7]基于自主開發(fā)的CFD求解器，采用切片模型，對(duì)恒定張力和變張力的柔性圓柱的渦激振動(dòng)進(jìn)行了仿真模擬，結(jié)果表明變張力條件會(huì)顯著增加圓柱體在振動(dòng)過(guò)程中的位移響應(yīng)。當(dāng)軸向張力以一階固有頻率變化時(shí)，位移增幅約為一倍，當(dāng)軸向張力以二階固有頻率變化時(shí)，位移增幅約為25%。李效民等[8]采用向量式有限元方法建立管動(dòng)力學(xué)模型，編寫了立管動(dòng)力分析程序，對(duì)頂張力立管的動(dòng)力響應(yīng)特性進(jìn)行了分析，通過(guò)與Euler-Bernouli梁模型的數(shù)值模擬結(jié)果進(jìn)行對(duì)比，證明向量式有限元方法可以應(yīng)用于立管的動(dòng)力行為分析。典型的Euler-Bernouli梁模型通常忽略剪切變形，同時(shí)不考慮到軸向位移和軸向力的變化，對(duì)于處于垂直狀態(tài)的頂張力立管，軸向張力扮演著重要角色，其對(duì)立管動(dòng)力特性的影響不可忽略。

目前關(guān)于頂張力立管的研究大部分假設(shè)立管兩端鉸接，沒(méi)有考慮頂端運(yùn)動(dòng)和頂端張力變化的影響。本文考慮軸向變形，采用彈簧模型模擬立管張緊器對(duì)立管作用，基于一種新的有限元方法——向量式有限元方法，以物理模式建立了立管三維非線性振動(dòng)模型，該模型將立管離散為有限數(shù)目的滿足牛頓第二定律的質(zhì)點(diǎn)，以質(zhì)點(diǎn)為控制對(duì)象建立立管控制方程，計(jì)算流程是逐點(diǎn)并逐步循環(huán)的，不需要集成剛度矩陣，避免了矩陣不收斂或出現(xiàn)病態(tài)的情況，通過(guò)顯示中央差分法求解質(zhì)點(diǎn)位移得到結(jié)構(gòu)單元變形，通過(guò)材料力學(xué)理論求解單元內(nèi)力，沒(méi)有作位置變化和變形的分離，對(duì)于求解大變形大變位的立管模型，該方法簡(jiǎn)單可靠同時(shí)求解方便。

1 基于向量式有限元的立管模型

1.1 立管模型控制方程

如圖1所示，假定該頂張力立管長(zhǎng)度為L(zhǎng)，在頂部張力與重力的共同作用下處于垂直位形，作用在立管頂部的張力由剛度為Ks的彈簧來(lái)施加。彈簧用來(lái)替代真實(shí)情況下張緊器的功能，整體坐標(biāo)系的原點(diǎn)設(shè)置在立管的底部，X軸的方向與海流來(lái)流方向平行；Z軸與初始時(shí)刻的立管軸線相重合，正方向指向水面；Y軸與另外兩個(gè)軸的方向相互垂直。假設(shè)立管初始長(zhǎng)度為L(zhǎng)，將立管均分成N段，形成N+1個(gè)質(zhì)點(diǎn)(從下向上依次編號(hào)為1，2，3，…，i，j，…，N，N+1)，假設(shè)任意點(diǎn)的坐標(biāo)為Xi=(xi,yi,zi)T，根據(jù)向量式有限元的點(diǎn)值描述，則立管的初始狀態(tài)可以定義如下：

(1)

式中：i=1,2,3，…，N+1。

質(zhì)點(diǎn)的控制方程如下:

(2)

(3)

本文模型中，初始時(shí)刻只考慮作用在立管上的預(yù)加頂部張力、重力和浮力，后期水動(dòng)力的計(jì)算采用向量式有限元斜坡加載函數(shù)進(jìn)行加載。張力通過(guò)彈簧直接作用到立管頂部的質(zhì)點(diǎn)上，某一時(shí)刻的張力可以表達(dá)為：

T=T0+KsΔtop。

(4)

式中：T0為初始時(shí)刻的預(yù)加頂部張力；Ks為頂部彈簧剛度系數(shù)；Δtop為立管模型頂部質(zhì)點(diǎn)的豎向位移。在初始時(shí)刻，Δtop=0，T=T0。

本文立管模型中用有限數(shù)量的質(zhì)點(diǎn)作為控制變量，通過(guò)質(zhì)點(diǎn)的位置變化和運(yùn)動(dòng)軌跡描述立管結(jié)構(gòu)的運(yùn)動(dòng)狀態(tài)，根據(jù)結(jié)構(gòu)為物理連續(xù)，結(jié)構(gòu)相鄰質(zhì)點(diǎn)之間點(diǎn)的位置使用標(biāo)準(zhǔn)內(nèi)插函數(shù)計(jì)算，結(jié)構(gòu)行為和性質(zhì)對(duì)空間的變化是點(diǎn)值的差，對(duì)時(shí)間的變化是一組常微分方程，從而回避了立管結(jié)構(gòu)控制上使用偏微分方程。利用質(zhì)點(diǎn)描述立管結(jié)構(gòu)的空間位置和立管結(jié)構(gòu)運(yùn)動(dòng)軌跡，在描述結(jié)構(gòu)力學(xué)狀態(tài)方面較為方便，尤其在處理立管這種具有大幾何變化的三維空間柔性結(jié)構(gòu)時(shí)，立管結(jié)構(gòu)的描述不需要使用微分方程控制，因不需要結(jié)構(gòu)分析中的許多簡(jiǎn)化，該模型優(yōu)勢(shì)更為明顯。

1.2 途徑單元

當(dāng)桿件受到外力作用時(shí)，將立管結(jié)構(gòu)上空間點(diǎn)的位置向量看作是一個(gè)有關(guān)時(shí)間的函數(shù)，稱之為途徑或時(shí)間軌跡。如圖2所示，假設(shè)計(jì)算開始及結(jié)束時(shí)間分別為t0和tn，將整個(gè)計(jì)算時(shí)間t0～tn離散為一組時(shí)間點(diǎn)t0，t1，t2，……,tn，其中任一時(shí)段ta～tb被定義為一個(gè)途徑單元。在途徑單元ta≤t≤tb內(nèi)，空間點(diǎn)ta、t、tb的位置向量分別為xa、x、xb，途徑單元解決的問(wèn)題就是空間點(diǎn)從xa轉(zhuǎn)換到xb的問(wèn)題，在轉(zhuǎn)換過(guò)程中，空間點(diǎn)滿足變形和運(yùn)動(dòng)力學(xué)準(zhǔn)則，軌跡則是時(shí)間的連續(xù)函數(shù)。結(jié)構(gòu)運(yùn)動(dòng)滿足的力學(xué)定律為運(yùn)動(dòng)定律和胡克定律，即空間點(diǎn)位移遵循運(yùn)動(dòng)公式，桿件元變形遵循應(yīng)力和應(yīng)變的關(guān)系公式，對(duì)于大變形問(wèn)題，經(jīng)過(guò)簡(jiǎn)化后可以用大變位小變理論來(lái)解決。

1.3 逆向運(yùn)動(dòng)

向量式有限元通過(guò)逆向運(yùn)動(dòng)解決空間桿件元的變形問(wèn)題，如圖3所示，將桿件元從A到B位置先經(jīng)逆向平移，再經(jīng)過(guò)逆向轉(zhuǎn)動(dòng)至途徑單元初始時(shí)刻ta位置，桿件元虛擬形態(tài)與ta時(shí)形態(tài)之間的差異是一個(gè)小變形和小剛體運(yùn)動(dòng)行為，以ta時(shí)形態(tài)作為計(jì)算變形參考，根據(jù)微應(yīng)變和工程應(yīng)力空間計(jì)算桿件元的小變形和小剛體運(yùn)動(dòng)。得到桿件元純變形之后，進(jìn)而可以計(jì)算桿件元節(jié)點(diǎn)內(nèi)力。

1.4 節(jié)點(diǎn)內(nèi)力計(jì)算

根據(jù)材料力學(xué)相關(guān)理論及計(jì)算得到節(jié)點(diǎn)的變形和轉(zhuǎn)角，對(duì)于軸力計(jì)算采用胡克定律，彎矩計(jì)算采用非對(duì)稱截面的撓曲理論，扭矩和剪力計(jì)算采用剪切胡克定律，則節(jié)點(diǎn)內(nèi)力和彎矩計(jì)算如下：

(5)

每一個(gè)桿件元共有12個(gè)力，由于桿件元的質(zhì)量都集中到質(zhì)點(diǎn)，桿件元沒(méi)有質(zhì)量，必須滿足結(jié)構(gòu)的靜力平衡條件，根據(jù)力平衡方程，即可得到其它6個(gè)力分量。

1.5 流體力計(jì)算

立管服役期間的環(huán)境載荷包含眾多載荷，考慮到本文重點(diǎn)研究立管渦激振動(dòng)的動(dòng)力響應(yīng)以及內(nèi)流對(duì)立管渦激振動(dòng)的影響，在本文中主要考慮海流載荷。

考慮海流的流速和流向比較穩(wěn)定而且海流的速度場(chǎng)剖面隨水深的變化是緩慢的，在計(jì)算中近似的認(rèn)定海流載荷為定常載荷。如圖4所示，水平拖曳力FD與水質(zhì)點(diǎn)的運(yùn)動(dòng)方向有關(guān)，由于渦激振動(dòng)會(huì)引起圓柱體的水平振動(dòng)，此時(shí)，海流方向U與圓柱體結(jié)構(gòu)運(yùn)動(dòng)方向V產(chǎn)生相對(duì)速度UR=U-V，渦激振動(dòng)引起的升力垂直于相對(duì)速度方向，順流向的水平拖曳力FD和橫向FL的升力表達(dá)式為：

(6)

(7)

(8)

式中：I=i×i+j×j+k×k為單位張量，i,j,k分別為沿X，Y，Z軸的單位方向矢量。

順流向曳力系數(shù)和橫向升力系數(shù)表達(dá)式為

(9)

(10)

式中：CD0為平均曳力系數(shù)；CflD為彈性支撐運(yùn)動(dòng)圓柱體波動(dòng)曳力系數(shù)；CflD0為固定靜止圓柱體波動(dòng)曳力系數(shù)；CL0為固定靜止圓柱體波動(dòng)升力系數(shù)；ω和q分別為圓柱體順流向和橫向的無(wú)量綱尾流振子系數(shù)。

相對(duì)速度表達(dá)式為：

(11)

(12)

對(duì)于相對(duì)速度絕對(duì)值|UR|，采用泰勒公式，并忽略高階項(xiàng)的方式得到

(13)

可得曳力在X軸方向和Y軸方向分力分別為

(14)

(15)

同理可得升力在X軸和Y軸分力，則立管在X軸是流體力為升力和拖曳力在該方向上的分力之和，同理可得Y軸上的流體力。

1.6 尾流振子模型

本文渦激振動(dòng)模型采用Modified-Rayleigh尾流振子模型，該尾流振子方程為

(16)

(17)

1.7 立管模型數(shù)值求解

本文采用中央差分法求解運(yùn)動(dòng)公式，作為一種常用的顯式積分法，避免了隱式積分法的迭代問(wèn)題和復(fù)雜的收斂問(wèn)題，在求解結(jié)構(gòu)的瞬時(shí)高度非線性問(wèn)題時(shí)，顯式積分法的優(yōu)越性更加明顯。中央差分公式為

(18)

(19)

2 立管渦激振動(dòng)分析

為驗(yàn)證本文渦激振動(dòng)模型的準(zhǔn)確性，采用Chaplin[9]中的試驗(yàn)作為參考，該試驗(yàn)布置簡(jiǎn)圖如圖5所示。立管參數(shù)為管長(zhǎng)13.12 m，管徑0.028 m，長(zhǎng)細(xì)比為467，質(zhì)量線密度為1.85 kg/m，質(zhì)量比為3.0，結(jié)構(gòu)阻尼因子為0.33%，抗彎剛度為29.9 N·m。該試驗(yàn)通過(guò)大氣壓原理巧妙地設(shè)計(jì)了一個(gè)試驗(yàn)裝置，在6.5 m的水槽中成功的進(jìn)行了13.12 m長(zhǎng)的立管渦激振動(dòng)實(shí)驗(yàn)，頂端55%長(zhǎng)度處于靜水中，其余立管被置于運(yùn)動(dòng)水流中，流速最高為1 m/s，雷諾數(shù)范圍為2 800～28 000,頂部彈簧剛度系數(shù)為38.5 kN/m,立管兩端鉸接，頂端施加預(yù)張力。

通過(guò)選取Chaplin等[10]中流速為0.31和0.95 m/s兩種工況進(jìn)行對(duì)比，驗(yàn)證本文模型的有效性。其對(duì)應(yīng)的頂部張力分別為457和1 002 N，數(shù)值模擬中，設(shè)置單元數(shù)為400，時(shí)間步長(zhǎng)為h=2×10-5。

如圖6(a)、7(a)所示，最大順流向無(wú)量綱位移為5.189和10.260，其對(duì)應(yīng)的試驗(yàn)結(jié)果為3.960和11.350，因此計(jì)算結(jié)果同試驗(yàn)結(jié)果較相近。Chaplin的預(yù)報(bào)結(jié)果中僅有2種經(jīng)驗(yàn)?zāi)Ｐ徒o出了順流向的位移預(yù)報(bào)，而且預(yù)報(bào)結(jié)果比試驗(yàn)結(jié)果小很多，相比之下本文模型預(yù)報(bào)結(jié)果在準(zhǔn)確率方面達(dá)到80%以上，在可接受范圍內(nèi)。

圖6(c)、(e)和圖7(c)、(e)皆為立管順流向和橫向的位移時(shí)程圖，立管的橫向振動(dòng)表現(xiàn)為駐波現(xiàn)象，但是順流向發(fā)現(xiàn)行波的存在，此時(shí)，行波從結(jié)構(gòu)振動(dòng)幅度較大的區(qū)域向振動(dòng)幅度較小的區(qū)域傳播，也就是從能量輸入?yún)^(qū)向能量輸出區(qū)傳播，這也符合模型中立管上方水流靜止，下方為運(yùn)動(dòng)水流的情況。隨著流速的增加，行波現(xiàn)象更加明顯(見(jiàn)圖6(b)和圖7(b))，但是這一現(xiàn)象在試驗(yàn)中并未發(fā)現(xiàn)，這可能是試驗(yàn)中的立管具有較小長(zhǎng)細(xì)比而預(yù)加張力較大、干擾因素較多。

0.95 m/s工況的流速雖然要遠(yuǎn)大于0.31 m/s，但是由于對(duì)應(yīng)的頂部預(yù)加張力遠(yuǎn)大于0.31 m/s工況的流速對(duì)應(yīng)的頂部預(yù)加張力，且其順流向位移均方根反而遠(yuǎn)小于0.31 m/s時(shí)的順流向位移均方根。這說(shuō)明相對(duì)流速，頂部張力對(duì)順流向的振幅作用更加明顯，可以通過(guò)增加頂部張力的方式來(lái)減弱橫向渦激振動(dòng)。圖6(b)和圖7(b)顯示順流向位移的均方根要比橫向位移均方根要小很多，只有橫向位移的15%～20%左右，這也是之前的大部分渦激振動(dòng)研究都關(guān)注橫向振動(dòng)的原因，但是考慮到立管渦激振動(dòng)的順流向振動(dòng)頻率接近橫向振動(dòng)頻率的兩倍，對(duì)于立管渦激振動(dòng)的分析，尤其是渦激振動(dòng)疲勞分析，需要重視立管順流向渦激振動(dòng)可能造成的疲勞破壞。

通過(guò)與試驗(yàn)對(duì)比可知，本文建立的模型同時(shí)考慮立管的順流向和橫向渦激振動(dòng)，在預(yù)報(bào)立管渦激振動(dòng)的位移、模態(tài)等方面的準(zhǔn)確程度是較高的，因此，對(duì)立管的渦激振動(dòng)預(yù)報(bào)是有效的。

本文模型采用的時(shí)域分析發(fā)現(xiàn)了行波現(xiàn)象，由于行波效應(yīng)的存在，立管渦激振動(dòng)的動(dòng)力響應(yīng)和導(dǎo)致的疲勞規(guī)律同駐波效應(yīng)有著明顯的區(qū)別，尤其對(duì)于分析行波效應(yīng)起主導(dǎo)作用的大長(zhǎng)細(xì)比立管時(shí)，采用頻域分析方法預(yù)報(bào)立管的疲勞具有一定的局限性，因此在復(fù)雜海洋環(huán)境中需要長(zhǎng)時(shí)間安全服役的立管，在分析其渦激振動(dòng)導(dǎo)致的疲勞時(shí)，采用時(shí)域分析方法更加合理。

3 立管軸向變形和軸向張力分析

向量式有限元采用空間彎曲桿件元來(lái)模擬立管各個(gè)質(zhì)點(diǎn)間的相互作用。質(zhì)點(diǎn)是廣義質(zhì)點(diǎn)，對(duì)于空間問(wèn)題，質(zhì)點(diǎn)具有6個(gè)位移分量即3個(gè)平移分量和3個(gè)轉(zhuǎn)動(dòng)分量，點(diǎn)之間的內(nèi)力有軸力、扭矩和兩組彎矩及剪力，并不像大多數(shù)研究者那樣忽略剪切和扭轉(zhuǎn)變形甚至軸向變形的影響，本文基于向量式有限元建立的模型是一種真正的三維模型，能夠更加真實(shí)模擬立管的行為。

從圖8、9可以看出立管軸向位移和軸向張力是周期性振動(dòng)的，而且頻率相同，其頻率與立管順流向振動(dòng)頻率接近，這說(shuō)明立管的軸向變形主要是由順流向振動(dòng)引起的，這主要是因?yàn)轫樍飨虻钠骄妨σh(yuǎn)大于渦激振動(dòng)升力和波動(dòng)曳力，從而導(dǎo)致立管順流向位移遠(yuǎn)大于橫向位移。盡管立管的軸向位移變化相對(duì)立管整體長(zhǎng)度很小，但是根據(jù)胡克定律，當(dāng)立管材料的勁度系數(shù)很大，較小的軸向位移也會(huì)引起較大的軸向力的變化。立管初始頂端張力為1 002 N，9.84 m處的初始張力為958 N，從圖9中可以看出，立管渦激振動(dòng)達(dá)到穩(wěn)定狀態(tài)時(shí)，9.84 m處的軸向張變化范圍為987.5～1 197 N，比軸向力初始值高出3.08%～24.95%，因此軸向變形引起的軸向張力的增大是不可以忽略的。

從圖10可以看出立管軸向變形最大處接近于立管順流向位移最大處，說(shuō)明軸向變形主要是由于順流向位移引起。從圖11可以看出立管軸向張力從立管頂部到底部是逐漸較小的，立管頂部的張力變化范圍最大，其最大張力超出預(yù)加頂部張力27.3%，立管底部的張力的變化范圍最小，但是其張力仍大于其初始張力值，因?yàn)轫敳繌埩?duì)于頂張力立管的動(dòng)力響應(yīng)有重要影響，因此不能忽略軸向變形引起的軸向張力的變化。本文立管模型的一大優(yōu)點(diǎn)是考慮立管軸向變形和軸向張力變化，是一種真正的三維模型。

4 結(jié)語(yǔ)

本文采用向量式有限元方法，使用彈簧模型模擬立管張緊器對(duì)立管作用，將立管離散為有限數(shù)目的滿足牛頓第二定律的質(zhì)點(diǎn)，以質(zhì)點(diǎn)為控制對(duì)象，給出立管控制方程，建立了立管三維非線性振動(dòng)模型，考慮立管和流體之間相對(duì)速度，采用Modified-Rayleigh尾流振子模型，推導(dǎo)了作用在立管上的流體力計(jì)算公式，并采用顯示中央差分法對(duì)立管模型進(jìn)行數(shù)值求解。

通過(guò)對(duì)立管順流向、橫向和垂向三維耦合渦激振動(dòng)動(dòng)力響應(yīng)進(jìn)行計(jì)算，并將計(jì)算結(jié)果同試驗(yàn)結(jié)果對(duì)比可知，本文模型能較準(zhǔn)確地預(yù)報(bào)立管的渦激振動(dòng)動(dòng)力響應(yīng)。通過(guò)對(duì)立管軸向變形以及軸向張力變化的分析，發(fā)現(xiàn)立管軸向變形會(huì)引起立管軸向張力的明顯增大，影響立管渦激振動(dòng)的動(dòng)力響應(yīng)，在立管的分析中，軸向變形引起的軸向張力的變化是不可忽略的。