剪切流下海洋立管渦激振動的三維數(shù)值模擬

羅冬冬，朱仁慶

(江蘇科技大學(xué) 艦舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江212003)

0 引 言

近年來，由于地球人口急劇增加和能源的不斷消耗，各國能源危機(jī)越發(fā)明顯。同時(shí)伴隨陸上石油資源的日益枯竭，使得世界油氣的開發(fā)重點(diǎn)逐步向海洋轉(zhuǎn)移，尤其是深海區(qū)域。海洋立管在深海海洋平臺中應(yīng)用最為廣泛，是連接海底的資源與海上作業(yè)平臺的關(guān)鍵設(shè)備，但也是薄弱易損的構(gòu)件之一。一旦采油立管在渦激振動作用下發(fā)生斷裂，帶來的不僅僅是經(jīng)濟(jì)上的巨大損失，嚴(yán)重的是將會造成巨大的環(huán)境污染?？偠灾?，目前對海洋立管進(jìn)行VIV 問題的研究已成為國內(nèi)外的熱點(diǎn)之一。

經(jīng)過數(shù)十年來對渦激振動的研究，學(xué)者們雖然還未完全把握渦激振動的機(jī)理，但取得了許多階段性的成果，構(gòu)成了當(dāng)今渦激振動研究的基石。無論是實(shí)驗(yàn)研究還是數(shù)值模擬，大多數(shù)的研究主要針對在均勻來流情況下的渦激振動，如張建僑［7］研究質(zhì)量比對柔性細(xì)長立管渦激振動的影響。在同一流速下，質(zhì)量比大的立管模型響應(yīng)位移小，其中空管的渦激振動響應(yīng)一直處于大振幅的鎖定狀態(tài)下。在相同流速下，質(zhì)量比大的立管模型所激起的模態(tài)更高。王嘉松等［4］對1 根長502.92 m，直徑為0.534 m 的細(xì)長立管的進(jìn)行模擬，指出了當(dāng)立管發(fā)生大的、非均勻彎曲變形時(shí)，多模態(tài)的振動和線性動力響應(yīng)不應(yīng)該被忽視。同時(shí)表明預(yù)張力的增加，能有效的減小響應(yīng)振幅，從而有效地減小疲勞損傷。但真實(shí)海洋來流情況往往比較復(fù)雜，不可能出現(xiàn)這種均勻來流的情況，相反更類似于具有一定斜率的剪切流的情況。因此剪切流下的立管渦激振動的模擬更加貼近海洋中的真實(shí)情況，符合實(shí)際的工程需要。目前，對于剪切流下的立管渦激振動的研究還十分少，尤其在國內(nèi)無論是實(shí)驗(yàn)研究還是數(shù)值模擬都十分稀少。在國外，2003 年ExxonMobil［1］公司在MARINTEK 的深水水池中利用旋轉(zhuǎn)的架子進(jìn)行了均勻流和剪切流的一系列的立管渦激振動系列實(shí)驗(yàn)，說明了只有在低階模態(tài)參與的響應(yīng)中，軸向疲勞損傷十分重要。近年來，Remi Bourguet［3］研究了立管在剪切流下的鎖定情況。結(jié)果表明在所有模擬條件中，在渦的發(fā)放與垂向振動局部同步時(shí)，在高流速區(qū)域并且長度至少超過圓柱長度的30%，才會發(fā)生鎖定現(xiàn)象。在鎖定情況下，振幅與在相似雷偌數(shù)下受迫振動的剛性圓柱的局部頻率有關(guān)。梁勇［12］基于Matlab 編寫了預(yù)測渦激振動響應(yīng)的計(jì)算程序，來模擬剪切流下的立管。結(jié)果表明約束條件和立管密度對立管的振動響應(yīng)有著十分重要的影響。

