基于IVCBC渦方法和動(dòng)剛度矩陣的深海立管的三維計(jì)算模型

龐建華，宗 智，鄭向遠(yuǎn)，周 力

（1.清華大學(xué) 海洋工程技術(shù)中心，廣東 深圳 518021；2.大連理工大學(xué) 工業(yè)裝備與結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024）

0 引 言

渦激振動(dòng)是深海立管普遍存在的一種現(xiàn)象，容易誘發(fā)立管疲勞損傷，降低其使用壽命，影響海洋平臺(tái)的安全，一直受到學(xué)術(shù)界和工程界的廣泛關(guān)注[1]。立管具有較大的長徑尺度比，全尺度的三維數(shù)值模擬對計(jì)算平臺(tái)的要求高，計(jì)算效率低[2]。為此本文基于切片思想降低對計(jì)算平臺(tái)的要求，通過有限體積法將立管離散為多個(gè)有限體積單元，分別計(jì)算每個(gè)單元的流體力，從而達(dá)到了提高效率的目的。

在傳統(tǒng)的立管渦激振動(dòng)中，立管剛度矩陣系數(shù)通過靜態(tài)平衡獲得，所有的系數(shù)在整個(gè)時(shí)域計(jì)算中恒定不變。然而在立管的索模型系統(tǒng)中，切向剛度矩陣是由應(yīng)變能求導(dǎo)獲得，并應(yīng)變能具有時(shí)變性，立管的切向剛度矩陣也應(yīng)具有時(shí)變性，為了提高立管渦激振動(dòng)計(jì)算的準(zhǔn)確性，本文在立管的索模型中引入了動(dòng)剛度矩陣。

本文針對深海立管的獨(dú)特的結(jié)構(gòu)特征，視立管為質(zhì)量集中的多自由度索模型系統(tǒng)[11]，采用有限體積方法將立管離散為多個(gè)有限體單元，通過瞬時(shí)渦量守恒渦方法（IVCBC）[3-10]計(jì)算每個(gè)單元的外載荷，通過應(yīng)變能獲得立管的動(dòng)剛度矩陣，建立計(jì)算立管渦激振動(dòng)的三維數(shù)值計(jì)算新方法。應(yīng)用該方法探索了立管耦合前和耦合后的特性，分析了振型和尾流模型，對立管的設(shè)計(jì)制造有重要的參考價(jià)值和指導(dǎo)意義。

1 三維立管渦激振動(dòng)數(shù)值計(jì)算模型

1.1 立管的自由振動(dòng)分析

結(jié)構(gòu)振動(dòng)時(shí)產(chǎn)生的動(dòng)力響應(yīng)和結(jié)構(gòu)系統(tǒng)的固有特性有著緊密的關(guān)系，所以首先分析立管在自由振動(dòng)時(shí)的動(dòng)力特性將有助于進(jìn)一步分析立管在復(fù)雜動(dòng)力荷載—渦激振動(dòng)作用下的動(dòng)力響應(yīng)。

對于小阻尼結(jié)構(gòu)體系，阻尼力對系統(tǒng)的動(dòng)力特性(振型和頻率)影響很小，因此在進(jìn)行動(dòng)力特性分析時(shí)，通常可以忽略阻尼力的影響。本文采用的是有限體積法，考慮到恢復(fù)力（節(jié)點(diǎn)力）的非線性特性，可以將其寫成增量形式如下：

式中：ω為結(jié)構(gòu)的固有頻率；｛｝A為相應(yīng)的特征向量；本文采用子空間迭代法，根據(jù)（2）式便可容易地求出立管的低階固有頻率和固有振型。

1.2 數(shù)值逐步積分法[12]

立管的位移變化是非線性的，阻尼和剛度與瞬時(shí)位移相互耦合，因此立管的振動(dòng)是一個(gè)非線性過程，其振動(dòng)過程不便采用模態(tài)分析方法，數(shù)值逐步積分法是分析任意非線性系統(tǒng)響應(yīng)方程的唯一普遍適用的方法。其數(shù)值計(jì)算過程如下：

