海底振動管道局部泥沙沖刷數(shù)值研究

殷 俊，劉名名，冀 昊，金 鑫，郭曉玲

(1.國家電投集團(tuán)廣東電力有限公司，廣州 510130；2.成都理工大學(xué) 能源學(xué)院，成都 610059;3.水電水利規(guī)劃設(shè)計總院，北京 100120;4.國家知識產(chǎn)權(quán)局專利局專利審查協(xié)作天津中心，天津 300304)

海底管道作為一種高效、經(jīng)濟(jì)的海上輸運方式，在海洋油氣資源開發(fā)中得到廣泛應(yīng)用。海底管道周圍的泥沙顆粒在水流作用下將產(chǎn)生輸運，泥沙顆粒的運移會在管道周圍形成沖刷坑。沖坑的存在將嚴(yán)重影響海底管道的在位穩(wěn)定性。

1 前言

長期以來研究者主要以物理實驗和數(shù)值模擬兩種方法對海底管道局部沖刷問題進(jìn)行研究。受限于計算機(jī)模擬技術(shù)，早期的研究主要集中在物理實驗研究方面，如Mao[1]，Ibrahim和Nalluri[2]，Sumer等[3]，F(xiàn)reds?e等[4]以及Chiew等[5]的研究工作。上述試驗結(jié)果表明：海底管道局部沖刷特性與來流流速、管道直徑、水深以及泥沙粒徑密切相關(guān)，并建立了相關(guān)經(jīng)驗公式預(yù)報海底管道的沖刷特性。這些開拓性研究工作，極大地加深了人們對海底管線局部沖刷的認(rèn)識。

20世紀(jì)80年代起，眾多學(xué)者開始采用數(shù)值模擬對海底管道局部沖刷問題開展研究工作。Hansen等[6]基于無粘、無旋假設(shè)，建立了海底管道局部沖刷勢流數(shù)值分析模型，由于無法考慮流體的粘性效應(yīng)，導(dǎo)致數(shù)值模擬結(jié)果與試驗結(jié)果存在較大差異。為克服勢流理論數(shù)值分析模型缺點，Beek和Wind[7]建立了基于N-S方程的數(shù)值分析模型，結(jié)合標(biāo)準(zhǔn)k-ε湍流封閉模型，實現(xiàn)了對海底管道附近泥沙輸運的模擬，模擬結(jié)果表明：數(shù)值模擬得到的數(shù)據(jù)在沖刷趨勢上與實驗室結(jié)果吻合較好，但是數(shù)值結(jié)果與實驗室結(jié)果在管道周圍底床剖面形態(tài)上存在較大差異。此后，Li和Chengi[8]采用N-S方程和大渦模擬的方法，對海底管道局部泥沙沖刷問題進(jìn)行了數(shù)值研究。在后續(xù)研究工作中，Li和Cheng[9]將上述數(shù)值模型拓展到管道自沉的數(shù)值模擬研究工作，相關(guān)數(shù)值分析模型得到了令人滿意的結(jié)果。此后，Liu等[10]以及劉名名等[11]成功通過求解利用SSTk-ω湍流模型封閉的Navier-Stokes方程組，實現(xiàn)了對海底管道在波浪作用下產(chǎn)生的局部泥沙沖刷的模擬。相關(guān)的研究結(jié)果表明，SSTk-ω湍流模型可以實現(xiàn)對渦流場及近壁區(qū)高流速梯度的準(zhǔn)確模擬，進(jìn)而保障海底管道局部泥沙沖刷的準(zhǔn)確性。

但是實際的情況為：海底管道發(fā)生局部沖刷后，沖坑會沿著管道軸向方向不斷擴(kuò)展，從而形成懸空管道。在周期性交變流體力作用下，懸空管道會發(fā)生渦激振動，尤其當(dāng)管道后方的渦脫落頻率接近于懸跨段的固有頻率時，會發(fā)生“鎖定”(Lock-in)現(xiàn)象。目前，針對振動管道局部泥沙沖刷的研究工作比較少，并且研究對象均為實驗室小直徑管道。開拓性研究工作為Sumer等[3]的試驗研究工作。在Sumer等[3]的試驗中，管道只允許發(fā)生橫流向振動，研究其周圍泥沙沖刷狀況。在此基礎(chǔ)上，Shen等[12]對具有x、y兩個自由度的海底管道在水流作用下的局部泥沙沖刷問題進(jìn)行了實驗研究。Yang等[13]通過實驗手段對管道振動與泥沙局部沖刷的耦合作用進(jìn)行了研究。Gao等[14]和沙勇等[15-16]同樣利用實驗手段對振動管道在水流作用下的局部泥沙沖刷問題進(jìn)行了研究。張瑋等[17]和孫國民等[18]也對水流作用下的海底管道沖刷問題開展了研究。