本文利用Workbench 平臺，建立三維流固耦合模型，自編了進(jìn)口剪切流的程序，采用了k － ω 模型，結(jié)合基于薄殼理論的有限元方法，并通過一種新的方法System Coupling 實(shí)現(xiàn)流－固耦合交界面的數(shù)據(jù)交換，從而實(shí)現(xiàn)剪切流下三維流固耦合模擬。與相關(guān)文獻(xiàn)進(jìn)行對比，研究立管的振動變形響應(yīng)及尾渦形態(tài)特征等問題。

1 數(shù)值模擬方法概述

1.1 兩場控制方程

1.1.1 流體場控制方程

流場海水可看做粘性不可壓縮的流體，流場滿足以下微分方程:

1)連續(xù)方程

連續(xù)方程是質(zhì)量守恒定律在流體力學(xué)中的表現(xiàn)形式，該方程可以表述為:單位時(shí)間內(nèi)流入該微元的凈質(zhì)量等于單位時(shí)間流出該微元的凈質(zhì)量。對于三維流動，其平面直角坐標(biāo)形式可表示如下:

式中u，v，w 分別為流體在x，y，z 方向上的速度，m/s。

2)動量方程

動量方程是動量守恒定律在流體力學(xué)中的表現(xiàn)形式，其本質(zhì)是滿足牛頓第二定律。該方程可以表述為:對于一給定的流體微元，其具有的動量對時(shí)間的導(dǎo)數(shù)等于作用在該微元上的各種力之和。其數(shù)學(xué)表達(dá)形式為Navier －Stokes 方程(簡稱N －S)方程，對于不可壓縮流體，寫成Cartesian 坐標(biāo)形式為:

式中:f 為流體在x，y，z 各方向上所受質(zhì)量力，m/s2;p 為流體壓強(qiáng)，Pa;ρ 為密度，kg/m3;ν 為運(yùn)動粘性系數(shù)，m2/s。

1.1.2 固體場控制方程

對海洋立管的動力響應(yīng)采用基于三維實(shí)體單元的有限元方法進(jìn)行模擬。對運(yùn)動方程進(jìn)行有限元離散后可得到有限元方程:

2 算例分析

2.1 模型的建立

本文模型參數(shù)如表1 所示。

表1 立管模型參數(shù)Tab.1 Parameters of riser model

進(jìn)口速度為剪切流v=(0.03 × Z +0.2)m/s，Z為立管長度，立管軸線與流場速度入口邊界距離為10 D，與出口邊界距離為30 D，與兩側(cè)對稱邊界距離為10 D?？紤]的現(xiàn)有的計(jì)算機(jī)能力，使用較精密的平面網(wǎng)格體系，如圖1 所示。

圖1 二維網(wǎng)格拓?fù)鋱DFig.1 Two dimensional grid

流場采用k － ω 模型，迭代方法選用具有高穩(wěn)定性的SIMPLEC 算法。選取立管的圓柱面為流固耦合面，受到動網(wǎng)格的限制，時(shí)間的離散方式選用了一階隱式(1nd－Order－Implicit)。時(shí)間步長為t=0.002 s，計(jì)算總時(shí)長為6 s。

圖2 結(jié)構(gòu)有限元網(wǎng)格圖Fig.2 Finite element of structure

對于結(jié)構(gòu)模型的網(wǎng)格劃分如圖2 所示，同樣結(jié)構(gòu)中，設(shè)定立管的圓柱面為流固耦合面。同時(shí)，為了保持迭代耦合的一致性，時(shí)間步長與時(shí)長同流場相一致。

2.2 海洋立管的運(yùn)動分析

2.2.1 順流向的運(yùn)動分析