在t0時(shí)刻在慣性力、阻尼力、恢復(fù)力以及外力作用下立管系統(tǒng)達(dá)到平衡狀態(tài)，其平衡狀態(tài)方程為：

式中：fI0、fD0、fS0和Rt0分別為t0時(shí)刻的慣性力、阻尼力、恢復(fù)力和外力。t0+Δt時(shí)刻，系統(tǒng)達(dá)到新的平衡狀態(tài)：

方程式（4）減去方程式（3）可得運(yùn)動(dòng)方程

式中：

在（6）式中 ［C （t ）］和 ［K （t）］為時(shí)間增量期間的阻尼和剛度特性的平均值。 如圖1所示的平均斜率，由于是非線性，只能通過迭代計(jì)算。

把（7）式代入（5）式，得到的增量運(yùn)動(dòng)方程為

圖1 響應(yīng)系數(shù)的表示Fig.1 Presentations of the response coefficients

圖2 基于線性變化的加速度增量Fig.2 Linear variation of acceleration for a time increment

假設(shè)在一個(gè)時(shí)間步長內(nèi)的加速度是線性變化的，且系統(tǒng)的動(dòng)力參數(shù)為常量如圖2所示。由此可得加速度、速度和位移的增量形式：

由（10）式可得到

將（11）式代入（9）式可得

將（11）式和（12）式代入方程式（8）得到位移增量方程

式中：

當(dāng)位移增量Δu（t）確定以后，容易獲得速度增量

該時(shí)間增量最終時(shí)刻的位移和速度向量表示為

該時(shí)間增量最終時(shí)刻的加速度向量可表示為

式中： fD（t+Δ t ）和 fS（t+Δ t）分別代表了t+Δt時(shí)刻的速度和位移條件。

1.3 計(jì)算步驟

A.初始化網(wǎng)格，計(jì)算立管所有控制節(jié)點(diǎn)的原始坐標(biāo)；

B.進(jìn)入時(shí)間步循環(huán)（1）-（6）

（1） 時(shí)間增量 ti+1=ti+Δt；

（2）計(jì)算該時(shí)刻的各控制節(jié)點(diǎn)的外荷載；

（3）計(jì)算立管的平衡位置，即立管在自重和外荷載作用下的各控制節(jié)點(diǎn)的坐標(biāo)；

（4）形成立管的總質(zhì)量陣和總阻尼陣；

（5） 計(jì)算該時(shí)刻立管的切線剛度陣［K（ ti）］，循環(huán)①-⑩

①判斷節(jié)點(diǎn)是典型節(jié)點(diǎn)（內(nèi)部節(jié)點(diǎn)）還是邊界節(jié)點(diǎn)；

②對內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)分別采用適當(dāng)?shù)奶幚矸绞?，組建單元節(jié)點(diǎn)的內(nèi)力矢量；

③對內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)分別采用適當(dāng)?shù)奶幚矸绞?，組建單元切線剛度陣；

④ 組裝總體節(jié)點(diǎn)內(nèi)力向量｛F｝和總體切線剛度陣［K］，并計(jì)算出各節(jié)點(diǎn)的不平衡力｛g｝=｛F｝-｛R｝；

⑤計(jì)算每個(gè)節(jié)點(diǎn)的原始位置向量Δold；

⑥加入邊界條件，進(jìn)行邊界處理；

⑦ 對總體切線剛度矩陣［K］進(jìn)行Crout因式分解，并求解非線性方程（5）得到坐標(biāo)增量矢量｛ΔX ｝；

⑧ 計(jì)算新坐標(biāo)矢量｛X｝new=｛X｝old+｛ΔX ｝；

⑨用方程（7）確定每個(gè)節(jié)點(diǎn)的位置矢量迭代變化值是否收斂。若所有節(jié)點(diǎn)都已經(jīng)滿足β≤βd則表示收斂，否則需要進(jìn)行進(jìn)一步迭代；

⑩若方程不收斂且最大的迭代次數(shù)還沒有達(dá)到，則跳到第①步，否則就結(jié)束計(jì)算；

（9）計(jì)算速度、加速度增量

（10）計(jì)算ti+1時(shí)刻的位移，速度和加速度

