中缅油气管道龙陵段管沟回填土冲蚀过程模拟

葛华　方迎潮　赵飞　孔志岗　戴国文　张垚

摘要：基于Navier-Stokes方程和牛顿第二定律构建的颗粒离散元流-固耦合数值模拟方法，以滇西龙陵县中缅油气管道穿越的强烈风化花岗岩区的土体为研究对象，构建了冲沟数值模型，研究了边坡回填土的冲蚀过程。模拟结果显示：① 当模型中水流速度较大时，土坎被迅速冲毁，冲毁的土体呈现散体状，并且沟底的部分土体颗粒也被流水挟带走;随着水流速度降低，土坎被冲毁的速度降低，冲毁的土体局部黏结在一起，流水对沟底的侵蚀现象不明显。② 提高模型中颗粒间的黏结强度，在较低流速的流水冲刷下土坎未出现冲蚀破坏现象;较高流速的流水作用下，被冲走的土坎呈整体块状，沟底未出现侵蚀的现象。③ 降低坡面流水的动力条件和改善土体本身的强度特征将有助于防止土体流水侵蚀现象的发生。

关键词：油气管道边坡; 回填土冲蚀; 水流速度; 离散元方法; 流-固耦合方法

中图法分类号： P66

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2021.11.020

0引 言

中缅油气管道呈线状自西向东分别穿越滇西横断山脉、滇中红色高原两大地貌格局，沿线属于“三高四活跃”的不良地质作用发育地区（高地震烈度、高地应力、高地热，活跃的新构造运动、活跃的地热水环境、活跃的外动力地质条件、活跃的岸坡再造过程），地质复杂，立体气候明显，单点暴雨活动强烈。如此背景之下，油气管道作业边坡的开挖等人类工程活动破坏了原始地貌、地质和生态平衡，开挖回填后，回填土植被恢复差，造成管道铺设区域坡面滑塌、水土流失等一系列环境地质问题。有关水土流失及土体浅表层冲刷的研究，前人针对黄土地区、红土地区、丘陵地带的水土流失做了大量的研究[1-3]。以往研究表明，边坡土层冲刷主要受控于水动力条件和土体本身的物理力学性质。降雨强度和边坡坡度影响着水动力强度，降雨强度、坡度越大，边坡侵蚀越强烈。另外，土体级配越好，边坡越容易发生侵蚀;经过改良后的土体或有植被防护的土体通常具有较强的抗冲刷能力[4-5]。冯秀等[6]通过室内人工模拟降雨试验，研究坡度、降雨强度和地表覆盖等因素对坡面侵蚀过程的影响，认为坡面径流率随降雨强度和坡度增大而增大;随着降雨强度的增大，产沙率明显增大，且降雨强度越大，坡度对产沙率的影响越明显。杜婷婷等[7]认为在黄土中增加改良剂能够提高边坡的抗冲刷的能力。裴向军等[8]研究发现，经过改性钠羧甲基纤维素胶结固化后的土质边坡，能有效防止坡面冲沟发育，减缓边坡水土流失，增强其抗冲刷性能力。

近年来，基于数值模拟探索土体冲蚀形成机理成为水土流失研究新的热点。由于颗粒离散元方法在描述固体介质材料破坏、运动等方面的优势，被广泛应用于岩土工程、地质灾害成因研究中[9-14]。土体冲蚀过程是土体颗粒在水流作用下破坏并被挟带运移的过程，是非常典型的流固耦合作用过程，基于颗粒离散元的流固耦合方法将能够揭示这一过程中的机理。吴谦等[15]利用三维颗粒流软件结合室内降雨模型试验分析了倾角为70°黄土邊坡的冲刷过程，基于颗粒运动路径、土体孔隙率变化等参数揭示了降雨对边坡侵蚀能力的分布规律。张雁等[16]，柯云斌等[17]基于二维颗粒流方法分别探索了黄土和残坡积土边坡的冲刷破坏机理。离散元数值模拟在水土流失研究方面有了初步成果，但是从细观尺度揭示土颗粒与流水相互作用的过程的研究尚少，尤其是尚未开展流水直接冲刷坡面的三维模拟研究。中缅油气管道龙陵段管道作业边坡回填土冲蚀过程其实就是流水夹带土壤颗粒形成的颗粒流流失过程，本文基于颗粒离散元构建的流固耦合算法，采用数值模拟方法，从细观尺度研究水动力条件与土体力学性质对土体冲蚀破坏过程的影响，分析流水侵蚀的起动条件，为该类型土体坡面水土流失治理提供理论依据。