在順流向方向，立管因?yàn)槭艿酵弦妨Φ挠绊懀l(fā)生彎曲變形，如圖3 所示，在t = 2.57 s 時(shí)，彎曲達(dá)到最大，立管的最大彎曲變形處發(fā)生在約z/L=0.6 處，最大變形為1.2 D，此后，流場趨于穩(wěn)定，立管變形很小，只是在最大位置處產(chǎn)生小幅振動，振動幅值為0.1 D。如圖4 所示。與文獻(xiàn)［11］對比發(fā)現(xiàn):其運(yùn)動趨勢與均勻流下立管的運(yùn)動趨勢相同，并都多出現(xiàn)了順流向的不可忽略的振動，只有最大彎曲變形位置向上移動。該結(jié)果說明由于剪切流流速遞增的特性，立管受到水流的拖曳力的中心向上移動，使得立管的最大彎曲變形的位置上移，但順流向上對立管的運(yùn)動趨勢的影響與均勻流下相同。同時(shí)說明了System Coupling 實(shí)現(xiàn)流固耦合具有一定的可行性。

圖3 順流向運(yùn)動效果圖Fig.3 In-line motion of riser

圖4 最大位移處點(diǎn)的位移隨時(shí)間的變化圖Fig.4 Scatter diagram of displacement varying with time of point with largest displacement

2.2.2 橫流向的運(yùn)動分析

從圖5 可看出，在橫流向上，立管的振動在1階、2 階振動模態(tài)之間往復(fù)，立管的彎曲變形呈非對稱性，當(dāng)流體趨于穩(wěn)定時(shí)，最終鎖定在2 階振動模態(tài)，同時(shí)出現(xiàn)了3 階振動模態(tài)的情況。

而根據(jù)文獻(xiàn)［11］得到立管自振頻率公式為:

式中:md為單位長度隔水管排開海水的質(zhì)量;CA為附加質(zhì)量系數(shù)，一般取1.0;I 為截面慣性矩;T 為預(yù)緊力。同時(shí)根據(jù)“鎖定”的定義:

圖5 立管不同時(shí)刻橫流向的振動效果圖Fig.5 Cross-flow vibration shape of riser at different moment

式中:f 為隔水管實(shí)際振動頻率;fn為隔水管理論自振頻率;fst為渦脫落頻率。

由此得到前3 階振動模態(tài)的激發(fā)流速:1 階模態(tài)的激發(fā)流速0.188 m/s;2 階和3 階模態(tài)的激發(fā)流速為0.389 m/s 和0.607 m/s。而本文的最高流速為0.488 9 m/s，最低流速為0.2 m/s，并沒有達(dá)到3 階模態(tài)的激發(fā)流速。本文激發(fā)2 階振動模態(tài)的來流約占立管長度的40%，符合上文Remi Bourguet 的指出:在渦的發(fā)放與垂向振動局部同步時(shí)，在高流速區(qū)域并且長度至少超過圓柱長度的30%，鎖定才會發(fā)生。所以立管的振動最終鎖定在2 階振動模態(tài)。同時(shí)激發(fā)1階振動模態(tài)的來流約占立管長度的60%，對立管的影響也很大，因此可能導(dǎo)致了3 階振動模態(tài)的出現(xiàn)。該現(xiàn)象表明:剪切流增加了立管振動的復(fù)雜性，立管出現(xiàn)了比剪切流最大流速下的振動模態(tài)更高1 階的3 階模態(tài)振動，但最終鎖定模態(tài)處于2 階振動模態(tài)。

2.3 海洋立管的流場分析

圖6 為立管在前期和后期不同截面的尾流區(qū)渦的三維特性圖。圖中5 個(gè)截面分別對應(yīng)Z/L=0.1，0.3，0.5，0.7，0.9 的渦量圖。在t = 1.2 時(shí)，由于剪切流的關(guān)系，在Z/L=0.7 和0.9 的截面上都形成了渦的發(fā)放，而其他3 個(gè)界面渦還沒有形成。當(dāng)最終整個(gè)立管的VIV 穩(wěn)定后，渦的發(fā)放形式主要是“2S”形的渦。從圖中可以看出:由于立管的運(yùn)動和剪切流的作用，渦在立管不同的位置有明顯的區(qū)別，在Z/L=0.1，0.5，0.9 處，立管在橫流向的運(yùn)動范圍比較小;兩渦在橫流向上的距離明顯小于Z/L=0.3，0.7 處。

