龐建華,宗 智,鄭向遠(yuǎn),周 力
(1.清華大學(xué) 海洋工程技術(shù)中心,廣東 深圳 518021;2.大連理工大學(xué) 工業(yè)裝備與結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連 116024)
渦激振動(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)意義。
結(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)式便可容易地求出立管的低階固有頻率和固有振型。
立管的位移變化是非線性的,阻尼和剛度與瞬時(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í)刻的速度和位移條件。
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為基本未知量。
在本文中采用無量綱時(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
為了驗(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
立管計(jì)算參數(shù)如表2所示。
表2 立管的基本參數(shù)Tab.2 Basic parameters of the riser
當(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
圖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
圖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
為了進(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
圖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
針對深海立管長細(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)處理。