李宏偉，王立忠，國 振，袁 峰

(1. 浙江大學 濱海和城市巖土工程研究中心,浙江 杭州 310058；2. 塔里木大學 水利與建筑工程學院， 新疆 阿拉爾 843300)

海底泥流沖擊懸跨管道拖曳力系數(shù)分析

李宏偉1,2，王立忠1，國 振1，袁 峰1

(1. 浙江大學 濱海和城市巖土工程研究中心,浙江 杭州 310058；2. 塔里木大學 水利與建筑工程學院， 新疆 阿拉爾 843300)

海底滑坡是海洋油氣工程最危險的地質(zhì)災害之一，直接影響海底管線運營安全?；麦w失穩(wěn)滑動過程中，由于海水摻和作用逐漸加速轉變成快速滑動泥流，沖擊海底懸跨管道。基于當前國際通用的計算流體動力學(CFD)方法，采用赫巴模型描述快速滑動泥流，計算分析海底滑坡沖擊懸跨管道的受力特性，重點分析懸跨高度對管道法向拖曳力系數(shù)的影響。研究發(fā)現(xiàn)，管道法向拖曳力系數(shù)隨懸跨高度的增長而增大，當達到某一懸跨高度時，管道拖曳力系數(shù)保持穩(wěn)定。

海底滑坡；計算流體動力學(CFD)；泥流；拖曳力；非牛頓流體；海底管線

海底管道具有輸送連續(xù)、管理方便、效率高、成本低等優(yōu)點，是深水油氣輸運的重要生命線，其直徑一般在0.1～1.0 m之間。由于海底管道長距離輸運，管道經(jīng)歷的海底地形崎嶇，起伏變化大，容易形成管道局部懸跨段。而這些懸跨段所在的海溝、海槽等低洼區(qū)域往往是海底滑坡泥流的主要輸運通道。海底滑坡泥流的流速極快，可達30 m/s，其影響面積也很廣。因此，對于懸跨段管道而言，極易遭受滑坡泥石流的快速沖擊作用而破壞[1]，有必要對懸跨段管道遭受快速泥流沖擊時的受力特點展開研究。

Vanneste等[2]指出，目前對海底滑坡方面的研究主要針對滑坡形成的濁流階段，將滑坡土體簡化成黏塑性流體，通過水動力學方法模擬滑坡土體的長距離滑移[3-5]。Coussot等[6]討論了泥漿流溫度、PH值、電解質(zhì)濃度、固體顆粒含量和黏土含量等參數(shù)對泥流流變性能的影響，指出赫巴模型中屈服應力能夠很好估計泥流真實的屈服應力狀態(tài)。Locat等[7]和Jeanjean等[8]對海底滑坡及其影響區(qū)域內(nèi)的海底管道安全性進行了論述。海底滑坡泥流對管道的作用可分為快速泥流的沖擊作用和慢速泥流的推力作用，已有研究大多針對慢速泥流。王立忠等[9]分析了慢速泥流作用下管道受力的一般規(guī)律。Parker等[10]、Yuan等[11]直接將海底滑坡沖擊荷載簡化為定值，分析了海底滑坡對管道整體受力的影響，簡化了海底滑坡快速泥流的沖擊作用。Zakeri等[12-15]開展了海底滑坡泥流沖擊海底管道的模型實驗研究和數(shù)值模擬計算，得出了管道拖曳力系數(shù)的經(jīng)驗公式。Randolph等[16-17]開展了海底滑坡沖擊海底管線的離心機模型試驗研究，提出了基于巖土力學和流體力學方法的管道拖曳力計算方法。

基于計算流體動力學(CFD)方法對海底滑坡泥流正面沖擊懸跨管道的過程進行數(shù)值模擬。在計算中，采用赫巴模型描述滑坡泥流，著重分析了滑坡泥流沖擊不同懸跨高度下管道的拖曳力系數(shù)，得出了拖曳力系數(shù)與懸跨高度(H)之間的關系，提出了更為合理的拖曳力系數(shù)取值方法，為海底管線設計提供參考依據(jù)。

1 滑坡沖擊荷載計算方法

