武 磊，趙偉文，萬德成

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點(diǎn)實(shí)驗(yàn)室 船海計(jì)算水動(dòng)力學(xué)研究中心(CMHL)，上海 200240)

深海柔性立管因其長細(xì)比大，易產(chǎn)生高階模態(tài)渦激振動(dòng)，進(jìn)而造成立管結(jié)構(gòu)疲勞損傷乃至破壞，成為海洋工程領(lǐng)域長期以來研究的熱點(diǎn)問題之一。區(qū)別于剛性圓柱渦激振動(dòng)，柔性立管在海洋來流的作用下更容易激發(fā)出高階模態(tài)的渦激振動(dòng)，并表現(xiàn)出多模態(tài)振動(dòng)的特性。關(guān)于柔性單立管渦激振動(dòng)問題的試驗(yàn)研究以及數(shù)值研究已廣泛開展，Wu等[1]對細(xì)長柔性單立管渦激振動(dòng)問題的研究進(jìn)展進(jìn)行了總結(jié)。隨著深海中海洋平臺立管數(shù)量的逐漸增加，多立管渦激振動(dòng)問題越發(fā)常見。相較于單立管渦激振動(dòng)情形，立管間的尾流干擾以及立管間距的影響使得多立管渦激振動(dòng)響應(yīng)更為復(fù)雜。

由于多立管渦激振動(dòng)研究的開展難度較大，目前關(guān)于多立管渦激振動(dòng)的研究主要集中在雙立管情形，且多為試驗(yàn)研究。柔性雙立管渦激振動(dòng)試驗(yàn)研究主要集中在串列、并列以及錯(cuò)列雙立管布置形式，主要關(guān)注點(diǎn)包括不同立管間距下雙立管的渦激振動(dòng)響應(yīng)特性[2-6]以及上游立管固定情況下下游立管的振動(dòng)響應(yīng)特性[7-8]。Assi[8]通過水槽試驗(yàn)發(fā)現(xiàn)，在上游圓柱固定且雙圓柱順流向間距固定為2倍圓柱直徑時(shí)，在橫流向間距小于2倍圓柱直徑的工況中，下游圓柱振動(dòng)表現(xiàn)為典型的流致振動(dòng)(wake-induced vibration，簡稱WIV)，在橫流向間距為3倍圓柱直徑的工況中則表現(xiàn)為明顯的渦激振動(dòng)(vortex-induced vibration,簡稱VIV)。

目前關(guān)于柔性雙立管渦激振動(dòng)的數(shù)值研究較為有限。González等[9]通過求解雷諾平均納維爾-斯托克斯(Reynolds-averaged Navier Stokes，簡稱RANS)方程和結(jié)構(gòu)動(dòng)力學(xué)控制方程，對Huera Huarte等[7]的試驗(yàn)進(jìn)行了數(shù)值驗(yàn)證計(jì)算。計(jì)算結(jié)果表明，流固耦合計(jì)算的數(shù)值穩(wěn)定性受單個(gè)時(shí)間步內(nèi)流場與結(jié)構(gòu)場耦合次數(shù)的影響較大，而流場網(wǎng)格質(zhì)量、內(nèi)迭代次數(shù)以及時(shí)間步的影響不大。Lin等[10]參考Huera Huarte等[3]的試驗(yàn)，采用離散渦方法(Discrete vortex method，簡稱DVM)模擬階梯流中柔性雙立管的渦激振動(dòng)，數(shù)值結(jié)果與試驗(yàn)結(jié)果吻合良好。Chen等[11]采用重疊網(wǎng)格方法對Huera Huarte等[3]的試驗(yàn)工況進(jìn)行了數(shù)值模擬，研究發(fā)現(xiàn)立管間距為3倍立管直徑時(shí)，下游立管表面出現(xiàn)了明顯的尾流重附著現(xiàn)象。Wang等[12]通過數(shù)值模擬研究了雷諾數(shù)為500時(shí)立管間距對串列柔性雙立管渦激振動(dòng)響應(yīng)的影響，研究結(jié)果表明當(dāng)立管間距足夠大使得上游立管的尾渦可以從上游立管分離時(shí)，上游立管表現(xiàn)為典型的VIV，而下游立管表現(xiàn)為WIV。