1边坡回填土物理力学性质及冲蚀特征

研究区位于滇西龙陵县，地处怒江、龙川江两江之间，滇西纵谷南段，高黎贡山南延部分，地势呈中部高而东西斜势，海拔高程一般在1 600 m左右，属低-中高山区。年平均降水量2 105.7 mm，月最大降水量726.6 mm，日最大降水量134.7 mm，最大单点降雨量56.9 mm。雨季主要集中于5～10月，降水量占全年降水量的82.6%。

研究区范围内地层岩性比较单一，油气管道沿线出露主要地层由新至老主要有冲洪积层（Q4al+pl）、残坡积黏性土（Q4el+dl）;中生界（T、J、K）地层，为一套紫红色、杂色的泥岩、长石石英砂岩、粉砂岩、夹煤层。燕山期花岗岩、混合岩、粗面岩在研究区附近大面积出露。管道穿越段主要为花岗岩分布区。坡面表层为网纹状全风化花岗岩层（粗砂土层），全风化花岗岩中矿物成分最高为石英，其次是长石、高岭土等，原岩中的长石、云母已经分解，石英因抗风化能力强而保留下来。全风化花岗岩物理性质与砂土类似，但颗粒粒度更大，孔隙比大，裂隙发育，具有一定崩解性。管道作业边坡回填土厚度一般为3 m，局部地段有变化。回填土主要以原地的全风化花岗岩回填，以砂砾土为主，土体中夹杂着大量的花岗岩砾石和石英砾石。砾石粒径介于2.0～20.0 cm，含量3%～8%;砂土的颗粒粒径一般为1～2 mm，2.0 cm>d>0.075 mm粒径的颗粒含量65.6%～80.2%，小于0.075 mm颗粒含量16.6%～26.2%。通过对研究区边坡回填土取样并进行室内试验，得到土体的物理力学参数，如表1所列。

输油气管道的铺设过程中，需开挖15～25 m宽的作业边坡，管沟铺设管道后进行回填，回填土为早期开挖的土或就地取材。管沟一般深2～3 m，因作业边坡的开挖和回填，不仅直接剥离了抗冲蚀性较强的本来就比较薄的土壤层，使砂土层直接出露于地表，而且破坏了原有的周边环境和生态。这种天然的原生岩土体条件和生态环境一旦破坏，不容易恢复。同时，回填土由于各种粒径颗粒混杂，固结性差，颗粒间黏结性差，容易被冲刷，使冲蚀从天然状态下处于相对平衡状态的较为缓慢的自然侵蚀过程，在极短的时间迅速转变为强烈的工程性冲蚀过程。根据现场调查，龙陵油气管道穿越区坡面冲蚀类型主要包括雨滴溅蚀、面蚀、细沟冲蚀、浅沟冲蚀、冲沟冲蚀等（见图1）。

2基于颗粒流的流-固耦合算法

2.1颗粒流方法

Cundall提出了可以分析颗粒材料力学性态的离散单元法，从而可以对颗粒团粒体的稳定、变形及本构关系进行分析，对固体力学大变形问题具有良好的适应性[18]。它是将颗粒介质的运动及其相互作用通过圆形（或异型）离散单元来模拟。每一时刻颗粒在平面内的位置和速度由平动和转动运动方程来确定。

PFC3D（Particle Flow Code 3 Dimension）即三维颗粒流程序，是以球形颗粒介质的运动及其相互作用为研究对象的一种离散单元方法。该方法通过将实际的材料离散成由若干个刚性颗粒组成的模型，来研究实际材料的各种力学特性。该程序最初被用来研究颗粒介质特性，它利用有代表性的大量颗粒单元来表示实际物体，并利用这种局部的模拟结果来研究连续介质边值问题。PFC 3D颗粒流程序基于细观力学通过离散单元方法来模拟球形颗粒介质的运动及其相互作用，以牛顿第二定律和力–位移定律为基础，采用显式差分算法，循环应用牛顿定律分析颗粒的运动特征、力–位移定律分析颗粒的接触特征，从而实现对颗粒介质的运动及其相互作用过程的模拟[11-12，14]，控制方程如公式（1）所示。

3冲蚀过程数值模拟

3.1数值模型与参数