海底滑坡泥流沖擊管道主要體現(xiàn)在拖曳力的大小，而拖曳力系數(shù)是衡量拖曳力大小的關鍵參數(shù)。已有研究主要基于水動力學方法，采用Morison公式計算細長結構物上流體作用力，這是一種半經(jīng)驗半理論的描述。White[18]提出了法向拖曳力系數(shù)(CD)的計算方法：

式中：FD是法向拖曳力，U是滑體上游自由流體流速，ρ是滑坡體流體的密度，A為障礙物的特征面積。

Zakeri等[14，19]采用冪律模型通過數(shù)值(CFD)模擬手段研究了海底滑坡對海底管線的沖擊荷載，對室內(nèi)試驗結果進行了反演分析，提出了海底滑坡泥流正面沖擊管道時產(chǎn)生的法向拖曳力系數(shù)如式(2)、(3)。

管底距離海床1D處的管道：

管底距離海床2 mm處的管道：

式中：Renon-Newtonian是非牛頓流體雷諾數(shù)，如式(4)所示。

其中，τ是滑坡泥流剪切應力。

2 數(shù)值模型

2.1水和泥流控制方程

采用ANSYS CFX不可壓縮兩相流(水和泥流)進行模擬計算，不同相采用小寫希臘字母α和β表示，總相數(shù)為Np=2。其控制方程如下：

1)連續(xù)性方程

2)動量方程

式中：Pα和μα是α相的壓強和黏性系數(shù)；SMα是外部質(zhì)量力引起的動量源相；Mα是由于其他相引起的作用在α相上的總界面力，由下式可得

包括拖曳力、升力、壁面潤濕力、虛擬質(zhì)量力和湍流耗散力等[14]。

2.2泥流流變本構

目前，研究滑坡泥流流變本構的模型有賓漢模型、赫巴模型和冪律模型，各模型的剪切應力與應變率之間的關系如式(8)、(9)所示。

赫巴模型(herschel-bulkley model):

冪律模型(power-law model):

2.3模型設定

2.3.1 計算域設定

Zakeri等[12]進行了泥流沖擊管道的常重力模型試驗，其中管道直徑采用22.2 mm和28.6 mm兩組類型。由于采用無量綱參數(shù)描述拖曳力系數(shù)，因而這里統(tǒng)一采用計算模型管道直徑為25 mm。圖2所示為一個典型的計算模型，在計算域中，管道距離入口邊界為9.5D，距離出口邊界為13.5D，距離海床邊界為4D，距離計算域上邊界為5D，水與泥流交界面距離計算域上邊界為1D，泥流從左邊入口邊界進入計算域，從右邊出口邊界流出，海底管道懸跨高度為H，H為管道中心到底部邊界的距離。

圖1 流變模型Fig. 1 Rheological model

圖2 幾何模型和數(shù)值模型Fig. 2 Domain geometry and simulation model

采用圖2模型進行數(shù)值模擬，將滑坡泥流正面沖擊懸跨管道模擬成兩相(水和泥流)并考慮浮力連續(xù)相的自由液面流?；履嗔鞑捎每紤]浮力相的層流模型，其本構關系采用赫巴模型；水中產(chǎn)生的湍流采用k-ε模型。在兩相流中，由于水和泥流的密度不同，從而影響流體計算的控制方程(式(5)、式(6))。浮力相SM,Buoyancy,α=(ρα-ρref)g引入動量方程式(6)中，其中g是重力加速度，ρα是滑坡泥流的密度，ρref是浮力參考密度，為水的密度997kg/m3。由于兩相(水和泥流)的速度場不同，因而兩相自由液面間存在界面拖曳力，其通過無量綱參數(shù)界面拖曳力系數(shù)來描述界面拖曳力，其界面拖曳力系數(shù)為2～3[14]。

2.3.2 邊界設定

計算域的邊界如圖3所示，管道表面為無滑移粗糙壁面邊界，其粗糙系數(shù)ks=0.001 5 mm[14]，出口邊界為開放邊界，計算域前后邊界為對稱邊界，其中，對稱邊界認為是無剪切應力的平面應變邊界，計算域上部邊界為自由滑移壁面邊界，下部邊界為無滑移粗糙壁面邊界，其粗糙系數(shù)ks=0.5 mm[14]。

圖3 流體計算域邊界Fig. 3 Fluid domain boundary conditions