從以上分析可知，目前關(guān)于柔性雙立管渦激振動(dòng)的試驗(yàn)研究主要關(guān)注不同布置形式下立管間距對雙立管渦激振動(dòng)響應(yīng)的影響，未關(guān)注來流剖面沿立管展向分布的影響，而數(shù)值研究主要針對模型試驗(yàn)工況進(jìn)行驗(yàn)證計(jì)算且相關(guān)研究工作比較有限。

因此，采用viv-FOAM-SJTU求解器，研究階梯流中立管浸沒長度對串列柔性雙立管渦激振動(dòng)響應(yīng)的影響，求解器的有效性已被Duan等[13]、Fu等[14]和Wu等[15]驗(yàn)證。第一部分介紹采用的數(shù)值方法；第二部分介紹數(shù)值模擬的計(jì)算模型及工況設(shè)置；第三部分給出數(shù)值結(jié)果并進(jìn)行必要的分析，該部分首先參考Huera Huarte等[3]的試驗(yàn)進(jìn)行驗(yàn)證計(jì)算，然后改變立管浸沒長度并給出各工況下的數(shù)值計(jì)算結(jié)果；最后對全文的數(shù)值計(jì)算結(jié)果進(jìn)行總結(jié)。

1 數(shù)值方法

1.1 切片理論

圖1 切片模型

考慮到深海柔性立管的軸向尺度較大，若對其渦激振動(dòng)進(jìn)行全三維數(shù)值模擬將消耗大量計(jì)算資源，因此采用切片理論對流場進(jìn)行簡化。切片理論即使用二維流場切片沿立管展向進(jìn)行均勻劃分，并認(rèn)為將切片所受水動(dòng)力進(jìn)行插值可得到模型整體受力，切片理論的有效性已被Herfjord等[16]證明。如圖1所示為立管的切片模型。

1.2 流場控制方程

針對每個(gè)切片，通過求解不可壓縮RANS方程得到切片所受的流體力。RANS方程表述如下：

(1)

(2)

1.3 結(jié)構(gòu)動(dòng)力學(xué)控制方程

結(jié)構(gòu)動(dòng)力學(xué)計(jì)算部分，將立管視為兩端簡支的小位移歐拉-伯努利彎曲梁模型，忽略頂張力隨時(shí)間的變化。來流中均質(zhì)立管在順流向及橫流向的受力平衡方程可表示為式(3)、(4)。

(3)

(4)

式中：EI為立管彎曲剛度；T(z)=T-ωs(L-z)為立管軸向張力，其中T為立管頂張力，ωs為立管在水中的單位長度重力，已考慮浮力影響，L為立管長度；m、c、k分別為單位長度立管的質(zhì)量、阻尼以及剛度。

考慮到立管兩端邊界條件為簡支，即位移和彎矩為0，通過有限元方法將式(3)、(4)離散，可得到立管的結(jié)構(gòu)運(yùn)動(dòng)控制方程，即式(5)、(6)中的二階常微分方程，詳細(xì)推導(dǎo)過程可參考文獻(xiàn)[17]，其求解采用紐馬克-貝塔(Newmark-β)[18]法。

(5)

(6)

式中：M、C、K分別為質(zhì)量矩陣、阻尼矩陣和剛度矩陣；x、y分別為立管結(jié)構(gòu)順流向和橫流向的位移響應(yīng)，變量上方的圓點(diǎn)代表對時(shí)間的導(dǎo)數(shù)；Fx、Fy分別為立管順流向和橫流向所受的水動(dòng)力載荷。阻尼矩陣C基于瑞利阻尼模型獲取，表述如下：

C=αM+βK

(7)

圖2 流固耦合求解流程

(8)

流場部分和結(jié)構(gòu)動(dòng)力學(xué)部分的耦合計(jì)算通過流固耦合模塊完成，如圖2所示。在每個(gè)時(shí)間步內(nèi)，首先將流場計(jì)算模塊求解得到的各切片所受流體力，插值映射到各結(jié)構(gòu)節(jié)點(diǎn)上；然后，依此計(jì)算結(jié)構(gòu)位移響應(yīng)；進(jìn)而將結(jié)構(gòu)節(jié)點(diǎn)位移映射到流場，以完成流場網(wǎng)格的更新；最后完成時(shí)間步的向前推進(jìn)。其中，流場網(wǎng)格的變形操作采用OpenFOAM自帶的Laplace方法完成。