边坡坡面的冲蚀是一个复杂的过程，涉及到降雨、土体饱和、雨水汇集等连续的环节。在降雨过程中，坡面土体通常会经历雨滴溅蚀、片蚀、细沟侵蚀、切沟侵蚀等阶段。这些阶段的产生与降雨强度、坡面形态及坡面土体本身的物理力学性质相关。降雨强度及坡面形态主要影响坡面侵蚀中的水动力条件。本次模拟研究中仅考虑水力条件和土体力学性质对坡面土体冲刷过程的影响，不讨论水动力条件成因。构建一个100 mm×100 mm×100 mm的流固模型（见图2），包括流体计算网格和固体颗粒。模型中流体计算网格规模为10×10×10，每个网格单元大小为10 mm×10 mm×10 mm;土體整体厚度为50 mm，右侧有个高为20 mm土坎。模型中土体由球颗粒构成，其粒径介于2.0～3.2 mm，颗粒数量为26 850个，颗粒之间采用黏结模型进行连接。模型中的参数如表2所列。

3.2水流冲蚀模拟结果

模拟过程中，固体部分首先在重力作用下达到力学平衡状态，形成初始的土体物质基础;基于模型试验及前人经验，设定模型左侧入水口初始水流速度分别为2.0，1.0 m/s和0.5 m/s，模型右侧排水口流体压力为0 Pa，形成初始的稳定流场;之后，开启流固耦合计算。图3～5给出了3种水流速度作用下模型中土体颗粒运移距离、流体速度分布与流体压力分布随时间的演化过程。由图3～5可见，凸起的土坎在3种流速的水流侵蚀作用下逐渐被冲毁。当入口水流速度较大时（流速为2.0 m/s），在较大的流体拖拽力作用下土坎被迅速冲毁，呈现散开状被挟带走，并且沟底的部分土体颗粒也被流水刮起;随入口水流速度降低，土坎被冲毁的速度降低，冲毁的土体局部还黏结在一起，流水对沟底的侵蚀现象不明显。由此可见，水动力条件是影响土体冲蚀过程的关键参数之一，降低水动力条件是防止边坡土体水毁的重要手段之一。

除水流速度之外，土体冲蚀破坏与其自身的强度特征关系密切。在模拟时，通过提高模型中颗粒间的黏结强度来改善土体的力学性能，进而分析土体力学特征对冲蚀过程的影响。图6展示了土体黏结强度提高10倍之后的模拟结果。与图3～5相比发现，当入口流

速为2.0 m/s时，土体黏结强度增大后土坎整体被流水冲走，但是未出现沟底侵蚀的现象;当入口流速降为0.5 m/s时，土坎未产生冲蚀破坏现象。由此可见，改善土体本身的强度特征有助于防止土体水毁现象的产生。

上述模拟从细观尺度探讨水动力条件与土体力学性质对土坎冲蚀破坏过程的影响。当模拟水流速度大时，土坎以散开形式被迅速冲毁，沟底的部分土体颗粒也被流水刮擦带走;降低水流速度后，土坎以团块的形式被冲走，并没有被冲毁。另外，提高模型中颗粒间的黏结强度能够改善土体的力学性能，进而提高土体的抗冲蚀的能力。

4结 语

本文基于Navier-Stokes方程和牛顿第二定律构建的颗粒离散元流固耦合数值模拟方法，以强风化花岗岩区的土体为研究对象构建了冲沟数值模型，研究了边坡回填土的冲蚀过程。研究表明本文中构建的数值模拟方法能够较好地再现土体冲蚀过程。研究还发现降低坡面流水的动力条件和改善土体本身的强度特征将有助于防止土体流水侵蚀现象的发生，对于边坡水土流失治理具有一定的启示意义。本文关于土体冲蚀的研究获得了一些定性规律的认识，下一步将结合具体降雨特征、汇水区面积、边坡的土体特征、坡长和坡高等几何特征，获得边坡土体冲蚀的启动流速等定量数据，以指导边坡水土流失的防治工作。

参考文献：

[1]李宗善，杨磊，王国梁，等.黄土高原水土流失治理现状、问题及对策[J].生态学报，2019，39（20）：7398-7409.

[2]白皓.红土工程堆积体坡面水蚀过程及泥沙颗粒分布特性研究[D].杨凌：西北农林科技大学，2018.

[3]王玉宽，文安邦，张信宝.长江上游重点水土流失区坡耕地土壤侵蚀的～（137）Cs法研究[J].水土保持学报，2003，17（2）：77-80.

