中緬油氣管道龍陵段管溝回填土沖蝕過程模擬

葛華　方迎潮　趙飛　孔志崗　戴國文　張垚

摘要：基于Navier-Stokes方程和牛頓第二定律構建的顆粒離散元流-固耦合數(shù)值模擬方法，以滇西龍陵縣中緬油氣管道穿越的強烈風化花崗巖區(qū)的土體為研究對象，構建了沖溝數(shù)值模型，研究了邊坡回填土的沖蝕過程。模擬結果顯示：① 當模型中水流速度較大時，土坎被迅速沖毀，沖毀的土體呈現(xiàn)散體狀，并且溝底的部分土體顆粒也被流水挾帶走;隨著水流速度降低，土坎被沖毀的速度降低，沖毀的土體局部黏結在一起，流水對溝底的侵蝕現(xiàn)象不明顯。② 提高模型中顆粒間的黏結強度，在較低流速的流水沖刷下土坎未出現(xiàn)沖蝕破壞現(xiàn)象;較高流速的流水作用下，被沖走的土坎呈整體塊狀，溝底未出現(xiàn)侵蝕的現(xiàn)象。③ 降低坡面流水的動力條件和改善土體本身的強度特征將有助于防止土體流水侵蝕現(xiàn)象的發(fā)生。

關鍵詞：油氣管道邊坡; 回填土沖蝕; 水流速度; 離散元方法; 流-固耦合方法

中圖法分類號： P66

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2021.11.020

0引 言

中緬油氣管道呈線狀自西向東分別穿越滇西橫斷山脈、滇中紅色高原兩大地貌格局，沿線屬于“三高四活躍”的不良地質(zhì)作用發(fā)育地區(qū)（高地震烈度、高地應力、高地熱，活躍的新構造運動、活躍的地熱水環(huán)境、活躍的外動力地質(zhì)條件、活躍的岸坡再造過程），地質(zhì)復雜，立體氣候明顯，單點暴雨活動強烈。如此背景之下，油氣管道作業(yè)邊坡的開挖等人類工程活動破壞了原始地貌、地質(zhì)和生態(tài)平衡，開挖回填后，回填土植被恢復差，造成管道鋪設區(qū)域坡面滑塌、水土流失等一系列環(huán)境地質(zhì)問題。有關水土流失及土體淺表層沖刷的研究，前人針對黃土地區(qū)、紅土地區(qū)、丘陵地帶的水土流失做了大量的研究[1-3]。以往研究表明，邊坡土層沖刷主要受控于水動力條件和土體本身的物理力學性質(zhì)。降雨強度和邊坡坡度影響著水動力強度，降雨強度、坡度越大，邊坡侵蝕越強烈。另外，土體級配越好，邊坡越容易發(fā)生侵蝕;經(jīng)過改良后的土體或有植被防護的土體通常具有較強的抗沖刷能力[4-5]。馮秀等[6]通過室內(nèi)人工模擬降雨試驗，研究坡度、降雨強度和地表覆蓋等因素對坡面侵蝕過程的影響，認為坡面徑流率隨降雨強度和坡度增大而增大;隨著降雨強度的增大，產(chǎn)沙率明顯增大，且降雨強度越大，坡度對產(chǎn)沙率的影響越明顯。杜婷婷等[7]認為在黃土中增加改良劑能夠提高邊坡的抗沖刷的能力。裴向軍等[8]研究發(fā)現(xiàn)，經(jīng)過改性鈉羧甲基纖維素膠結固化后的土質(zhì)邊坡，能有效防止坡面沖溝發(fā)育，減緩邊坡水土流失，增強其抗沖刷性能力。

