葛華　方迎潮　趙飛　孔志崗　戴國(guó)文　張垚

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

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

中圖法分類(lèi)號(hào)： P66

文獻(xiàn)標(biāo)志碼： A

DOI：10.16232/j.cnki.1001-4179.2021.11.020

0引 言

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

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

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

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

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

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

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

2.1顆粒流方法

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

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

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

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

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

3.2水流沖蝕模擬結(jié)果

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

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

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

上述模擬從細(xì)觀尺度探討水動(dòng)力條件與土體力學(xué)性質(zhì)對(duì)土坎沖蝕破壞過(guò)程的影響。當(dāng)模擬水流速度大時(shí)，土坎以散開(kāi)形式被迅速?zèng)_毀，溝底的部分土體顆粒也被流水刮擦帶走;降低水流速度后，土坎以團(tuán)塊的形式被沖走，并沒(méi)有被沖毀。另外，提高模型中顆粒間的黏結(jié)強(qiáng)度能夠改善土體的力學(xué)性能，進(jìn)而提高土體的抗沖蝕的能力。

4結(jié) 語(yǔ)

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

參考文獻(xiàn)：

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

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

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

[4]王玲.陡坡地水蝕過(guò)程與泥沙搬運(yùn)機(jī)制[D].北京：中國(guó)科學(xué)院大學(xué)，2016.

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

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

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

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

[9]王學(xué)良，孫娟娟，周書(shū)，等.尾礦庫(kù)潰壩運(yùn)動(dòng)特征模擬研究[J].工程地質(zhì)學(xué)報(bào)，2019，27（1）：144-151.

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

[11]費(fèi)建波，介玉新，張丙印，等.土的三維破壞準(zhǔn)則在顆粒流模型中的應(yīng)用[J].巖土力學(xué)，2016，37（6）：1809-1817.

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

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

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

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

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

[17]柯云斌，符元帥，闕云，等.降雨誘發(fā)殘積土陡坡坡面沖刷破壞顆粒流模擬[J].濟(jì)南大學(xué)學(xué)報(bào)（自然科學(xué)版），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]李識(shí)博，王常明，王鋼城，等.松散堆積物壩基滲透淤堵試驗(yàn)及顆粒流模擬[J].水利學(xué)報(bào)，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]宿崇，侯俊銘，朱立達(dá)，等.基于流固耦合算法的單顆磨粒切削仿真研究[J].系統(tǒng)仿真學(xué)報(bào)，2008，20（19）：5250-5253.

[22]方志，陳育民，何森凱.基于單相流的減飽和砂土流固耦合改進(jìn)算法[J].巖土力學(xué)，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