[4]王玲.陡坡地水蚀过程与泥沙搬运机制[D].北京：中国科学院大学，2016.

[5]汤珊珊.覆沙坡面水蚀产沙动力过程与微地貌变化的响应研究[D].西安：西安理工大学，2018.

[6]冯秀，查轩，黄少燕.人工模拟降雨条件下花岗岩红壤坡面侵蚀过程与特征分析[J].中国水土保持科学，2014，12（1）：19-23.

[7]杜婷婷，李志清，王晓明，等.黄土边坡降雨冲刷模型试验研究[J].工程地质学报，2018，26（3）：732-740.

[8]裴向军，杨晴雯，许强，等.改性钠羧甲基纤维素胶结固化土质边坡机制与抗冲蚀特性研究[J].岩石力学与工程学报，2016，35（11）：2316-2327.

[9]王学良，孙娟娟，周书，等.尾矿库溃坝运动特征模拟研究[J].工程地质学报，2019，27（1）：144-151.

[10]眭靜，姜元俊，樊晓一，等.碎屑流冲击刚性挡墙的力学模型研究[J].岩石力学与工程学报，2019，38（1）：121-132.

[11]费建波，介玉新，张丙印，等.土的三维破坏准则在颗粒流模型中的应用[J].岩土力学，2016，37（6）：1809-1817.

[12]朱涵成，韩文喜，陈超.砂岩常规三轴的颗粒流数值模拟[J].地质灾害与环境保护，2013，24（3）：104-108.

[13]徐金明，谢芝蕾，贾海涛.石灰岩细观力学特性的颗粒流模拟[J].岩土力学，2010，31（增2）：390-395.

[14]杨庆华，姚令侃，杨明.地震作用下松散堆积体崩塌的颗粒流数值模拟[J].西南交通大学学报，2009，44（4）：580-584.

[15]吴谦，王常明，宋朋燃，等.黄土陡坡降雨冲刷试验及其三维颗粒流流-固耦合模拟[J].岩土力学，2014，35（4）：977-985.

[16]张雁，高树增，闫超群，等.降雨对黄土路基边坡的冲刷规律[J].中国地质灾害与防治学报，2017，28（4）：34-39.

[17]柯云斌，符元帅，阙云，等.降雨诱发残积土陡坡坡面冲刷破坏颗粒流模拟[J].济南大学学报（自然科学版），2019，33（5）：453-459.

[18]CUNDALL P A.A computer model for simulating progressive large scale movements in Blocky Rock Systems[C]∥Proceedings of the Symposium of the International Society for Rock Mechanics，Society for Rock Mechanics（ISRM），France，1971，II-8.

[19]李识博，王常明，王钢城，等.松散堆积物坝基渗透淤堵试验及颗粒流模拟[J].水利学报，2012，43（10）：1163-1170.

[20]SADEK M A，CHEN Y.Feasibility of using PFC3D to simulate soil flow resulting from a simple soil-engaging tool[J].Transactions of the Asabe，2015：987-996.

[21]宿崇，侯俊铭，朱立达，等.基于流固耦合算法的单颗磨粒切削仿真研究[J].系统仿真学报，2008，20（19）：5250-5253.

[22]方志，陈育民，何森凯.基于单相流的减饱和砂土流固耦合改进算法[J].岩土力学，2018，39（5）：1851-1857.

（编辑：黄文晋）

Abstract：According to the Navier-Stokes equation and Newtons second law，a fluid-mechanical coupling numerical simulation method was constructed based on the particle discrete element method.The soil in the strongly weathered granite area where the Sino-Burmese Oil and Gas Pipeline crossed in Longling County，western Yunnan Province，was employed as our research object，and a numerical model of a gully with a ridge at the bottom was constructed.The results showed that： ① when the fluid velocity in the model was large，the soil ridges were quickly washed away，and the washed-out soil was loose.Some part of the soil particles at the bottom of the gully was also carried away.As the fluid speed decreased，the washed-out soil was partially bonded together，and the erosion of the gully bottom was not apparent.② When the bonded strength between the particles increased in the model，the erosion of ridges did not reproduce at a lower flow rate.When the ridges were washed away under a higher fluid flow rate as a whole block，no erosion occurred at the bottom of the gully.③ The simulation results indicate that reducing the fluid speed on the slope and improving the soils strength will help prevent soil erosion.

Key words：oil and gas pipeline slope;backfill soil erosion;water flow speed;discrete element method;flow-solid coupling method