2 問題描述

文中的數(shù)值計(jì)算模型參考Huera Huarte等[3]開展的階梯流中串列雙立管渦激振動(dòng)系列試驗(yàn)。由于模型試驗(yàn)只關(guān)注下游立管的振動(dòng)響應(yīng)，只在下游立管上安裝測量儀器，因此雙立管的結(jié)構(gòu)參數(shù)略有不同，詳細(xì)參數(shù)如表1所示(上標(biāo)a代表上游立管)。

模型試驗(yàn)中，雙立管串列布置，立管間距為3D。雙立管處于階梯流中，即立管下部40%長度處于均勻流中，流速為0.55 m/s，立管上部60%長度處于空氣中。本文的數(shù)值模擬中，首先對試驗(yàn)工況進(jìn)行驗(yàn)證計(jì)算；然后將流速設(shè)為0.3 m/s，并改變立管浸沒長度，計(jì)算浸沒長度在40%～70%(間隔15%)立管長度下的雙立管振動(dòng)響應(yīng)。

表1 立管模型的主要結(jié)構(gòu)參數(shù)

由于立管上部處于空氣中的部分長度所受的升阻力較小，參考Chen等[11]研究忽略立管處于空氣中部分長度所受的升阻力載荷。因此雙立管處于空氣中的長度部分不設(shè)置切片，網(wǎng)格總量得以減小。圖3為階梯流中串列雙立管的切片布置示意，切片沿立管軸向均勻布置，其數(shù)量隨立管浸沒長度改變。所有切片的計(jì)算域以及網(wǎng)格劃分一致，如圖4所示，網(wǎng)格y+值約為3。流體沿x正向流動(dòng)，流動(dòng)入口距上游立管15D，流動(dòng)出口距下游立管30D，上下邊界分別距立管中心15D?？拷⒐苤行膮^(qū)域的網(wǎng)格劃分較密，靠近計(jì)算域邊界的網(wǎng)格較為稀疏。流體控制方程的邊界條件設(shè)置為：流動(dòng)入口的速度邊界條件與來流一致，流動(dòng)出口設(shè)置為相對值為0的壓力出口，雙立管表面設(shè)置為無滑移邊界，垂直于切片上下側(cè)面以及立管展向邊界方向的速度設(shè)為0。

結(jié)構(gòu)動(dòng)力學(xué)計(jì)算部分，將立管模型均勻劃為80個(gè)結(jié)構(gòu)單元，總共81個(gè)結(jié)構(gòu)節(jié)點(diǎn)。其中，處于均勻流場中的結(jié)構(gòu)單元數(shù)隨浸沒長度均勻變化。

圖3 串列雙立管的切片布置

圖4 切片網(wǎng)格劃分

3 計(jì)算結(jié)果及分析

3.1 試驗(yàn)工況驗(yàn)證計(jì)算

為確保數(shù)值計(jì)算結(jié)果的有效性，首先參考試驗(yàn)工況開展了驗(yàn)證計(jì)算。驗(yàn)證計(jì)算工況選取立管間距S=3D，頂張力T=110 N，流速U=0.55 m/s。針對試驗(yàn)工況的驗(yàn)證計(jì)算共選取3套網(wǎng)格，網(wǎng)格1、網(wǎng)格2以及網(wǎng)格3中單個(gè)切片的網(wǎng)格量分別為26 000、45 000以及80 000。圖5給出了3套網(wǎng)格中的上下游立管橫流向以及順流向的振動(dòng)位移均方根曲線，可以看出網(wǎng)格2與網(wǎng)格3的計(jì)算結(jié)果基本吻合，而網(wǎng)格1由于過分稀疏，計(jì)算結(jié)果稍差。整體而言，求解器的網(wǎng)格收斂性良好，因此文中選取網(wǎng)格2作為計(jì)算網(wǎng)格。

圖5 立管橫流向及順流向位移均方根

