李宏伟，王立忠，国 振，袁 峰

(1. 浙江大学 滨海和城市岩土工程研究中心,浙江 杭州 310058；2. 塔里木大学 水利与建筑工程学院， 新疆 阿拉尔 843300)

海底泥流冲击悬跨管道拖曳力系数分析

李宏伟1,2，王立忠1，国 振1，袁 峰1

(1. 浙江大学 滨海和城市岩土工程研究中心,浙江 杭州 310058；2. 塔里木大学 水利与建筑工程学院， 新疆 阿拉尔 843300)

海底滑坡是海洋油气工程最危险的地质灾害之一，直接影响海底管线运营安全。滑坡体失稳滑动过程中，由于海水掺和作用逐渐加速转变成快速滑动泥流，冲击海底悬跨管道。基于当前国际通用的计算流体动力学(CFD)方法，采用赫巴模型描述快速滑动泥流，计算分析海底滑坡冲击悬跨管道的受力特性，重点分析悬跨高度对管道法向拖曳力系数的影响。研究发现，管道法向拖曳力系数随悬跨高度的增长而增大，当达到某一悬跨高度时，管道拖曳力系数保持稳定。

海底滑坡；计算流体动力学(CFD)；泥流；拖曳力；非牛顿流体；海底管线

海底管道具有输送连续、管理方便、效率高、成本低等优点，是深水油气输运的重要生命线，其直径一般在0.1～1.0 m之间。由于海底管道长距离输运，管道经历的海底地形崎岖，起伏变化大，容易形成管道局部悬跨段。而这些悬跨段所在的海沟、海槽等低洼区域往往是海底滑坡泥流的主要输运通道。海底滑坡泥流的流速极快，可达30 m/s，其影响面积也很广。因此，对于悬跨段管道而言，极易遭受滑坡泥石流的快速冲击作用而破坏[1]，有必要对悬跨段管道遭受快速泥流冲击时的受力特点展开研究。

Vanneste等[2]指出，目前对海底滑坡方面的研究主要针对滑坡形成的浊流阶段，将滑坡土体简化成黏塑性流体，通过水动力学方法模拟滑坡土体的长距离滑移[3-5]。Coussot等[6]讨论了泥浆流温度、PH值、电解质浓度、固体颗粒含量和黏土含量等参数对泥流流变性能的影响，指出赫巴模型中屈服应力能够很好估计泥流真实的屈服应力状态。Locat等[7]和Jeanjean等[8]对海底滑坡及其影响区域内的海底管道安全性进行了论述。海底滑坡泥流对管道的作用可分为快速泥流的冲击作用和慢速泥流的推力作用，已有研究大多针对慢速泥流。王立忠等[9]分析了慢速泥流作用下管道受力的一般规律。Parker等[10]、Yuan等[11]直接将海底滑坡冲击荷载简化为定值，分析了海底滑坡对管道整体受力的影响，简化了海底滑坡快速泥流的冲击作用。Zakeri等[12-15]开展了海底滑坡泥流冲击海底管道的模型实验研究和数值模拟计算，得出了管道拖曳力系数的经验公式。Randolph等[16-17]开展了海底滑坡冲击海底管线的离心机模型试验研究，提出了基于岩土力学和流体力学方法的管道拖曳力计算方法。

基于计算流体动力学(CFD)方法对海底滑坡泥流正面冲击悬跨管道的过程进行数值模拟。在计算中，采用赫巴模型描述滑坡泥流，着重分析了滑坡泥流冲击不同悬跨高度下管道的拖曳力系数，得出了拖曳力系数与悬跨高度(H)之间的关系，提出了更为合理的拖曳力系数取值方法，为海底管线设计提供参考依据。

1 滑坡冲击荷载计算方法

海底滑坡泥流冲击管道主要体现在拖曳力的大小，而拖曳力系数是衡量拖曳力大小的关键参数。已有研究主要基于水动力学方法，采用Morison公式计算细长结构物上流体作用力，这是一种半经验半理论的描述。White[18]提出了法向拖曳力系数(CD)的计算方法：