2.3.3 網(wǎng)格剖分

采用ANSYS ICEM軟件對整個計算域進行非結構化網(wǎng)格分區(qū)劃分。如圖4所示，對管道周圍網(wǎng)格進行局部細化，這部分區(qū)域是本文研究泥流流動狀態(tài)及流體與管道相互作用的重要區(qū)域。計算域整體最大網(wǎng)格尺寸為25 mm，管道表面網(wǎng)格大小為2 mm逐漸增長到50 mm處的7 mm，管道表面有覆面層網(wǎng)格，需對其分五層，覆面層網(wǎng)格總厚度為2 mm[19]。覆面層內(nèi)為三棱柱網(wǎng)格，以減小管道表面的網(wǎng)格扭曲率。

圖4 管道直徑為25 mm的計算域網(wǎng)格與覆面層網(wǎng)格Fig. 4 The domain grid size and inflated boundary, with 25 mm pipe diameter

2.3.4 求解設定

在數(shù)值模擬中，初始時刻(t=0)需對計算域進行初始化，計算域中水的體積分數(shù)為1，各速度分量為0。時間步的選擇由式(10)確定：

φip=φup+β

3 數(shù)值分析算例

SI[22]的實驗研究了由高嶺土、石英砂和水混合配置而成的泥流的一般特性，給出了冪律模型與赫巴模型的流變參數(shù)，如表1所示。利用表1中六種不同成分的泥漿流體的流變參數(shù)進行模擬計算。

表1 泥流的組成和特性Tab. 1 Slurry flow composition and properties

注：C*表示黏土；W*表示水；S*表示石英砂。

3.1驗證與對比

Zakeri[19]采用表1中的冪律模型，研究了管底距離水槽底部1D、2 mm時的管道拖曳力系數(shù)變化。采用表1中的赫巴模型計算的CD值與Zakeri的實驗及計算結果進行對比。從圖5可以看出，采用赫巴模型計算結果與實驗結果更為接近，尤其是在非牛頓流體雷諾數(shù)小于10的情況。非牛頓流體雷諾數(shù)增大，法向拖曳力系數(shù)減小。擬合公式如式(12)、(13)，兩者擬合R2值分別為99.00%和98.40%。

管底距離海床1D處的管道：

管底距離海床2 mm處的管道：

通過式(12)與(13)可以計算得到每延米管道上產(chǎn)生的法向拖曳力。

由圖5觀察到，當Renon-Newtonian<10時，本文與Zakeri計算結果都高于實驗值。這主要由于在試驗中，管道表面粗糙，凸凹不平，局部存在光滑區(qū)域，泥流流過管道時，兩者接觸面之間存在少量的水泡，使得測出的管道拖曳力較小。而在數(shù)值模擬中，管道表面采用無滑移粗糙表面，在低雷諾數(shù)下，黏性增大，管道上下表面處泥流流速大，管道的拖曳力較大。

圖5 結果比對Fig. 5 Comparing with the results of Zakeri

3.2參數(shù)分析

3.2.1 層流邊界層

首先對不放置管道近海床區(qū)域泥流流動的邊界層發(fā)展情況進行了模擬計算。采用入口流速U=1.0 m/s，泥流為35%Clay，海床粗糙系數(shù)為ks=0.5 mm[14]的情況下，得到了泥流計算域速度矢量圖，如圖6所示?？拷４擦魉贋榱悖∮嬎阌蛑辛魉龠_到0.99U∞的高度作為邊界層的充分發(fā)展高度。

圖6 邊界層流速矢量圖Fig. 6 The velocity vector of boundary layer

由于近海床區(qū)域泥流受到黏性影響，毗鄰海床附近的泥流停滯下來，速度為零，而在垂直于流動方向上形成速度梯度，并使靠近海床的泥流速度減慢，隨著黏滯泥流向前移動，邊界層逐漸增加，更多黏性泥流速度減緩，逐漸形成一穩(wěn)定的層流邊界層。計算域入口流速、泥流密度及海床粗糙度對層流邊界層充分發(fā)展高度有影響，如圖7所示。其中，海床粗糙度主要取決于海床起伏(比如砂脊、砂波的影響)和沉積物粒徑大小，這里未考慮海床起伏的影響，只研究了粒徑大小在粉土(ks=0.001 mm)到粗砂(ks=1 mm)范圍內(nèi)變化對層流邊界層的影響[23]。