圖6為下游立管中間節(jié)點(diǎn)的橫流向振動(dòng)位移時(shí)歷曲線，可以看出中間節(jié)點(diǎn)的橫流向振幅最大值可達(dá)0.8D。圖7給出了下游立管的瞬時(shí)振動(dòng)輪廓疊加圖。其中圖7(a)、7(b)、7(c)分別表示模型試驗(yàn)中下游立管相對平衡位置的橫流向振動(dòng)輪廓、順流向振動(dòng)輪廓以及相對于平衡位置的順流向振動(dòng)輪廓；圖7(d)、7(e)、7(f)分別是與試驗(yàn)結(jié)果相對應(yīng)的數(shù)值模擬結(jié)果。

圖6 下游立管中間節(jié)點(diǎn)橫流向振動(dòng)時(shí)歷曲線

圖7 下游立管的瞬時(shí)振動(dòng)輪廓疊加圖((a)、(b)、(c)為試驗(yàn)結(jié)果，(d)、(e)、(f)為數(shù)值結(jié)果)

可以看出，下游立管橫流向及順流向的主振模態(tài)均為一階，數(shù)值計(jì)算結(jié)果與模型試驗(yàn)結(jié)果基本一致。模型試驗(yàn)中，下游立管橫流向振幅范圍為-0.75D～0.68D，順流向振幅范圍為-0.03D～0.18D；與試驗(yàn)結(jié)果相對應(yīng)的數(shù)值模擬結(jié)果分別為，橫流向振幅范圍-0.74D～0.68D，順流向振幅范圍-0.04D～0.17D，數(shù)值誤差在可接受范圍內(nèi)。

圖8給出了下游立管中間節(jié)點(diǎn)的位移功率譜密度(Power spectral density，簡稱PSD)，其中圖8(a)、8(b)分別為模型試驗(yàn)中橫流向及順流向的位移功率譜密度，圖8(c)、8(d)分別為與試驗(yàn)結(jié)果相對應(yīng)的數(shù)值結(jié)果?？梢钥闯?，模型試驗(yàn)中，下游立管的橫流向振動(dòng)為單一頻率的一階振動(dòng)模態(tài)，頻率約為6 Hz；順流向振動(dòng)的主振模態(tài)仍為一階，主振頻率約為6 Hz，同時(shí)存在較為微弱的二階模態(tài)成分，頻率約為12.5 Hz，可以看出，下游立管順流向表現(xiàn)出了多模態(tài)振動(dòng)特性。數(shù)值模擬中，下游立管的橫流向振動(dòng)響應(yīng)與試驗(yàn)結(jié)果一致，為單一頻率的一階振動(dòng)模態(tài)，主振頻率為6 Hz；順流向主振頻率仍為6 Hz，并且位移功率譜中同樣發(fā)現(xiàn)了較微弱的二階模態(tài)成分，而頻率約為14 Hz，與試驗(yàn)結(jié)果有一定的誤差。引起順流向二階振動(dòng)頻率誤差的原因可能是，數(shù)值模擬中采用二維切片模擬流場，忽略了尾渦的三維效應(yīng)，并且在將水動(dòng)力載荷映射到結(jié)構(gòu)單元時(shí)對每個(gè)切片所代表的立管長度范圍采用均勻插值，難以精確模擬實(shí)際水動(dòng)力載荷的分布，因此造成數(shù)值計(jì)算所得的振動(dòng)頻率與試驗(yàn)結(jié)果存在誤差。但總體而言，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果基本一致，求解器的有效性較為可信。

圖8 下游立管中間節(jié)點(diǎn)的位移功率譜密度

3.2 浸沒長度對立管渦激振動(dòng)的影響

基于對試驗(yàn)工況進(jìn)行的驗(yàn)證計(jì)算，文中進(jìn)一步研究了浸沒長度對立管渦激振動(dòng)的影響。選取的計(jì)算工況是來流速度為0.3 m/s，串列雙立管的浸沒長度分別為0.40L、0.55L、0.70L，分別記為工況1、工況2、工況3，其余的計(jì)算設(shè)置與3.1節(jié)中的驗(yàn)證計(jì)算工況一致。本節(jié)將分別介紹不同浸沒長度下上游立管及下游立管的渦激振動(dòng)響應(yīng)特性。

3.2.1 上游立管振動(dòng)響應(yīng)