C.完成全部時(shí)間步長的迭代，獲得各時(shí)刻的位移、速度、加速度及張力結(jié)果。

上述求解方法和求解步驟均以位移u為基本未知量。

2 方法的驗(yàn)證

2.1 IVCBC渦方法的驗(yàn)證

在本文中采用無量綱時(shí)間步長Δt=aD/U，其中a是修正值，D是立管半徑，U是均勻來流的速度，為了研究時(shí)間步長和渦數(shù)對IVCBC渦方法的計(jì)算精度的影響，修正值設(shè)定a=0.1，0.05，圓柱表面產(chǎn)生新生渦的個(gè)數(shù)分別為32，64，128，由此構(gòu)建出六組工況。圖3展示六種工況在雷諾數(shù)為Re=1×105的壓力分布和Blevins等人[14]的實(shí)驗(yàn)結(jié)果。從這個(gè)圖中能觀察到組合為N=128、Δt=0.05的結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。這表明：渦數(shù)越多，時(shí)間步長越小，計(jì)算精度越高。

為了捕捉到高雷諾數(shù)下較小的渦旋，表面渦的數(shù)目應(yīng)該設(shè)置越大越好，并且時(shí)間步長應(yīng)該越小越好。然而，數(shù)值計(jì)算表明渦數(shù)越大，時(shí)間步長越小，計(jì)算量相應(yīng)地增大。綜合考慮，在下面的數(shù)值模擬中，當(dāng)雷諾數(shù)大于5.5×104我們采用的渦數(shù)N=128，Δt=0.05時(shí)，當(dāng)雷諾數(shù)小于5.5×104，渦數(shù)采用 64，時(shí)間步長為Δt=0.05。

為了驗(yàn)證IVCBC渦方法引起的渦激振動(dòng)，本文采用IVCBC渦方法進(jìn)行渦激振動(dòng)計(jì)算，考慮了二維平面上流體與結(jié)構(gòu)的耦合作用。選用Stappenbelt[15]的實(shí)驗(yàn)來檢驗(yàn)IVCBC渦程序，參數(shù)如表1所示。

圖4給出了不同折合速度下無因次的振動(dòng)頻率和泄渦頻率，其中圖4（a）是Stappenbelt[15]的實(shí)驗(yàn)結(jié)果；圖4（b）是IVCBC渦方法計(jì)算結(jié)果。圖4（a）中，在折合速度5.5-8之間出現(xiàn)振動(dòng)頻率與泄渦頻率相當(dāng)?shù)那闆r，這就是渦激振動(dòng)中有名的鎖定現(xiàn)象（不同于普通共振，鎖定發(fā)生在較大的范圍內(nèi)）。從圖4（a）中明顯看到與試驗(yàn)結(jié)果類似，也有一定振動(dòng)頻率與泄渦頻率相等。IVCBC渦方法計(jì)算出了與實(shí)驗(yàn)一致的鎖定現(xiàn)象，但本文出現(xiàn)在折合速度4-7之間，比試驗(yàn)稍微提前。

圖3 平均壓力系數(shù)（Re=1×105時(shí)間步長為Δt=0.1和Δt=0.05）Fig.3 Mean surface pressure coefficients at Re=1×105with Δt=0.1 and Δt=0.05 by the present method

表1 計(jì)算參數(shù)Tab.1 Computing parameters

圖4 無因次振動(dòng)頻率和泄渦頻率Fig.4 Vibration and vortex shedding frequencies

2.2 三維模型的驗(yàn)證

為了驗(yàn)證計(jì)算深海立管索模型的有效性，本文采用2005年Chaplin[16]大長細(xì)比立管渦激振動(dòng)試驗(yàn)的相關(guān)數(shù)據(jù)進(jìn)行計(jì)算，并將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行了比較。試驗(yàn)是在Delta水池中進(jìn)行，該水池長230 m，寬5 m，拖車速度0.85 m/s。實(shí)驗(yàn)數(shù)據(jù)如下：立管直徑28 mm，質(zhì)量比3.0，立管總長13.12 m，下端45%浸沒在均勻流中，其他部分浸沒在靜水中，彎曲剛度29.2 N/m2，頂端張力939 N，長細(xì)比L/D=466。