圖7 邊界層厚度δFig. 7 Boundary layer thickness

從圖7可以看出，Clay含量越高，泥流的黏度越大，層流邊界層充分發(fā)展的高度越大；泥流速度越大，層流邊界層充分發(fā)展的高度越小；層流邊界層充分發(fā)展高度不受海床粗糙系數(shù)的影響。從圖7(a)、(b)可以看出，計算域入口流速為1.0 m/s、海床邊界粗糙系數(shù)為ks=0.5 mm[14]的情況下，不同Clay含量(10%Clay、15%Clay、20%Clay、25%Clay、30%Clay、35%Clay)的泥流，隨著Clay含量的增高，層流邊界層厚度也逐漸發(fā)展，其厚度在0.4D～1.4D之間變化；圖7(c)、(d)揭示了在35%Clay含量的泥流中，在海床邊界粗糙系數(shù)為ks=0.5 mm[14]的情況下，隨著計算域入口流速的增大，層流邊界層厚度逐漸減小，其厚度在0.4D～1.8D之間變化；從圖7(e)、(f)可以看出，在計算域入口流速為1.0 m/s、35%Clay含量的泥流情況下，邊界層隨海床粗糙系數(shù)的增大而保持在1.24D。

圖8 拖曳力系數(shù)與雷諾數(shù)之間的關系曲線Fig. 8 Drag coefficient vs. Reynolds number curves

3.2.2 懸跨高度比H/D

考慮層流邊界層的影響，對不同懸跨高度(H)的海底管道進行了數(shù)值模擬分析，得到了不同懸跨高度下的管道拖曳力系數(shù)，如圖8所示。從圖8中可以看出，管道拖曳力系數(shù)隨非牛頓流體雷諾數(shù)的增大而減小，隨管道懸跨高度的增高而增大。當非牛頓流體雷諾數(shù)較大時，層流邊界層較小，H/D從0開始增長很小范圍后CD才保持穩(wěn)定；當非牛頓流體雷諾數(shù)較小時，層流邊界層較大，H/D在從0開始增長很大范圍后CD才保持不變。為了方便比較，將H/D=5時的拖曳力系數(shù)CD作為所有工況的C。

圖9 不同雷諾數(shù)下的管周圍流態(tài)示意(H/D=5)Fig. 9 The flow regime around the pipe under different Reynolds number at H/D=5

圖9揭示了海底管道懸跨高度為H=5D的情況下，不同非牛頓流體雷諾數(shù)下管道周圍流態(tài)示意圖。從圖9中可以看出，管道周圍流態(tài)穩(wěn)定，并未受到層流邊界層的影響。當Renon-Newtonian<10時，管道側面的泥流流線不存在分離點；當Renon-Newtonian>10時，管道側面的泥流流線分離點逐漸向泥流上游移動，管道后面區(qū)域出現(xiàn)回流現(xiàn)象。

圖10 海底管道不同懸跨高度速度矢量圖Fig. 10 The velocity vectors of different suspended heights of the submarine pipeline

圖10揭示了在U=1.0 m/s、20%Clay、ks=0.5 mm[14]和Renon-Newtonian=21.3的工況下，管道在不同懸跨高度下泥流的流速矢量圖。由于管道對泥流的阻流作用，泥流在管道迎流面位置速度減緩，繞流過管道。當H/D<2時，管道周圍的流速分布受到海床邊界的影響較大；當H/D≥2時，管道逐漸遠離邊界層，其周圍泥流流速分布不受海床邊界的影響。

由圖11可知，管道拖曳力系數(shù)隨著管道懸跨高度的增加而增大，當懸跨高度比(H/D)增加到2之后，拖曳力系數(shù)保持穩(wěn)定。這里將此懸跨高度比定義為影響拖曳力系數(shù)的臨界懸跨高度比(H/D)Critical。由于滑坡泥流正面沖擊管道時管道周邊只形成很薄的邊界層，影響管道拖曳力系數(shù)的主要為近海床層流邊界層的厚度。因此，臨界懸跨高度比可通過下式簡單計算：

