基于IVCBC渦方法和動剛度矩陣的深海立管的三維計算模型

龐建華，宗 智，鄭向遠，周 力

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

0 引 言

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

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

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

1 三維立管渦激振動數(shù)值計算模型

1.1 立管的自由振動分析

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

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

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

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

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

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

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

方程式（4）減去方程式（3）可得運動方程

式中：

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

把（7）式代入（5）式，得到的增量運動方程為

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

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

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

由（10）式可得到

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

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

式中：

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

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

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

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

1.3 計算步驟

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

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

（1） 時間增量 ti+1=ti+Δt；

（2）計算該時刻的各控制節(jié)點的外荷載；

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

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

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

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

②對內(nèi)部節(jié)點和邊界節(jié)點分別采用適當?shù)奶幚矸绞剑M建單元節(jié)點的內(nèi)力矢量；

③對內(nèi)部節(jié)點和邊界節(jié)點分別采用適當?shù)奶幚矸绞剑M建單元切線剛度陣；

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

⑤計算每個節(jié)點的原始位置向量Δold；

⑥加入邊界條件，進行邊界處理；

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

⑧ 計算新坐標矢量｛X｝new=｛X｝old+｛ΔX ｝；

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

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

（9）計算速度、加速度增量

（10）計算ti+1時刻的位移，速度和加速度

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

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

2 方法的驗證

2.1 IVCBC渦方法的驗證

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

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

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

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

圖3 平均壓力系數(shù)（Re=1×105時間步長為Δ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 計算參數(shù)Tab.1 Computing parameters

圖4 無因次振動頻率和泄渦頻率Fig.4 Vibration and vortex shedding frequencies

2.2 三維模型的驗證

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

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

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

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

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

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

3.1 立管耦合前后尾流特征

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

圖 6 （a）和（b）給出了當雷諾數(shù)為 6×104時，耦合前后立管在5個方向不同水平層次的渦量分布圖。從圖中可以看出，未耦合前，沿立管軸向不同水平層處，尾流均表現(xiàn)為2S模態(tài)。耦合后，沿立管軸向不同，在立管橫向振幅較小的水平層處，尾流呈現(xiàn)為2S模態(tài)；在立管橫向振幅較大的水平層處，尾流表現(xiàn)為P+S模態(tài)如圖7所示。同時看出，在不同切片層處，渦元的速度分布不一樣。這表明，沿立管軸向的尾流具有明顯的三維特征。由圖7進一步可知，在耦合前，立管的渦街寬度較小，渦脫落后衰減速度很快，耦合后，尾流渦街寬度變大，尾流分離位置發(fā)生變化，脫落后渦對的衰減速度減小，三維特征更明顯。由此表明，在立管與海流發(fā)生耦合的過程中，立管尾渦的動力特性會受到立管振動變形的顯著影響。這不僅導致尾流中渦的強度發(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ù)的時歷特性。由圖可知，計算統(tǒng)計平均收斂后，升力系數(shù)的幅值趨于一個常數(shù)，表現(xiàn)為單頻周期性的變化，該頻率為立管泄渦頻率。通過傅里葉分析，可得到渦的脫落頻率為fs=0.54 Hz和斯托哈爾數(shù)為St=0.21。研究表明：當雷諾數(shù)在亞臨界區(qū)域內(nèi)時，一個剛性圓柱繞流的斯托哈爾數(shù)為0.17～0.21，本文的計算結(jié)果與實驗基本吻合，但稍微高于實驗值，這是因為本方法中采用切片渦方法，尾流的三維特性對立管動力特性的影響沒有考慮，水平切片層之間的渦的影響被忽略，從而加快了渦的脫落。由此可見，由切片渦的方法得到泄渦頻率時在立管和繞流耦合前是滿足斯托哈爾數(shù)分布的。由圖可知，耦合后，升力系數(shù)明顯的增大，而且出現(xiàn)多頻率振動現(xiàn)象。采用傅里葉分析，Re=60 000時，泄渦頻率為fs=0.54，0.50 Hz，而斯托哈爾數(shù)為St=0.20和St=0.21，由此可見，耦合后泄渦頻率不再滿足斯托哈爾數(shù)分布，而且立管的振動變形對泄渦頻率有顯著影響。不僅僅減緩了泄渦頻率，而且出現(xiàn)多頻泄渦現(xiàn)象。

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

3.3 振動頻域分析

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

圖9 耦合后立管橫向振動位移時歷Fig.9 Distributions of transverses displacement of riser respect to coupling

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

3.4 鎖定分析

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

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

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

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

圖11 振幅達到最大值時，立管的橫向振動變形Fig.11 The deformation of the riser at maximum of the transverse displacement

3.5 立管的振動模型

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

圖12 節(jié)點橫向和流向運動位移Fig.12 Stream-wise displacements and transverse displacements of nodes

4 結(jié) 論

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