圖5中實(shí)線為使用本文所介紹的索模型計(jì)算得到的結(jié)果。從圖中可以看出計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好。圖5中虛線為文獻(xiàn)[16]中使用梁模型計(jì)算結(jié)果。與試驗(yàn)比較，很顯然，其結(jié)果不如索模型的計(jì)算結(jié)果。通過以上比較證明了采用索模型計(jì)算大長細(xì)比立管優(yōu)于梁模型。

圖5 本文計(jì)算結(jié)果及文獻(xiàn)結(jié)果與Chaplin試驗(yàn)結(jié)果的比較Fig.5 Comparisons of results of the present method and Ref.[16]versus experimental results of Chaplin

3 數(shù)值計(jì)算結(jié)果

立管計(jì)算參數(shù)如表2所示。

表2 立管的基本參數(shù)Tab.2 Basic parameters of the riser

3.1 立管耦合前后尾流特征

當(dāng)均勻來流繞過立管時(shí)，在立管兩側(cè)形成周期性的渦脫落，周期性的渦脫落會(huì)對立管產(chǎn)生橫向脈動(dòng)壓力，導(dǎo)致立管在橫向方向上產(chǎn)生變形，立管的振動(dòng)變形改變了尾流中的渦量場，導(dǎo)致作用在立管上的形變發(fā)生變化，因此這種尾流特征和立管的振動(dòng)變形是相互耦合，相互影響的。研究表明：立管的尾流形態(tài)特征主要有三種，分別是P模態(tài)、S模態(tài)和P+S模態(tài)。

圖 6 （a）和（b）給出了當(dāng)雷諾數(shù)為 6×104時(shí)，耦合前后立管在5個(gè)方向不同水平層次的渦量分布圖。從圖中可以看出，未耦合前，沿立管軸向不同水平層處，尾流均表現(xiàn)為2S模態(tài)。耦合后，沿立管軸向不同，在立管橫向振幅較小的水平層處，尾流呈現(xiàn)為2S模態(tài)；在立管橫向振幅較大的水平層處，尾流表現(xiàn)為P+S模態(tài)如圖7所示。同時(shí)看出，在不同切片層處，渦元的速度分布不一樣。這表明，沿立管軸向的尾流具有明顯的三維特征。由圖7進(jìn)一步可知，在耦合前，立管的渦街寬度較小，渦脫落后衰減速度很快，耦合后，尾流渦街寬度變大，尾流分離位置發(fā)生變化，脫落后渦對的衰減速度減小，三維特征更明顯。由此表明，在立管與海流發(fā)生耦合的過程中，立管尾渦的動(dòng)力特性會(huì)受到立管振動(dòng)變形的顯著影響。這不僅導(dǎo)致尾流中渦的強(qiáng)度發(fā)生明顯的變化，還使海流與立管分離的位置也發(fā)生改變。

圖6 耦合前后的渦粒子及其速度分布圖Fig.6 Distribution of vorticity and velocity respect to non-coupling and coupling

圖7 耦合后的2S和P+S尾流模態(tài)Fig.7 ‘2S’ and ‘P+S’ modal wake at coupling

3.2 頻域分析

圖8給出了立管耦合前和耦合后升力系數(shù)的時(shí)歷特性。由圖可知，計(jì)算統(tǒng)計(jì)平均收斂后，升力系數(shù)的幅值趨于一個(gè)常數(shù)，表現(xiàn)為單頻周期性的變化，該頻率為立管泄渦頻率。通過傅里葉分析，可得到渦的脫落頻率為fs=0.54 Hz和斯托哈爾數(shù)為St=0.21。研究表明：當(dāng)雷諾數(shù)在亞臨界區(qū)域內(nèi)時(shí)，一個(gè)剛性圓柱繞流的斯托哈爾數(shù)為0.17～0.21，本文的計(jì)算結(jié)果與實(shí)驗(yàn)基本吻合，但稍微高于實(shí)驗(yàn)值，這是因?yàn)楸痉椒ㄖ胁捎们衅瑴u方法，尾流的三維特性對立管動(dòng)力特性的影響沒有考慮，水平切片層之間的渦的影響被忽略，從而加快了渦的脫落。由此可見，由切片渦的方法得到泄渦頻率時(shí)在立管和繞流耦合前是滿足斯托哈爾數(shù)分布的。由圖可知，耦合后，升力系數(shù)明顯的增大，而且出現(xiàn)多頻率振動(dòng)現(xiàn)象。采用傅里葉分析，Re=60 000時(shí)，泄渦頻率為fs=0.54，0.50 Hz，而斯托哈爾數(shù)為St=0.20和St=0.21，由此可見，耦合后泄渦頻率不再滿足斯托哈爾數(shù)分布，而且立管的振動(dòng)變形對泄渦頻率有顯著影響。不僅僅減緩了泄渦頻率，而且出現(xiàn)多頻泄渦現(xiàn)象。