表2給出了不同雷諾數(shù)下層流邊界層的厚度變化。對比圖12，可以發(fā)現(xiàn)，式(14)給出的臨界懸跨高度比較好地預測了管道拖曳力系數(shù)的變化。

表2 不同雷諾數(shù)下邊界層厚度δTab. 2 Boundary layer thickness under different Reynolds Number

圖11 管道拖曳力系數(shù)與懸跨高度的關系(Renon-Newtonian=21.3)Fig. 11 The relationship of the drag coefficient and the suspended height at Renon-Newtonian=21.3

圖12 拖曳力系數(shù)與管道懸空高度比之間的關系Fig. 12 The relationship of drag coefficient and suspended height ratio

基于圖12中拖曳力系數(shù)數(shù)據(jù)，給出了管道拖曳力系數(shù)與懸跨高度之間擬合關系，兩者擬合公式：

式中：a，b，c是與Renon-Newtonian有關的擬合系數(shù)，二者關系如表3和圖13所示。

表3 不同雷諾數(shù)下a、b、c擬合系數(shù)值Tab. 3 The values of a, b and c under different Reynolds Number

圖13 a、b、c與Renon-Newtonian之間的關系Fig. 13 The relationships of a, b, c and non-Newtonian Reynolds number

4 結 語

1)文中數(shù)值計算中泥流的模擬采用赫巴模型，其計算結果與Zakeri相比更接近于實驗結果，證明了本計算模型的合理性。

2)基于數(shù)值計算，進一步優(yōu)化了無量綱拖曳力系數(shù)CD與非牛頓流體雷諾數(shù)Renon-Newtonian之間的關系，得到了管道拖曳力系數(shù)隨非牛頓流體雷諾數(shù)的增大而減小的規(guī)律。

3)研究不同懸跨高度下管道拖曳力系數(shù)的變化，發(fā)現(xiàn)管道拖曳力系數(shù)隨著H/D的增大而增大。在Renon-Newtonian一定的情況下，存在臨界懸跨高度比(H/D)Critical，在小于此值時，管道拖曳力系數(shù)隨之增大，大于此值時，管道拖曳力系數(shù)保持穩(wěn)定。

4)擬合得到了管道拖曳力系數(shù)與懸跨高度之間的簡單計算公式，基于此公式可計算不同雷諾數(shù)和懸跨高度下的管道拖曳力，為海底管線工程設計提供參考。

[1] 李家鋼，修宗祥，申宏，等. 海底滑坡塊體運動研究綜述[J]. 海岸工程，2012，31(4): 67-78. (LI Jiagang, XIU Zongxiang, SHEN Hong, et al. A review of the studies on submarine mass movement[J]. Coastal Engineering, 2012, 31(4): 67-78. (in Chinese))

[2] VANNESTE M, FORSBERG C F, GLIMSDAL S, et al. Submarine landslides and their consequences: what do we know, what can we do? [M]. Landslide Science and Practice, Springer Berlin Heidelberg, 2013: 5-17.

[3] HUNGR O. A model for the runout analysis of rapid flow slides, debris flows, and avalanches[J]. Canadian Geotechnical Journal, 1995, 32:610-623.

[4] HANCE J J. Development of a database and assessment of seafloor slope stability based on published literature[D]. University of Texas at Austin, 2003.

[5] BOUKPETI N, WHITE D J, Randolph M F. Analytical modelling of the steady flow of a submarine slide and consequent loading on a pipeline[J]. Géotechnique, 2012, 62(2):137-146.

[6] COUSSOT P, PIAU J M. On the behaviour of fine mud suspensions[J]. Rheological Acta, 1994, 33: 175-184.

[7] LOCAT J, LEE H J. Submarine landslides: advances and challenges[J]. Can. Geotech. J., 2002,39(1): 193-212.

[8] JEANJEAN P, HAMPSON K, EVANS T, et al. An operator’s perspective on offshore risk assessment and geotechnical design in geohazard-prone areas[C]//Proceedings of the International Symposium on Frontiers in Offshore Geotechnics (IS-FOG 2005). 2005: 115-144.