近年來，基于數(shù)值模擬探索土體沖蝕形成機理成為水土流失研究新的熱點。由于顆粒離散元方法在描述固體介質(zhì)材料破壞、運動等方面的優(yōu)勢，被廣泛應用于巖土工程、地質(zhì)災害成因研究中[9-14]。土體沖蝕過程是土體顆粒在水流作用下破壞并被挾帶運移的過程，是非常典型的流固耦合作用過程，基于顆粒離散元的流固耦合方法將能夠揭示這一過程中的機理。吳謙等[15]利用三維顆粒流軟件結合室內(nèi)降雨模型試驗分析了傾角為70°黃土邊坡的沖刷過程，基于顆粒運動路徑、土體孔隙率變化等參數(shù)揭示了降雨對邊坡侵蝕能力的分布規(guī)律。張雁等[16]，柯云斌等[17]基于二維顆粒流方法分別探索了黃土和殘坡積土邊坡的沖刷破壞機理。離散元數(shù)值模擬在水土流失研究方面有了初步成果，但是從細觀尺度揭示土顆粒與流水相互作用的過程的研究尚少，尤其是尚未開展流水直接沖刷坡面的三維模擬研究。中緬油氣管道龍陵段管道作業(yè)邊坡回填土沖蝕過程其實就是流水夾帶土壤顆粒形成的顆粒流流失過程，本文基于顆粒離散元構建的流固耦合算法，采用數(shù)值模擬方法，從細觀尺度研究水動力條件與土體力學性質(zhì)對土體沖蝕破壞過程的影響，分析流水侵蝕的起動條件，為該類型土體坡面水土流失治理提供理論依據(jù)。

1邊坡回填土物理力學性質(zhì)及沖蝕特征

研究區(qū)位于滇西龍陵縣，地處怒江、龍川江兩江之間，滇西縱谷南段，高黎貢山南延部分，地勢呈中部高而東西斜勢，海拔高程一般在1 600 m左右，屬低-中高山區(qū)。年平均降水量2 105.7 mm，月最大降水量726.6 mm，日最大降水量134.7 mm，最大單點降雨量56.9 mm。雨季主要集中于5～10月，降水量占全年降水量的82.6%。

研究區(qū)范圍內(nèi)地層巖性比較單一，油氣管道沿線出露主要地層由新至老主要有沖洪積層（Q4al+pl）、殘坡積黏性土（Q4el+dl）;中生界（T、J、K）地層，為一套紫紅色、雜色的泥巖、長石石英砂巖、粉砂巖、夾煤層。燕山期花崗巖、混合巖、粗面巖在研究區(qū)附近大面積出露。管道穿越段主要為花崗巖分布區(qū)。坡面表層為網(wǎng)紋狀全風化花崗巖層（粗砂土層），全風化花崗巖中礦物成分最高為石英，其次是長石、高嶺土等，原巖中的長石、云母已經(jīng)分解，石英因抗風化能力強而保留下來。全風化花崗巖物理性質(zhì)與砂土類似，但顆粒粒度更大，孔隙比大，裂隙發(fā)育，具有一定崩解性。管道作業(yè)邊坡回填土厚度一般為3 m，局部地段有變化?；靥钔林饕栽氐娜L化花崗巖回填，以砂礫土為主，土體中夾雜著大量的花崗巖礫石和石英礫石。礫石粒徑介于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%。通過對研究區(qū)邊坡回填土取樣并進行室內(nèi)試驗，得到土體的物理力學參數(shù)，如表1所列。

輸油氣管道的鋪設過程中，需開挖15～25 m寬的作業(yè)邊坡，管溝鋪設管道后進行回填，回填土為早期開挖的土或就地取材。管溝一般深2～3 m，因作業(yè)邊坡的開挖和回填，不僅直接剝離了抗沖蝕性較強的本來就比較薄的土壤層，使砂土層直接出露于地表，而且破壞了原有的周邊環(huán)境和生態(tài)。這種天然的原生巖土體條件和生態(tài)環(huán)境一旦破壞，不容易恢復。同時，回填土由于各種粒徑顆?；祀s，固結性差，顆粒間黏結性差，容易被沖刷，使沖蝕從天然狀態(tài)下處于相對平衡狀態(tài)的較為緩慢的自然侵蝕過程，在極短的時間迅速轉(zhuǎn)變?yōu)閺娏业墓こ绦詻_蝕過程。根據(jù)現(xiàn)場調(diào)查，龍陵油氣管道穿越區(qū)坡面沖蝕類型主要包括雨滴濺蝕、面蝕、細溝沖蝕、淺溝沖蝕、沖溝沖蝕等（見圖1）。

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

2.1顆粒流方法

Cundall提出了可以分析顆粒材料力學性態(tài)的離散單元法，從而可以對顆粒團粒體的穩(wěn)定、變形及本構關系進行分析，對固體力學大變形問題具有良好的適應性[18]。它是將顆粒介質(zhì)的運動及其相互作用通過圓形（或異型）離散單元來模擬。每一時刻顆粒在平面內(nèi)的位置和速度由平動和轉(zhuǎn)動運動方程來確定。