在立管渦的發(fā)放穩(wěn)定前，還伴隨其他不同形式的渦，如圖7 所示。圖7 為t = 2.4 s 時(shí)，Z/L=0.3、0.5 的截面的渦量圖。圖7(a)為Z/L=0.3 截面渦量圖，此時(shí)渦為“P +S”的形式。而到圖7(b)時(shí)即Z/L =0.5 截面，渦為“2P”的形式。

圖6 立管不同截面不同時(shí)刻的渦量圖Fig.6 Vorticity of riser at different sections for different moment

圖7 t=2.4 s 時(shí)渦量圖Fig.7 Vorticity of riser at t=2.4 s

3 結(jié) 語

本文基于計(jì)算流體力學(xué)(CFD)和計(jì)算結(jié)構(gòu)力學(xué)(CSD)方法，借助商業(yè)軟件Ansys 14.0，參考長徑比為482 的均勻流的立管模擬，實(shí)現(xiàn)了VIV 的三維數(shù)值模擬。結(jié)果表明，數(shù)值模擬與均勻流有部分相似;大長徑比柔性隔水管發(fā)生VIV 時(shí)流場受到明顯擾動，渦的脫落呈現(xiàn)如2P、2S、P+S 等模態(tài);立管最大位移處向上移動，立管在橫流向上出現(xiàn)了大的彎曲變形，出現(xiàn)了1 階、2 階模態(tài)之間的轉(zhuǎn)變，最終鎖定在2 階模態(tài)。

［1］TOGNARALLI M A．VIV response of a long fexible cylinder in uniform and linearly sheared currents［C］//Offshore Technology Conference，2004．

［2］FARSHIDIANFAR A，ZANGANEH H． A modified wake oscillator model for vortex-induced vibration of circular cylinders for a wide range of mass － damping ratio［J］．Journal of Fluids and Structures，2010．

［3］BAHMANI M H，AKBARI M H． Effects of mass and damping ratios on VIV of a circular cylinder［J］． Ocean Engineering，2010．

［4］REMI B，MICHAEL S T． Lock-in of the vortex-induced vibrations of a long tensioned beam in shear flow［J］．Journal of Fluids and Structures，2011．

［5］WANG Jia-song，ZHAN Lu-lu． Three dimensional numerical simulation of vortex induced vibration for a 500 － m － long marine riser［C］//ISOPE，2012．

［6］LIN Li-ming，ZHONG Xing-fu，WU Ying-xiang． Vortexinduced vibrations of a circular cylinder with different geometric disturbances［C］//ISOPE，2012．

［7］HUANG Zhi-yong，CUI Wei-cheng．Vortex-induced vibration of flexible riser in time domain［D］． Shanghai:Shanghai Jiao Tong University，2008．

［8］ZHANG Jian-qiao，SONG Ji-ning，LV Lin，et al．Analysis on the vortexoinduced vibrations experimental results of flexible risers［J］． China Offshore Platform，2009．

［9］TANG Guo-qiang，TENG Bin，LV Lin，et al．Development and application of an empirical model for vortex-induced vibration of flexible deep water riser［J］．China Offshore Platform，2010．

［10］CAI Jie． Investigation on numerical simulation of the fluidsolidInteraction characteristics for oceanic riser［D］．Shanghai:Shanghai Jiao Tong University，2010．

［11］CAI Jie，YOU Yun-xiang，LI Wei，et al． The VIV characteristics of deep-sea risers with high aspect ratio in a uniform current profile［J］．Chinese Journal of Hydrodynamics，2010．

［12］WANG Cheng-guan，TIAN Zhong-xu，QIAO Xin-qi． Study on three dimensional numerical simulation of VIV characteristics for marine riser［J］． Chinese Journal of Hydrodynamics，2011．

［13］LIANG Yong．Response analysis of flexible riserssubjeted to Vortex-induced vibration in sheared flow［D］． Harbin:Harbin Institute of Technology，2011．