式中：FD是法向拖曳力，U是滑体上游自由流体流速，ρ是滑坡体流体的密度，A为障碍物的特征面积。

Zakeri等[14，19]采用幂律模型通过数值(CFD)模拟手段研究了海底滑坡对海底管线的冲击荷载，对室内试验结果进行了反演分析，提出了海底滑坡泥流正面冲击管道时产生的法向拖曳力系数如式(2)、(3)。

管底距离海床1D处的管道：

管底距离海床2 mm处的管道：

式中：Renon-Newtonian是非牛顿流体雷诺数，如式(4)所示。

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

2 数值模型

2.1水和泥流控制方程

采用ANSYS CFX不可压缩两相流(水和泥流)进行模拟计算，不同相采用小写希腊字母α和β表示，总相数为Np=2。其控制方程如下：

1)连续性方程

2)动量方程

式中：Pα和μα是α相的压强和黏性系数；SMα是外部质量力引起的动量源相；Mα是由于其他相引起的作用在α相上的总界面力，由下式可得

包括拖曳力、升力、壁面润湿力、虚拟质量力和湍流耗散力等[14]。

2.2泥流流变本构

目前，研究滑坡泥流流变本构的模型有宾汉模型、赫巴模型和幂律模型，各模型的剪切应力与应变率之间的关系如式(8)、(9)所示。

赫巴模型(herschel-bulkley model):

幂律模型(power-law model):

2.3模型设定

2.3.1 计算域设定

Zakeri等[12]进行了泥流冲击管道的常重力模型试验，其中管道直径采用22.2 mm和28.6 mm两组类型。由于采用无量纲参数描述拖曳力系数，因而这里统一采用计算模型管道直径为25 mm。图2所示为一个典型的计算模型，在计算域中，管道距离入口边界为9.5D，距离出口边界为13.5D，距离海床边界为4D，距离计算域上边界为5D，水与泥流交界面距离计算域上边界为1D，泥流从左边入口边界进入计算域，从右边出口边界流出，海底管道悬跨高度为H，H为管道中心到底部边界的距离。

图1 流变模型Fig. 1 Rheological model

图2 几何模型和数值模型Fig. 2 Domain geometry and simulation model

采用图2模型进行数值模拟，将滑坡泥流正面冲击悬跨管道模拟成两相(水和泥流)并考虑浮力连续相的自由液面流。滑坡泥流采用考虑浮力相的层流模型，其本构关系采用赫巴模型；水中产生的湍流采用k-ε模型。在两相流中，由于水和泥流的密度不同，从而影响流体计算的控制方程(式(5)、式(6))。浮力相SM,Buoyancy,α=(ρα-ρref)g引入动量方程式(6)中，其中g是重力加速度，ρα是滑坡泥流的密度，ρref是浮力参考密度，为水的密度997kg/m3。由于两相(水和泥流)的速度场不同，因而两相自由液面间存在界面拖曳力，其通过无量纲参数界面拖曳力系数来描述界面拖曳力，其界面拖曳力系数为2～3[14]。

2.3.2 边界设定

计算域的边界如图3所示，管道表面为无滑移粗糙壁面边界，其粗糙系数ks=0.001 5 mm[14]，出口边界为开放边界，计算域前后边界为对称边界，其中，对称边界认为是无剪切应力的平面应变边界，计算域上部边界为自由滑移壁面边界，下部边界为无滑移粗糙壁面边界，其粗糙系数ks=0.5 mm[14]。

图3 流体计算域边界Fig. 3 Fluid domain boundary conditions

2.3.3 网格剖分

采用ANSYS ICEM软件对整个计算域进行非结构化网格分区划分。如图4所示，对管道周围网格进行局部细化，这部分区域是本文研究泥流流动状态及流体与管道相互作用的重要区域。计算域整体最大网格尺寸为25 mm，管道表面网格大小为2 mm逐渐增长到50 mm处的7 mm，管道表面有覆面层网格，需对其分五层，覆面层网格总厚度为2 mm[19]。覆面层内为三棱柱网格，以减小管道表面的网格扭曲率。

