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

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

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

1 前言

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

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

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

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

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

2 數(shù)值模型

2.1 流體運(yùn)動(dòng)控制方程

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

(1)

(2)

(3)

(4)

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

2.2 泥沙運(yùn)動(dòng)控制方程

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

(5)

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

(6)

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

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

(7)

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

(8)

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

(9)

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

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

(10)

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

2.3 管道運(yùn)動(dòng)方程

管道橫向振動(dòng)方程可表示為

(11)

表1 計(jì)算參數(shù)表

3 數(shù)值模型驗(yàn)證

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

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

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

圖1 管道局部沖刷剖面圖

Fig.1 Sectional scouring section of the pipeline

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

表2 計(jì)算參數(shù)

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

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

4.1 振動(dòng)位移

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

圖2 管道位移時(shí)間歷程線

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

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

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

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

4.2 沖刷深度及剖面

圖3 管道局部沖刷剖面

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

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

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

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

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

4.3 渦脫落形式各個(gè)階段的尾渦脫落形式

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

4-a min 4-b zero

4-c max 4-d zero

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

Fig.4 Local vortex shedding of fixed pipeline

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

5-a min 5-b zero

5-c max 5-d zero

圖5 振動(dòng)管道局部渦脫落情況

Fig.5 Local vortex shedding of vibrated pipeline

5 結(jié)論

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

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

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

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