頂張力立管三維渦激振動研究?

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

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

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

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

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

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

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

1.1 立管模型控制方程

如圖1所示，假定該頂張力立管長度為L，在頂部張力與重力的共同作用下處于垂直位形，作用在立管頂部的張力由剛度為Ks的彈簧來施加。彈簧用來替代真實(shí)情況下張緊器的功能，整體坐標(biāo)系的原點(diǎn)設(shè)置在立管的底部，X軸的方向與海流來流方向平行；Z軸與初始時(shí)刻的立管軸線相重合，正方向指向水面；Y軸與另外兩個(gè)軸的方向相互垂直。假設(shè)立管初始長度為L，將立管均分成N段，形成N+1個(gè)質(zhì)點(diǎn)(從下向上依次編號為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ù)加頂部張力、重力和浮力，后期水動力的計(jì)算采用向量式有限元斜坡加載函數(shù)進(jìn)行加載。張力通過彈簧直接作用到立管頂部的質(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)作為控制變量，通過質(zhì)點(diǎn)的位置變化和運(yùn)動軌跡描述立管結(jié)構(gòu)的運(yùn)動狀態(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ì)對空間的變化是點(diǎn)值的差，對時(shí)間的變化是一組常微分方程，從而回避了立管結(jié)構(gòu)控制上使用偏微分方程。利用質(zhì)點(diǎn)描述立管結(jié)構(gòu)的空間位置和立管結(jié)構(gòu)運(yùn)動軌跡，在描述結(jié)構(gòu)力學(xué)狀態(tài)方面較為方便，尤其在處理立管這種具有大幾何變化的三維空間柔性結(jié)構(gòu)時(shí)，立管結(jié)構(gòu)的描述不需要使用微分方程控制，因不需要結(jié)構(gòu)分析中的許多簡化，該模型優(yōu)勢更為明顯。

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，途徑單元解決的問題就是空間點(diǎn)從xa轉(zhuǎn)換到xb的問題，在轉(zhuǎn)換過程中，空間點(diǎn)滿足變形和運(yùn)動力學(xué)準(zhǔn)則，軌跡則是時(shí)間的連續(xù)函數(shù)。結(jié)構(gòu)運(yùn)動滿足的力學(xué)定律為運(yùn)動定律和胡克定律，即空間點(diǎn)位移遵循運(yùn)動公式，桿件元變形遵循應(yīng)力和應(yīng)變的關(guān)系公式，對于大變形問題，經(jīng)過簡化后可以用大變位小變理論來解決。

1.3 逆向運(yùn)動

向量式有限元通過逆向運(yùn)動解決空間桿件元的變形問題，如圖3所示，將桿件元從A到B位置先經(jīng)逆向平移，再經(jīng)過逆向轉(zhuǎn)動至途徑單元初始時(shí)刻ta位置，桿件元虛擬形態(tài)與ta時(shí)形態(tài)之間的差異是一個(gè)小變形和小剛體運(yùn)動行為，以ta時(shí)形態(tài)作為計(jì)算變形參考，根據(jù)微應(yīng)變和工程應(yīng)力空間計(jì)算桿件元的小變形和小剛體運(yùn)動。得到桿件元純變形之后，進(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)角，對于軸力計(jì)算采用胡克定律，彎矩計(jì)算采用非對稱截面的撓曲理論，扭矩和剪力計(jì)算采用剪切胡克定律，則節(jié)點(diǎn)內(nèi)力和彎矩計(jì)算如下：

(5)

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

1.5 流體力計(jì)算

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

考慮海流的流速和流向比較穩(wěn)定而且海流的速度場剖面隨水深的變化是緩慢的，在計(jì)算中近似的認(rèn)定海流載荷為定常載荷。如圖4所示，水平拖曳力FD與水質(zhì)點(diǎn)的運(yùn)動方向有關(guān)，由于渦激振動會引起圓柱體的水平振動，此時(shí)，海流方向U與圓柱體結(jié)構(gòu)運(yùn)動方向V產(chǎn)生相對速度UR=U-V，渦激振動引起的升力垂直于相對速度方向，順流向的水平拖曳力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)動圓柱體波動曳力系數(shù)；CflD0為固定靜止圓柱體波動曳力系數(shù)；CL0為固定靜止圓柱體波動升力系數(shù)；ω和q分別為圓柱體順流向和橫向的無量綱尾流振子系數(shù)。

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

(11)

(12)

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

(13)

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

(14)

(15)

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

1.6 尾流振子模型

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

(16)

(17)

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

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

(18)

(19)

2 立管渦激振動分析

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

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

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