图4 管道直径为25 mm的计算域网格与覆面层网格Fig. 4 The domain grid size and inflated boundary, with 25 mm pipe diameter

2.3.4 求解设定

在数值模拟中，初始时刻(t=0)需对计算域进行初始化，计算域中水的体积分数为1，各速度分量为0。时间步的选择由式(10)确定：

φip=φup+β

3 数值分析算例

SI[22]的实验研究了由高岭土、石英砂和水混合配置而成的泥流的一般特性，给出了幂律模型与赫巴模型的流变参数，如表1所示。利用表1中六种不同成分的泥浆流体的流变参数进行模拟计算。

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

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

3.1验证与对比

Zakeri[19]采用表1中的幂律模型，研究了管底距离水槽底部1D、2 mm时的管道拖曳力系数变化。采用表1中的赫巴模型计算的CD值与Zakeri的实验及计算结果进行对比。从图5可以看出，采用赫巴模型计算结果与实验结果更为接近，尤其是在非牛顿流体雷诺数小于10的情况。非牛顿流体雷诺数增大，法向拖曳力系数减小。拟合公式如式(12)、(13)，两者拟合R2值分别为99.00%和98.40%。

管底距离海床1D处的管道：

管底距离海床2 mm处的管道：

通过式(12)与(13)可以计算得到每延米管道上产生的法向拖曳力。

由图5观察到，当Renon-Newtonian<10时，本文与Zakeri计算结果都高于实验值。这主要由于在试验中，管道表面粗糙，凸凹不平，局部存在光滑区域，泥流流过管道时，两者接触面之间存在少量的水泡，使得测出的管道拖曳力较小。而在数值模拟中，管道表面采用无滑移粗糙表面，在低雷诺数下，黏性增大，管道上下表面处泥流流速大，管道的拖曳力较大。

图5 结果比对Fig. 5 Comparing with the results of Zakeri

3.2参数分析

3.2.1 层流边界层

首先对不放置管道近海床区域泥流流动的边界层发展情况进行了模拟计算。采用入口流速U=1.0 m/s，泥流为35%Clay，海床粗糙系数为ks=0.5 mm[14]的情况下，得到了泥流计算域速度矢量图，如图6所示。靠近海床流速为零，取计算域中流速达到0.99U∞的高度作为边界层的充分发展高度。

图6 边界层流速矢量图Fig. 6 The velocity vector of boundary layer

由于近海床区域泥流受到黏性影响，毗邻海床附近的泥流停滞下来，速度为零，而在垂直于流动方向上形成速度梯度，并使靠近海床的泥流速度减慢，随着黏滞泥流向前移动，边界层逐渐增加，更多黏性泥流速度减缓，逐渐形成一稳定的层流边界层。计算域入口流速、泥流密度及海床粗糙度对层流边界层充分发展高度有影响，如图7所示。其中，海床粗糙度主要取决于海床起伏(比如砂脊、砂波的影响)和沉积物粒径大小，这里未考虑海床起伏的影响，只研究了粒径大小在粉土(ks=0.001 mm)到粗砂(ks=1 mm)范围内变化对层流边界层的影响[23]。

图7 边界层厚度δFig. 7 Boundary layer thickness

从图7可以看出，Clay含量越高，泥流的黏度越大，层流边界层充分发展的高度越大；泥流速度越大，层流边界层充分发展的高度越小；层流边界层充分发展高度不受海床粗糙系数的影响。从图7(a)、(b)可以看出，计算域入口流速为1.0 m/s、海床边界粗糙系数为ks=0.5 mm[14]的情况下，不同Clay含量(10%Clay、15%Clay、20%Clay、25%Clay、30%Clay、35%Clay)的泥流，随着Clay含量的增高，层流边界层厚度也逐渐发展，其厚度在0.4D～1.4D之间变化；图7(c)、(d)揭示了在35%Clay含量的泥流中，在海床边界粗糙系数为ks=0.5 mm[14]的情况下，随着计算域入口流速的增大，层流边界层厚度逐渐减小，其厚度在0.4D～1.8D之间变化；从图7(e)、(f)可以看出，在计算域入口流速为1.0 m/s、35%Clay含量的泥流情况下，边界层随海床粗糙系数的增大而保持在1.24D。