圖9(a)和圖9(b)分別給出了上游立管橫流向及順流向沿立管展向的無量綱化位移均方根?？梢钥闯錾嫌瘟⒐艿臋M流向振動(dòng)在各工況下均為一階，而順流向振動(dòng)僅在工況1中為一階，在工況2和工況3中主振模態(tài)增大為二階。隨著浸沒長度的增加，立管整體所受的水動(dòng)力載荷有所增大。由于數(shù)據(jù)取樣時(shí)間的差異，盡管工況2和工況3的最大位移均方根基本相同，但立管橫流向的位移均方根整體上呈現(xiàn)出增大的趨勢。圖9(b)展示的順流向位移均方根結(jié)果則有所不同。在工況1中，立管順流向主振模態(tài)為一階，頻率較低，立管的順流向振動(dòng)得以充分發(fā)展，其最大位移均方根大于后兩種工況；而在工況2和工況3中，立管順流向的主振模態(tài)變?yōu)槎A，振動(dòng)頻率較高，立管的順流向周期性振動(dòng)來不及充分發(fā)展便已進(jìn)入下一周期，因此最大位移均方根較工況1有明顯的減小。

圖9 上游立管無量綱化位移均方根

上游立管橫流向和順流向的展向位移功率譜密度分別在圖10和圖11中給出。圖10為3種工況中上游立管橫流向位移功率譜密度?？梢钥闯?，3種工況中上游立管的橫流向主振模態(tài)均為一階，主振頻率分別為3.3 Hz、3.2 Hz、3.1 Hz，且工況3的振動(dòng)頻率包含較微弱的一階高頻率成分。圖11為3種工況中上游立管順流向位移功率譜密度。工況1是頻率為6.7 Hz的純一階振動(dòng)模態(tài)，工況2、工況3則是主振頻率為8.6 Hz的二階振動(dòng)模態(tài)。隨著浸沒長度的增加，立管整體承受的水動(dòng)力載荷增加，更容易激發(fā)出高階振動(dòng)模態(tài)。高頻率的振動(dòng)會直接影響立管的振動(dòng)位移能否得到充分發(fā)展，這解釋了圖9(b)中3種工況中立管的順流向最大位移均方根的差異。

圖10 上游立管橫流向的展向位移功率譜密度

圖11 上游立管順流向的展向位移功率譜密度

3.2.2 下游立管振動(dòng)響應(yīng)

圖12(a)和圖12(b)分別給出了下游立管橫流向及順流向的無量綱化位移均方根沿立管展向的分布。從圖12(a)可以看出，三種工況中下游立管的橫流向主振模態(tài)均為一階，且由于工況2和工況3中下游立管整體承受的水動(dòng)力載荷增加，立管的最大位移均方根均達(dá)到了0.15D，約為工況1的兩倍；關(guān)于順流向的振動(dòng)響應(yīng)結(jié)果，工況1中下游立管仍表現(xiàn)出一階振動(dòng)，而工況2和工況3則均呈現(xiàn)出二階主振模態(tài)。與上游立管順流向的位移響應(yīng)類似，工況2和工況3的順流向最大位移均方根均小于工況1，僅為0.025D，約為工況1的三分之一。這是由于工況2和工況3中下游立管的主振頻率較高，其順流向的周期性振動(dòng)來不及充分發(fā)展便已進(jìn)入下一周期。

圖12 下游立管無量綱化位移均方根

綜合比較圖9和圖12可以發(fā)現(xiàn)，3種工況中上下游立管在橫流向及順流向的主振模態(tài)均保持一致。這可能是由于計(jì)算工況中流速設(shè)置的較小，對應(yīng)的雷諾數(shù)僅為4 800，上下游立管的瀉渦差異未得到充分激發(fā)。同時(shí)，由于上下游立管的間距較小，下游立管同時(shí)經(jīng)受遠(yuǎn)方來流及上游立管尾流的作用，在3種工況中，下游立管的橫流向及順流向的位移均方根均大于上游立管。