需要指出的是：上述物理模型實驗由于難以得到詳細(xì)的流場流動結(jié)構(gòu)，所以未能給出流場、管道振動以及泥沙輸運之間的耦合作用機(jī)理。Zhao和Cheng[19]建立起一套通過求解不可壓縮粘性牛頓流體的Navier-Stokes方程、簡化的彈簧阻尼約束下的結(jié)構(gòu)運動方程以及泥沙顆粒運動方程的有限元模型。通過該模型對實驗室尺寸的振動管道在水流作用下泥沙局部沖刷問題展開了數(shù)值模擬。研究結(jié)果表明：管道振動引起的最大沖刷深度較固定管道情況下大25%。

從以上的內(nèi)容可以看出，振動管道局部泥沙沖刷較固定管道情況下局部泥沙沖刷有較大的差別。但是在現(xiàn)在已發(fā)表的文獻(xiàn)中僅有對實驗?zāi)Ｐ统叽绻艿谰植磕嗌硾_刷的研究。而對于原型尺寸振動管道局部泥沙沖刷的研究工作尚未見報道。本文將對原型尺寸下管道振動、泥沙顆粒運動以及底床變形耦合問題展開研究工作。

2 數(shù)值模型

2.1 流體運動控制方程

本文采用的流體運動控制方程為Navier-Stokes方程組。為描述由于管道振動以及底床變形引起的網(wǎng)格位置變化，本文采用ALE方法對變形網(wǎng)格位置進(jìn)行追蹤，ALE參考坐標(biāo)系下流體運動控制方程可以表示為

(1)

(2)

(3)

(4)

式中：υt為湍流引起的粘性系數(shù)，k為特定的湍動能，δij為Kronecker函數(shù)。本文采用SST(Shear-Stress Transport)k-ω湍流模型求解k。具體方法及參數(shù)設(shè)置可參考Menter[20]以及Menter[21]等的文獻(xiàn)。

2.2 泥沙運動控制方程

本文采用的泥沙沖刷模型，同時考慮推移質(zhì)和懸移質(zhì)輸沙對沖刷的影響。在懸移質(zhì)輸沙模型中，泥沙濃度的擴(kuò)散控制方程為

(5)

式中：c為泥沙在水體中的濃度，ws為泥沙在水中的沉降速度，σc為常數(shù)，本文取σc= 1.0。泥沙在水中的沉降速度通過下式計算

(6)

式中：ds為計算采用泥沙的中值粒徑，D*為無因次后泥沙顆粒的粒徑，可以寫為

D*=[g(s-1)/v2]1/3ds

(7)

推移質(zhì)輸沙率采用如下的計算公式，相關(guān)表達(dá)式為

(8)

式中：qb為計算域中單寬輸沙率，s=ρs/ρ為泥沙顆粒的密度與水流密度的比值，ρs為計算采用的泥沙顆粒的密度。Shields參數(shù)(θ)可以寫成如下的形式

(9)

式中：u*為底摩阻流速。

通過如下的方法對底床坐標(biāo)進(jìn)行更新

(10)

式中：yb為海床的垂向坐標(biāo)值，ps為計算采用的泥沙顆粒的松散度。系數(shù)C在本文中取為2.0。

2.3 管道運動方程

管道橫向振動方程可表示為

(11)

表1 計算參數(shù)表

3 數(shù)值模型驗證

Sumer等[3]對水流作用下振動海底管道的局部沖刷問題進(jìn)行了實驗研究。為驗證本文所建數(shù)值模型的正確性，本文首先對單向流作用下振蕩管道周圍的沖刷問題進(jìn)行計算，并與實驗結(jié)果進(jìn)行了比較。相關(guān)計算參數(shù)如表1所示。

圖1給出了不同時刻管道周圍泥沙沖刷剖面圖。圖1-a給出了t= 30 min時，管道周圍泥沙沖刷的剖面形態(tài)。從圖1-a的對比圖中可以發(fā)現(xiàn)，本文所建數(shù)值模型得到的數(shù)值結(jié)果與Sumer等[3]實驗得到的實驗結(jié)果具有較高的吻合性。這種吻合性不僅僅體現(xiàn)在沖刷深度上，同時還體現(xiàn)在沖坑的形態(tài)上。這表明本文所建立的模型能夠較好地預(yù)測振動管道周圍的泥沙沖刷情況，進(jìn)而說明本文所建立數(shù)值模型可以用來對振蕩管道周圍的泥沙沖刷問題進(jìn)行研究。

1-at=30 min 1-bt=243 min

圖1 管道局部沖刷剖面圖

Fig.1 Sectional scouring section of the pipeline

圖1-b給出了t= 243 min時，管道周圍的泥沙沖刷剖面形態(tài)。從圖中可以發(fā)現(xiàn)，數(shù)值模擬結(jié)果與實驗結(jié)果吻合良好。從而證明本文建立的模型可以用來對振蕩管道周圍的沖刷過程進(jìn)行預(yù)測。對比圖1-a和圖1-b可以發(fā)現(xiàn)，從t= 30 min到t= 243 min管道下方的最大沖深增加不多。但是管道周圍的沖刷剖面改變較大。這與要體現(xiàn)在(a)管道下游剖面更加平緩(b)管道上游和下游的過渡段增長。