圖8 耦合前后升力系數(shù)的時(shí)程和功率譜Fig.8 Distributions of lift coefficient and power spectra respect to non-coupling and coupling

3.3 振動(dòng)頻域分析

圖9給出兩節(jié)點(diǎn)（1/4L和3/4L處）的橫向振動(dòng)時(shí)歷變化特性。對圖的結(jié)果進(jìn)行傅里葉轉(zhuǎn)化可得如圖10的結(jié)果。在立管1/4L處，節(jié)點(diǎn)的振動(dòng)位移小，振動(dòng)頻率為fv=0.540 1 Hz和fv=0.447 5 Hz，在立管3/4L處節(jié)點(diǎn)的橫向振動(dòng)位移較大，振動(dòng)頻率為fv=0.540 1 Hz。由此可見在立管不同的位移處，立管節(jié)點(diǎn)橫向振動(dòng)的成分有明顯差異，在橫向振動(dòng)較大的節(jié)點(diǎn)處，其振動(dòng)頻率僅僅包含一個(gè)頻率，在振幅較小的節(jié)點(diǎn)處，其振動(dòng)頻率成分中包含兩個(gè)。與前面的泄渦頻率分析可知：每一個(gè)節(jié)點(diǎn)的振動(dòng)頻率都包含與泄渦頻率一致的主頻，在立管橫向振幅較小的節(jié)點(diǎn)處，不僅包含一個(gè)主頻還包含一個(gè)其他的頻率。

圖9 耦合后立管橫向振動(dòng)位移時(shí)歷Fig.9 Distributions of transverses displacement of riser respect to coupling

圖10 耦合后立管橫向振動(dòng)位移時(shí)歷的功率譜Fig.10 Distributions of power spectra of transverses displacement respect to coupling

3.4 鎖定分析

為了進(jìn)一步分析立管的鎖定現(xiàn)象，在表3中給出了立管從第一階到第十五階模態(tài)的自然頻率，根據(jù)前面的分析可知立管第二階模態(tài)的自然頻率與立管振動(dòng)的泄渦頻率接近，該頻率正好是立管的橫向振動(dòng)頻率。振動(dòng)中還存在其他的頻率，因此在立管振動(dòng)中存在一個(gè)主頻鎖定現(xiàn)象，在不同立管不同位置處的振動(dòng)頻率還存在其他頻率，導(dǎo)致在立管不同位置處出現(xiàn)反差很大的振動(dòng)變形現(xiàn)象。

表3 立管自振頻率Tab.3 The nature frequency of riser

圖11給出立管橫向振動(dòng)振幅達(dá)到最大值時(shí)，橫向振動(dòng)變形的計(jì)算結(jié)果。從圖可知，立管在該雷諾數(shù)下呈現(xiàn)以第三階模態(tài)為主的振動(dòng)特性，導(dǎo)致立管出現(xiàn)復(fù)雜的彎曲變形響應(yīng)，這樣大幅增加了立管疲勞損傷。進(jìn)一步觀察發(fā)現(xiàn)整個(gè)立管的變形出現(xiàn)了非對稱現(xiàn)象，導(dǎo)致非對稱復(fù)雜彎曲變形，這是立管不同深處發(fā)生鎖定現(xiàn)象和非主頻振動(dòng)共同產(chǎn)生的，也是深海立管特有的現(xiàn)象[1]。