[9] 王立忠, 繆成章. 慢速滑動泥流對海底管道的作用力研究[J]. 巖土工程學報, 2008, 30(7): 982-987. (WANG Lizhong, MIAO Chengzhang. Pressure on submarine pipelines under slowly sliding mud flows[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(7): 982-987. (in Chinese))

[10] PARKER E J, TRAVERSO C, MOORE R, et al. Evaluation of landslide impact on deepwater submarine pipelines[C]//Proceedings of Offshore Technology Conference. 2008: OTC19459.

[11] YUAN F, LI L L, GUO Z, et al. Landslide impact on submarine pipelines: analytical and numerical analysis[J]. Journal of Engineering Mechanics, 2014.

[12] ZAKERI A, H?EG K, NADIM F. Submarine debris flow impact on pipelines - PartⅠ: experimental investigation[J]. Coastal Engineering, 2008, 55: 1209-1218.

[13] ZAKERI A, CHI K, HAWLADER B. Centrifuge modeling of glide block and out-runner block impact on submarine pipelines[C]//Proceedings of the Offshore Technology Conference. 2011:OTC 21256.

[14] ZAKERI A, H?EG K, NADIM F. Submarine debris flow impact on pipelines—Part II: numerical analysis[J]. Coastal Engineering, 2009, 56(1):1-10.

[15] AKERI A, HAWLADER B, CHI K. Drag forces caused by submarine glide block or out-runner block impact on suspended (free-span) pipelines[J]. Ocean Engineering ,2012, 47: 50-57.

[16] SAHDI F, GAUDIN C, WHITE D J, et al. Centrifuge modeling of active slide-pipeline loading in soft clay[J]. Géotechnique, 2014, 64(1): 16-27.

[17] RANDOLPH M F, DAVID J W. Interaction forces between pipelines and submarine slides - A geotechnical viewpoint[J]. Ocean Engineering, 2012, 48:32-37.

[18] WHITE F M. Viscous fluid flow[M]. McGraw-Hill Higher Education, Boston (XXI), 2006:629.

[19] ZAKERI A. Submarine debris flow impact on suspended (free-span) pipelines: normal and longitudinal drag forces[J]. Ocean Eng., 2009, 36(10):489-499.

[20] HUANG X, GARCIA M H. A herschel - bulkley model for mud flow down a slope[J]. Journal of Fluid Mechanics, 1998, (374): 305-333.

[21] IMRAN J, HARFF P, PARKER G. A numerical model of submarine debris flow with graphical user interface[J]. Computers & Geosciences, 2001，27: 717-729.

[22] SI Guangying. Experimentally study the rheology of fine-grained slurries and some numerical simulations of downslopes slurry movements[D]. University of OSLO, 2007.

[23] WRIGHT J, COLLING A, PARK D. Waves, tides, and shallow-water processes[M]. Gulf Professional Publishing, 1999.

Drag force of submarine landslides mudflow impacting on a suspended pipeline

LI Hongwei1,2, WANG Lizhong1, GUO Zhen1, YUAN Feng1

(1. Research Center of Coastal and Urban Geotechnical Engineering, Zhejiang University, Hangzhou 310058, China; 2. College of Hydraulic and Civil Engineering, Tarim University; Aler 843300, China)

Submarine landslide is one of the most dangerous geological disasters for the subsea installation of petroleum production. Under the action of the seawater, submarine landslide will change to slurry and impact on a submarine suspended pipeline. Computational fluid dynamics (CFD) is one of numerical methods employed to analyze the landslides impacting on the suspended pipeline. Herschel-Bulkley model is adopted, which is more suitable for modeling viscous slurry, to simulate slurry impacting on submarine pipeline. The effects of the suspending height from seabed are also discussed. The study finds that normal drag coefficient increases with the rising of the suspending height. When a critical suspending height is reached, the drag coefficient is constant.

submarine landslide; computational fluid dynamics (CFD); slurry flow; drag force; non-newtonian fluids; submarine pipeline

P756.2

A

10.16483/j.issn.1005-9865.2015.06.002

1005-9865(2015)06-010-10

2015-01-26

國家杰出青年基金(51325901);國家自然科學基金(51279176, 51409228, 51209183);博士后科學基金特別資助(2014T70574)

李宏偉(1985-)，男，甘肅隴西人，助教，從事巖土工程方面的研究。E-mail:zjulhw@163.com