表2 計算參數(shù)

4 計算結(jié)果與分析

通過上文的驗證可知，本文所建數(shù)值模型能夠較為準(zhǔn)確地預(yù)測振動海底管道局部沖刷。本文對直徑為1.0 m的振動海底管道在1 m/s的流速下的局部泥沙沖刷問題進(jìn)行研究。具體的計算參數(shù)如表2所示。由于管道的振動會導(dǎo)致管道接觸到海床，在本文的計算過程中當(dāng)管道距離海床表面的最小距離小于5 cm時，認(rèn)為管道與海床接觸并強(qiáng)制管道的速度為0 m/s。

4.1 振動位移

圖2給出了振動位移時間歷程線。從圖中可以看出管道的最大振動位移大約為0.8倍的管道直徑。管道位移在時間上可以分為如下的4個階段：

圖2 管道位移時間歷程線

(1)大約0～100 s時，管道振動幅值較小并且逐步增大。

(2)大約100～450 s時，管道振幅逐步增大，最大振幅增大到約為0.8倍管道直徑，并且管道正向振幅明顯大于負(fù)向振幅，在此階段管道運動會接觸底床。

(3)大約450～1 500 s，管道正向振幅近乎不變而負(fù)向逐步增大到約0.7倍管道直徑。在此階段管道運動也會接觸底床。

(4)大于1 500 s，管道振動幅值近乎穩(wěn)定。

4.2 沖刷深度及剖面

圖3 管道局部沖刷剖面

圖3給出了固定管道以及振動管道平衡狀態(tài)下管道局部泥沙沖刷剖面對比圖。從圖中可以看出，在本文的計算范圍內(nèi)振動管道的沖坑形態(tài)較固定管道時有較大的不同，具體體現(xiàn)在以下4個方面。

(1)振動管道最大沖刷深度較固定管道情況下最大可高50%。

(2)振動管道的最大沖刷深度位于管道后方而非管道正下方。

(3)振動管道后方?jīng)_刷范圍和沖刷深度明顯高于固定管道的情況。

(4)振動管道前方?jīng)_刷深度和沖刷范圍較固定管道情況下小。

4.3 渦脫落形式各個階段的尾渦脫落形式

由上文的討論可見振動管道情況下的局部泥沙沖刷較固定管道情況不同。圖4給出了直徑為1.0 m的固定海底管道，在1 m/s的流速下的局部沖刷平衡時，固定管道局部渦脫落情況。圖4-a～4-d分別對應(yīng)著固定管道所受升力最小值、零、最大值和零時刻管道局部渦脫落情況。圖5給出了振動管道局部渦脫落情況。圖5-a～5-d分別對應(yīng)著振動管道位移最小值、零、最大值和零時刻管道局部渦脫落情況。

4-a min 4-b zero

4-c max 4-d zero

圖4 固定管道局部渦脫落情況

Fig.4 Local vortex shedding of fixed pipeline

從圖4中可以看出在一個升力周期內(nèi)，管道后方有一對渦脫落。并且總體上渦是沿著水流方向脫落的。尾渦的最大長度約為管道直徑的2.0倍。從圖5-a中可以發(fā)現(xiàn)，此時管道周圍有正向渦A1脫落，同時存在正向渦A和負(fù)向渦B。從圖5-b中可知，當(dāng)管道繼續(xù)向上運動時，負(fù)向渦B向下運動并變長，正向渦A從管道上脫落下來。當(dāng)管道運動到位移最大處時，負(fù)向渦B分裂成渦B和渦B1。如圖5-c和5-d所示，當(dāng)管道向下運動到原位置時，負(fù)向渦B脫落。從而在管道后方形成由渦A1和渦B組成的一排渦以及由渦A和渦B1組成的另一排渦。在本文的計算條件下，從整體上講，振動管道局部的渦長度較固定情況下短，這使得振動管道后方底床切應(yīng)力大于固定管道的情況，進(jìn)而導(dǎo)致振動管道后方?jīng)_刷深度較固定管道情況下深。

5-a min 5-b zero

5-c max 5-d zero

圖5 振動管道局部渦脫落情況

Fig.5 Local vortex shedding of vibrated pipeline

5 結(jié)論

本文通過求解不可壓縮粘性流體N-S方程、結(jié)合SSTk-ω湍流封閉模模型、泥沙輸運方程以及管道振動方程，建立起振動管道局部泥沙沖刷數(shù)值模型。通過與已發(fā)表數(shù)據(jù)的對比證明本文所建模型的正確性。利用本文所建數(shù)值模型對原型尺度管道局部泥沙沖刷問題進(jìn)行了數(shù)值研究。在本文計算條件下得到如下結(jié)論：

(1)管道的最大振動位移約為0.8倍的管道直徑。

(2)振動管道最大沖刷深度較固定管道情況下最大可高50%。

(3)振動管道后方有兩排旋渦。