由數(shù)值模擬結(jié)果可得，當(dāng)Re=60 000時(shí)，立管無因次最大橫向振幅為A/D=0.21，當(dāng)耦合時(shí)間t=23.2 s時(shí)，最大橫截面應(yīng)力發(fā)生在Z=1/4L處，為 2.75×106N。當(dāng)Re=100 000時(shí)，立管無因次最大橫向振幅為A/D=0.216，當(dāng)耦合時(shí)間t=23.9 s時(shí)，最大橫截面應(yīng)力發(fā)生在立管底部，為 1.134×107N。由此可知，隨著雷諾數(shù)的增加，渦激振動(dòng)強(qiáng)度增大，立管振幅和橫截面應(yīng)力隨著增大，渦激振動(dòng)對立管的危害增強(qiáng)。同時(shí)，不同的雷諾數(shù)在不同時(shí)刻，立管最大橫截面應(yīng)力所在位置也會(huì)發(fā)生相應(yīng)的變化。因此，在大長徑比深海立管的設(shè)計(jì)與應(yīng)用中，應(yīng)該對相關(guān)參數(shù)進(jìn)行評估，合理地確定立管危險(xiǎn)段的范圍，便于結(jié)構(gòu)加強(qiáng)處理。

圖11 振幅達(dá)到最大值時(shí)，立管的橫向振動(dòng)變形Fig.11 The deformation of the riser at maximum of the transverse displacement

3.5 立管的振動(dòng)模型

圖12展示了海流為U=1.4 m/s時(shí)，在立管0～1/3L處兩自由度的振動(dòng)圖，其中X軸為橫向振動(dòng)位移，Y軸為流向振動(dòng)位移。從圖中看出，立管的振動(dòng)，在幅值較大的節(jié)點(diǎn)處，以“∞”運(yùn)動(dòng)，在幅值較小的節(jié)點(diǎn)處，橫向位移較小，節(jié)點(diǎn)振動(dòng)表現(xiàn)為 “一”。

圖12 節(jié)點(diǎn)橫向和流向運(yùn)動(dòng)位移Fig.12 Stream-wise displacements and transverse displacements of nodes

4 結(jié) 論

針對深海立管長細(xì)比非常大的結(jié)構(gòu)特征，本文視立管為質(zhì)量集中的多自由度索模型系統(tǒng)。通過有限體積法將該模型系統(tǒng)離散為多個(gè)有限體單元，首次提出基于應(yīng)變能計(jì)算立管動(dòng)態(tài)剛度矩陣的算法，并采用IVCBC渦方法計(jì)算有限體的外載荷，構(gòu)建了一種三維數(shù)值研究深海立管渦激振動(dòng)的新方法，同時(shí)探索了深海立管耦合前后的振動(dòng)和尾流特征。結(jié)果表明，立管渦激振動(dòng)過程中，渦泄頻率不再滿足Strouhal規(guī)律，出現(xiàn)渦泄頻率顯著增大和多頻渦泄等現(xiàn)象，從而導(dǎo)致立管的渦激振動(dòng)出現(xiàn)多頻“鎖定”現(xiàn)象；立管振動(dòng)響應(yīng)頻率除了與渦泄頻率一致的主頻振動(dòng)外，還包含其它成分的次頻振動(dòng)；在立管不同橫截面處出現(xiàn)不同頻率的振動(dòng)響應(yīng)現(xiàn)象，導(dǎo)致立管出現(xiàn)多種高階模態(tài)振動(dòng)共存的非對稱復(fù)雜彎曲變形；深海立管渦激振動(dòng)中出現(xiàn)的多頻、高階模態(tài)和非對稱振動(dòng)響應(yīng)等現(xiàn)象，導(dǎo)致立管最大橫截面應(yīng)力所在位置發(fā)生變化。因此，在深海立管的設(shè)計(jì)與應(yīng)用中，必須合理地確定其危險(xiǎn)段的可能范圍，并采取措施對多個(gè)危險(xiǎn)區(qū)進(jìn)行結(jié)構(gòu)加強(qiáng)處理。