PFC3D（Particle Flow Code 3 Dimension）即三維顆粒流程序，是以球形顆粒介質(zhì)的運動及其相互作用為研究對象的一種離散單元方法。該方法通過將實際的材料離散成由若干個剛性顆粒組成的模型，來研究實際材料的各種力學特性。該程序最初被用來研究顆粒介質(zhì)特性，它利用有代表性的大量顆粒單元來表示實際物體，并利用這種局部的模擬結果來研究連續(xù)介質(zhì)邊值問題。PFC 3D顆粒流程序基于細觀力學通過離散單元方法來模擬球形顆粒介質(zhì)的運動及其相互作用，以牛頓第二定律和力–位移定律為基礎，采用顯式差分算法，循環(huán)應用牛頓定律分析顆粒的運動特征、力–位移定律分析顆粒的接觸特征，從而實現(xiàn)對顆粒介質(zhì)的運動及其相互作用過程的模擬[11-12，14]，控制方程如公式（1）所示。

3沖蝕過程數(shù)值模擬

3.1數(shù)值模型與參數(shù)

邊坡坡面的沖蝕是一個復雜的過程，涉及到降雨、土體飽和、雨水匯集等連續(xù)的環(huán)節(jié)。在降雨過程中，坡面土體通常會經(jīng)歷雨滴濺蝕、片蝕、細溝侵蝕、切溝侵蝕等階段。這些階段的產(chǎn)生與降雨強度、坡面形態(tài)及坡面土體本身的物理力學性質(zhì)相關。降雨強度及坡面形態(tài)主要影響坡面侵蝕中的水動力條件。本次模擬研究中僅考慮水力條件和土體力學性質(zhì)對坡面土體沖刷過程的影響，不討論水動力條件成因。構建一個100 mm×100 mm×100 mm的流固模型（見圖2），包括流體計算網(wǎng)格和固體顆粒。模型中流體計算網(wǎng)格規(guī)模為10×10×10，每個網(wǎng)格單元大小為10 mm×10 mm×10 mm;土體整體厚度為50 mm，右側(cè)有個高為20 mm土坎。模型中土體由球顆粒構成，其粒徑介于2.0～3.2 mm，顆粒數(shù)量為26 850個，顆粒之間采用黏結模型進行連接。模型中的參數(shù)如表2所列。

3.2水流沖蝕模擬結果

模擬過程中，固體部分首先在重力作用下達到力學平衡狀態(tài)，形成初始的土體物質(zhì)基礎;基于模型試驗及前人經(jīng)驗，設定模型左側(cè)入水口初始水流速度分別為2.0，1.0 m/s和0.5 m/s，模型右側(cè)排水口流體壓力為0 Pa，形成初始的穩(wěn)定流場;之后，開啟流固耦合計算。圖3～5給出了3種水流速度作用下模型中土體顆粒運移距離、流體速度分布與流體壓力分布隨時間的演化過程。由圖3～5可見，凸起的土坎在3種流速的水流侵蝕作用下逐漸被沖毀。當入口水流速度較大時（流速為2.0 m/s），在較大的流體拖拽力作用下土坎被迅速沖毀，呈現(xiàn)散開狀被挾帶走，并且溝底的部分土體顆粒也被流水刮起;隨入口水流速度降低，土坎被沖毀的速度降低，沖毀的土體局部還黏結在一起，流水對溝底的侵蝕現(xiàn)象不明顯。由此可見，水動力條件是影響土體沖蝕過程的關鍵參數(shù)之一，降低水動力條件是防止邊坡土體水毀的重要手段之一。

除水流速度之外，土體沖蝕破壞與其自身的強度特征關系密切。在模擬時，通過提高模型中顆粒間的黏結強度來改善土體的力學性能，進而分析土體力學特征對沖蝕過程的影響。圖6展示了土體黏結強度提高10倍之后的模擬結果。與圖3～5相比發(fā)現(xiàn)，當入口流

速為2.0 m/s時，土體黏結強度增大后土坎整體被流水沖走，但是未出現(xiàn)溝底侵蝕的現(xiàn)象;當入口流速降為0.5 m/s時，土坎未產(chǎn)生沖蝕破壞現(xiàn)象。由此可見，改善土體本身的強度特征有助于防止土體水毀現(xiàn)象的產(chǎn)生。