图8 拖曳力系数与雷诺数之间的关系曲线Fig. 8 Drag coefficient vs. Reynolds number curves

3.2.2 悬跨高度比H/D

考虑层流边界层的影响，对不同悬跨高度(H)的海底管道进行了数值模拟分析，得到了不同悬跨高度下的管道拖曳力系数，如图8所示。从图8中可以看出，管道拖曳力系数随非牛顿流体雷诺数的增大而减小，随管道悬跨高度的增高而增大。当非牛顿流体雷诺数较大时，层流边界层较小，H/D从0开始增长很小范围后CD才保持稳定；当非牛顿流体雷诺数较小时，层流边界层较大，H/D在从0开始增长很大范围后CD才保持不变。为了方便比较，将H/D=5时的拖曳力系数CD作为所有工况的C。

图9 不同雷诺数下的管周围流态示意(H/D=5)Fig. 9 The flow regime around the pipe under different Reynolds number at H/D=5

图9揭示了海底管道悬跨高度为H=5D的情况下，不同非牛顿流体雷诺数下管道周围流态示意图。从图9中可以看出，管道周围流态稳定，并未受到层流边界层的影响。当Renon-Newtonian<10时，管道侧面的泥流流线不存在分离点；当Renon-Newtonian>10时，管道侧面的泥流流线分离点逐渐向泥流上游移动，管道后面区域出现回流现象。

图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可知，管道拖曳力系数随着管道悬跨高度的增加而增大，当悬跨高度比(H/D)增加到2之后，拖曳力系数保持稳定。这里将此悬跨高度比定义为影响拖曳力系数的临界悬跨高度比(H/D)Critical。由于滑坡泥流正面冲击管道时管道周边只形成很薄的边界层，影响管道拖曳力系数的主要为近海床层流边界层的厚度。因此，临界悬跨高度比可通过下式简单计算：

表2给出了不同雷诺数下层流边界层的厚度变化。对比图12，可以发现，式(14)给出的临界悬跨高度比较好地预测了管道拖曳力系数的变化。

表2 不同雷诺数下边界层厚度δTab. 2 Boundary layer thickness under different Reynolds Number

图11 管道拖曳力系数与悬跨高度的关系(Renon-Newtonian=21.3)Fig. 11 The relationship of the drag coefficient and the suspended height at Renon-Newtonian=21.3

图12 拖曳力系数与管道悬空高度比之间的关系Fig. 12 The relationship of drag coefficient and suspended height ratio

基于图12中拖曳力系数数据，给出了管道拖曳力系数与悬跨高度之间拟合关系，两者拟合公式：

式中：a，b，c是与Renon-Newtonian有关的拟合系数，二者关系如表3和图13所示。

表3 不同雷诺数下a、b、c拟合系数值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)文中数值计算中泥流的模拟采用赫巴模型，其计算结果与Zakeri相比更接近于实验结果，证明了本计算模型的合理性。

2)基于数值计算，进一步优化了无量纲拖曳力系数CD与非牛顿流体雷诺数Renon-Newtonian之间的关系，得到了管道拖曳力系数随非牛顿流体雷诺数的增大而减小的规律。

3)研究不同悬跨高度下管道拖曳力系数的变化，发现管道拖曳力系数随着H/D的增大而增大。在Renon-Newtonian一定的情况下，存在临界悬跨高度比(H/D)Critical，在小于此值时，管道拖曳力系数随之增大，大于此值时，管道拖曳力系数保持稳定。

4)拟合得到了管道拖曳力系数与悬跨高度之间的简单计算公式，基于此公式可计算不同雷诺数和悬跨高度下的管道拖曳力，为海底管线工程设计提供参考。

[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