圖13和圖14分別給出了下游立管橫流向及順流向的展向位移功率譜密度。3種工況中下游立管的橫流向主振模態(tài)均為一階，主振頻率分別為3.4 Hz、3.2 Hz、3.1 Hz。從圖14可以看出，工況1中下游立管順流向表現(xiàn)為頻率為6.8 Hz的單一頻率一階振動(dòng)模態(tài)，工況2和工況3中下游立管順流向則體現(xiàn)出了多模態(tài)振動(dòng)特性。工況2和工況3的主振模態(tài)均為二階，且伴隨有較高比重的一階頻率成分。工況3中各階模態(tài)的頻率范圍較工況2更廣，這是由于立管浸沒長度的增加使得立管承受的水動(dòng)力載荷增加，且下游立管受到上游立管尾流的影響更大。

綜合分析圖10、圖11、圖13以及圖14，隨著浸沒長度的增加，上下游立管受到的水動(dòng)力載荷增加，立管順流向的主振模態(tài)也會增大。同時(shí)，下游立管受到上游立管尾流的影響，下游立管在順流向的振動(dòng)更容易表現(xiàn)出多模態(tài)振動(dòng)特性。

圖13 下游立管橫流向的展向位移功率譜密度

圖14 下游立管順流向的展向位移功率譜密度

為進(jìn)一步說明下游立管順流向振動(dòng)的多模態(tài)特性，圖15和圖16分別給出了工況2和工況3中下游立管在不同時(shí)段的順流向位移時(shí)空云圖，圖中實(shí)線等值線代表正值，虛線等值線代表負(fù)值。圖15(a)中下游立管表現(xiàn)為純二階振動(dòng)，圖15(b)中則表現(xiàn)為純一階振動(dòng)，而在圖15(c)中則同時(shí)出現(xiàn)了二階振動(dòng)與一階振動(dòng)?？梢钥闯?，工況2中下游立管的順流向振動(dòng)表現(xiàn)出明顯的多模態(tài)振動(dòng)特性，隨著時(shí)間的推移，立管的主振模態(tài)在二階與一階間切換。在圖16給出的工況3中下游立管的順流向位移時(shí)空云圖中，同樣觀察到了多模態(tài)振動(dòng)的特性。圖15和圖16給出的下游立管順流向振動(dòng)的多模態(tài)特性與圖14(b)、圖14(c)的結(jié)果一致。工況2和工況3中下游立管順流向的多模態(tài)振動(dòng)特性可以解釋為，浸沒長度的增加導(dǎo)致上游立管的尾流對下游立管的瀉渦影響較大，使得下游立管瀉渦頻率不穩(wěn)定，進(jìn)而造成了下游立管順流向的多模態(tài)振動(dòng)。

4 結(jié) 語

基于深海立管流固耦合求解器viv-FOAM-SJTU，對階梯流中串列雙立管渦激振動(dòng)開展數(shù)值模擬。首先針對試驗(yàn)工況開展驗(yàn)證計(jì)算，數(shù)值結(jié)果與試驗(yàn)結(jié)果吻合良好，驗(yàn)證了求解器的有效性。為研究不同浸沒長度下串列雙立管的渦激振動(dòng)響應(yīng)，選取了浸沒長度分別為0.40L、0.55L和0.70L的3種工況展開計(jì)算，計(jì)算模型與模型試驗(yàn)保持一致，流速設(shè)置為0.3 m/s。通過對3種工況中上下游立管的位移響應(yīng)及頻率響應(yīng)進(jìn)行分析，得出了以下結(jié)論：

1)隨著浸沒長度的增加，立管整體承受的水動(dòng)力載荷增加，上下游立管的順流向主振模態(tài)均從工況1中的一階轉(zhuǎn)變?yōu)楣r2和工況3中的二階。

2)在橫流向振動(dòng)均為一階的情況下，由于工況2和工況3中立管承受的水動(dòng)力載荷較大，上下游立管的橫流向位移均方根值均大于工況1；關(guān)于順流向振動(dòng)，由于工況2和工況3中立管的主振模態(tài)均為二階，立管的周期性振動(dòng)來不及充分發(fā)展便進(jìn)入下一振動(dòng)周期，因此工況2和工況3中上下游立管的順流向位移均方根值均小于工況1。

3)隨浸沒長度的增加，下游立管的瀉渦受上游立管的尾流作用加強(qiáng)，瀉渦頻率變得不穩(wěn)定，使得工況2和工況3中下游立管的順流向振動(dòng)表現(xiàn)出了明顯的多模態(tài)振動(dòng)特性。