上述模擬從細觀尺度探討水動力條件與土體力學性質(zhì)對土坎沖蝕破壞過程的影響。當模擬水流速度大時，土坎以散開形式被迅速沖毀，溝底的部分土體顆粒也被流水刮擦帶走;降低水流速度后，土坎以團塊的形式被沖走，并沒有被沖毀。另外，提高模型中顆粒間的黏結強度能夠改善土體的力學性能，進而提高土體的抗沖蝕的能力。

4結 語

本文基于Navier-Stokes方程和牛頓第二定律構建的顆粒離散元流固耦合數(shù)值模擬方法，以強風化花崗巖區(qū)的土體為研究對象構建了沖溝數(shù)值模型，研究了邊坡回填土的沖蝕過程。研究表明本文中構建的數(shù)值模擬方法能夠較好地再現(xiàn)土體沖蝕過程。研究還發(fā)現(xiàn)降低坡面流水的動力條件和改善土體本身的強度特征將有助于防止土體流水侵蝕現(xiàn)象的發(fā)生，對于邊坡水土流失治理具有一定的啟示意義。本文關于土體沖蝕的研究獲得了一些定性規(guī)律的認識，下一步將結合具體降雨特征、匯水區(qū)面積、邊坡的土體特征、坡長和坡高等幾何特征，獲得邊坡土體沖蝕的啟動流速等定量數(shù)據(jù)，以指導邊坡水土流失的防治工作。

參考文獻：

[1]李宗善，楊磊，王國梁，等.黃土高原水土流失治理現(xiàn)狀、問題及對策[J].生態(tài)學報，2019，39（20）：7398-7409.

[2]白皓.紅土工程堆積體坡面水蝕過程及泥沙顆粒分布特性研究[D].楊凌：西北農(nóng)林科技大學，2018.

[3]王玉寬，文安邦，張信寶.長江上游重點水土流失區(qū)坡耕地土壤侵蝕的～（137）Cs法研究[J].水土保持學報，2003，17（2）：77-80.

[4]王玲.陡坡地水蝕過程與泥沙搬運機制[D].北京：中國科學院大學，2016.

[5]湯珊珊.覆沙坡面水蝕產(chǎn)沙動力過程與微地貌變化的響應研究[D].西安：西安理工大學，2018.

[6]馮秀，查軒，黃少燕.人工模擬降雨條件下花崗巖紅壤坡面侵蝕過程與特征分析[J].中國水土保持科學，2014，12（1）：19-23.

[7]杜婷婷，李志清，王曉明，等.黃土邊坡降雨沖刷模型試驗研究[J].工程地質(zhì)學報，2018，26（3）：732-740.

[8]裴向軍，楊晴雯，許強，等.改性鈉羧甲基纖維素膠結固化土質(zhì)邊坡機制與抗沖蝕特性研究[J].巖石力學與工程學報，2016，35（11）：2316-2327.

[9]王學良，孫娟娟，周書，等.尾礦庫潰壩運動特征模擬研究[J].工程地質(zhì)學報，2019，27（1）：144-151.

[10]眭靜，姜元俊，樊曉一，等.碎屑流沖擊剛性擋墻的力學模型研究[J].巖石力學與工程學報，2019，38（1）：121-132.

[11]費建波，介玉新，張丙印，等.土的三維破壞準則在顆粒流模型中的應用[J].巖土力學，2016，37（6）：1809-1817.

[12]朱涵成，韓文喜，陳超.砂巖常規(guī)三軸的顆粒流數(shù)值模擬[J].地質(zhì)災害與環(huán)境保護，2013，24（3）：104-108.

[13]徐金明，謝芝蕾，賈海濤.石灰?guī)r細觀力學特性的顆粒流模擬[J].巖土力學，2010，31（增2）：390-395.

[14]楊慶華，姚令侃，楊明.地震作用下松散堆積體崩塌的顆粒流數(shù)值模擬[J].西南交通大學學報，2009，44（4）：580-584.

[15]吳謙，王常明，宋朋燃，等.黃土陡坡降雨沖刷試驗及其三維顆粒流流-固耦合模擬[J].巖土力學，2014，35（4）：977-985.

[16]張雁，高樹增，閆超群，等.降雨對黃土路基邊坡的沖刷規(guī)律[J].中國地質(zhì)災害與防治學報，2017，28（4）：34-39.

[17]柯云斌，符元帥，闕云，等.降雨誘發(fā)殘積土陡坡坡面沖刷破壞顆粒流模擬[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），F(xiàn)rance，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].系統(tǒng)仿真學報，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