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

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

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

1 前言

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

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

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

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

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

2 數(shù)值模型

2.1 流體運動控制方程

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

(1)

(2)

(3)

(4)

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

2.2 泥沙運動控制方程

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

(5)

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

(6)

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

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

(7)

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

(8)

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

(9)

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

通過如下的方法對底床坐標進行更新

(10)

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

2.3 管道運動方程

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

(11)

表1 計算參數(shù)表

3 數(shù)值模型驗證

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

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

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ù)值模擬結果與實驗結果吻合良好。從而證明本文建立的模型可以用來對振蕩管道周圍的沖刷過程進行預測。對比圖1-a和圖1-b可以發(fā)現(xiàn)，從t= 30 min到t= 243 min管道下方的最大沖深增加不多。但是管道周圍的沖刷剖面改變較大。這與要體現(xiàn)在(a)管道下游剖面更加平緩(b)管道上游和下游的過渡段增長。

表2 計算參數(shù)

4 計算結果與分析

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

4.1 振動位移

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

圖2 管道位移時間歷程線

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

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

(3)大約450～1 500 s，管道正向振幅近乎不變而負向逐步增大到約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分別對應著固定管道所受升力最小值、零、最大值和零時刻管道局部渦脫落情況。圖5給出了振動管道局部渦脫落情況。圖5-a～5-d分別對應著振動管道位移最小值、零、最大值和零時刻管道局部渦脫落情況。

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和負向渦B。從圖5-b中可知，當管道繼續(xù)向上運動時，負向渦B向下運動并變長，正向渦A從管道上脫落下來。當管道運動到位移最大處時，負向渦B分裂成渦B和渦B1。如圖5-c和5-d所示，當管道向下運動到原位置時，負向渦B脫落。從而在管道后方形成由渦A1和渦B組成的一排渦以及由渦A和渦B1組成的另一排渦。在本文的計算條件下，從整體上講，振動管道局部的渦長度較固定情況下短，這使得振動管道后方底床切應力大于固定管道的情況，進而導致振動管道后方?jīng)_刷深度較固定管道情況下深。

5-a min 5-b zero

5-c max 5-d zero

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

Fig.5 Local vortex shedding of vibrated pipeline

5 結論

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

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

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

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