不同浸沒長度下串列雙立管渦激振動數值模擬

武 磊，趙偉文，萬德成

(上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室 船海計算水動力學研究中心(CMHL)，上海 200240)

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

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

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

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

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

1 數值方法

1.1 切片理論

圖1 切片模型

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

1.2 流場控制方程

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

(1)

(2)

1.3 結構動力學控制方程

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

(3)

(4)

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

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

(5)

(6)

式中：M、C、K分別為質量矩陣、阻尼矩陣和剛度矩陣；x、y分別為立管結構順流向和橫流向的位移響應，變量上方的圓點代表對時間的導數；Fx、Fy分別為立管順流向和橫流向所受的水動力載荷。阻尼矩陣C基于瑞利阻尼模型獲取，表述如下：

C=αM+βK

(7)

圖2 流固耦合求解流程

(8)

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

2 問題描述

文中的數值計算模型參考Huera Huarte等[3]開展的階梯流中串列雙立管渦激振動系列試驗。由于模型試驗只關注下游立管的振動響應，只在下游立管上安裝測量儀器，因此雙立管的結構參數略有不同，詳細參數如表1所示(上標a代表上游立管)。

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

表1 立管模型的主要結構參數

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

結構動力學計算部分，將立管模型均勻劃為80個結構單元，總共81個結構節點。其中，處于均勻流場中的結構單元數隨浸沒長度均勻變化。

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

圖4 切片網格劃分

3 計算結果及分析

3.1 試驗工況驗證計算

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

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

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

圖6 下游立管中間節點橫流向振動時歷曲線

圖7 下游立管的瞬時振動輪廓疊加圖((a)、(b)、(c)為試驗結果，(d)、(e)、(f)為數值結果)

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

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

圖8 下游立管中間節點的位移功率譜密度

3.2 浸沒長度對立管渦激振動的影響

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

3.2.1 上游立管振動響應

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

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

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

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

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

3.2.2 下游立管振動響應

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

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

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

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

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

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

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

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

4 結 語

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

1)隨著浸沒長度的增加，立管整體承受的水動力載荷增加，上下游立管的順流向主振模態均從工況1中的一階轉變為工況2和工況3中的二階。

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

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