人工黃土坡面重力侵蝕數值模擬

金 鑫，劉 喆，王春振，宋 穎，趙華榮，王 迪，李 霞

(1.桂林理工大學廣西環境污染控制理論與技術重點實驗室，廣西 桂林 541004；2.桂林理工大學巖溶地區水污染控制與用水安全保障協同創新中心，廣西 桂林 541004

0 引 言

當地表土石在地震、降雨、徑流和海浪、風、凍融、人工采掘等任何一種應力作用時，在重力作用下，而產生遷移和堆積的一種侵蝕過程稱為重力侵蝕[1]。重力侵蝕是我國黃土高原地區土壤侵蝕的三種主要類型之一[2-5]，但由于影響重力侵蝕的因素較多，且其發生具有隨機性和突發性，因而較難準確的測定產沙量[6]。相較于對水力侵蝕過程機理的認識，對重力侵蝕的研究還較為薄弱，侵蝕機理尚未完全明確[7]。近年來，隨著研究技術和方法的不斷應用，對重力侵蝕的認識逐漸加深。例如，譚萬沛等[8]探究了不同的坡度條件對重力侵蝕的影響，指出重力侵蝕的發生形式和侵蝕產沙量在不同坡度條件下存在較大差異。Gabet[9]發現重力侵蝕的出現次數在不同下墊面植被類型條件下有較大波動，且重力侵蝕的發生次數在下墊面植物根系具有較深的深度及較高的強度時，而出現了明顯減少。鄭書彥[10,11]等對重力侵蝕的發生機理分別從滑坡等土壤侵蝕角度進行了探討。朱同新[12]等研究指出在不同的空間位置處和不同的時間階段內，重力侵蝕發生強度和類型是不同的。此外，王光謙等[13]研究了溝坡沖蝕條件下重力侵蝕的力學機制，并通過力學方法建立了坡溝重力侵蝕的概化力學模型，李容全[14]通過14C法計算了在重力侵蝕下溝谷后退速度，葉浩[15]通過GPS研究了砒砂巖地區重力侵蝕引起的溝緣線后退速度等。朱同新等[16-18]等對黃土高原重力侵蝕發生現象進行了觀察和統計分析，楊吉山[19]通過對橋溝小流域的研究指出不同的溝道發育階段，重力侵蝕的特點也不盡相同。在大多數學者的研究中，重力侵蝕主要是在不同的地形、地貌、植被覆蓋等條件下進行研究。但從溝坡系統力學穩定性的角度來看，由于溝坡系統重力侵蝕具有隨機性和多發性，基礎觀測資料相對較少，相關研究相對薄弱[20]。因此，本文通過Rhino和FLAC 3D數字模擬軟件的二次開發，建立人工黃土坡面侵蝕概化模型，并對黃土坡面重力侵蝕破壞區進行研究，為今后的水土流失相關研究提供參考。

1 實驗材料與方法

1.1 實驗概況

本實驗用土為陜西榆林岔巴溝黃土，屬于第四紀時期形成的土壤堆積物，含多量鈣質或黃土結核，在干燥時較堅硬，被水流潤濕后，更容易剝落和遭受侵蝕，持水性一般，土體垂直節理發育，平均容重為1.15 g/cm3。本實驗用土地區選自岔巴溝黃土侵蝕典型代表地區作為實驗用土采集地，采用機械挖掘表層土樣后，將土樣裝車運回。土樣采集回來放置于庫房將土壤經風干后，過1 cm的篩網，剔除天然雜草、石礫以及較大的土壤結塊。為保證土槽的透水性及透水狀況接近于天然坡面，先在土槽底部鋪上10 cm細沙，并蓋上透水紗布。根據土槽的體積及所需容重進行分層填土，以10 cm為一層并壓實，每填完一層把表面刮毛以減少各層間的差異性。完成一組實驗后，將表層10 cm土壤刨松，然后填入新土進行下一組實驗。

每次正式降雨實驗前一天，選用30 mm/h降雨強度進行前期預降雨，降雨至坡面有積水出現，并用塑料布覆蓋并靜置24 h。前期降雨的目的一是保持坡面下墊面有相同的前期土壤含水量，二是減少下墊面條件的空間變異性。

實驗中降雨系統采用全自動不銹鋼模擬降雨器，該降雨裝置配有旋轉下噴式噴頭，該設備采用潛水泵提供動力，其揚程為50 m，降雨高度為6 m，雨強可調節范圍為10～200 mm/h，雨滴直徑為0.1～6 mm。實驗采取兩個坡度(20°、25°)，一個雨強(60 mm/h)。每場次降雨持續時間為90 min，每組實驗分為4個場次進行連續性降雨，因此使坡面細溝得到充分的發育，便于觀測在重力侵蝕作用下坡面細溝不同發育階段的形態特征。

1.2 人工黃土坡面重力侵蝕數字分析原理及建模方法

1.2.1 FLAC 3D數值模擬基本原理

FLAC 3D(Fast Lagrangian Analysis of Continua)是由Itasca公司研發推出的連續介質力學分析軟件，內置有靜力、動力、蠕變、滲流、溫度5種計算模式，各種模式間可以相互耦合，以模擬計算復雜的工程力學行為。作為有限差分軟件，相對于其他有限元軟件，在算法上有以下優點：FLAC 3D采用混合離散法[21]來模擬材料的塑性破壞和塑性流動，較通常采用的離散集成法更為準確合理；FLAC 3D采用顯示差分法求解微分方程，可以方便地求得應力增加、不平衡力并跟蹤系統的演化過程[22]。

在FLAC 3D中采用混合離散方法，區域將被離散為常應變多面體單元[23]，在運算過程中，每個多面體又進一步離散為以該多面體頂點為頂點的常應變四面體，所有變量在四面體上計算，取多面體單元內四面體應力、應變的加權平均值作為多面體單元的應力、應變值。如圖1所示，對任意一個四面體，節點編號為1～4，第n面表示與節點n相對的面，設其內任一點的速率分量為vi，則可由高斯公式得

(1)

式中：V為四面體的體積；S為四面體的外表面；nj為外表面的單位法向向量分量。

對于常應變單元，vi為線性分布；nj在每個面上為常量，由式(2)可得：

(2)

式中：l為節點l的變量；(l)為面l的變量。

圖1 四面體Fig.1 Tetrahedron

在FLAC 3D中計算的對象以節點為主，將力與質量集中在節點上，通過運動方程進行求解，公式如下。

(3)

FLAC3D由速率求解某一步時長的單元應變增量，公式如下。

(4)

對于靜態問題，在不平衡力中加入非黏性阻尼，使得系統振動逐漸衰減至達到平衡狀態。此時式3變為：

(5)

其二，常見的書面體育文本主要包括體育新聞報道、體育教學資料、體育營銷文案、體育學術論文等，所涉體育項目眾多（羅永洲，2012）。盡管如此，體育類文本也呈現出一定共有特征。體育文本往往牽扯到國際賽事、國際組織、運動員姓名、比賽規則、技術要領等方面，具有專業性強，專有名詞、專業術語、縮略語、合成詞繁多等特征。因此，體育文本的翻譯過程需要處理大量的特殊詞匯。鑒于體育文本的這一特點，相對于其他類型的翻譯而言，體育類翻譯對譯員的專業知識和詞匯的積累要求更高。

(6)

式中：∝為阻尼系數，默認值為0.8。

由上可知FLAC 3D的計算過程，如圖2所示。

圖2 FLAC 3D計算流程Fig.2 Calculation process of FLAC 3D

1.2.2 Rhino三維建模軟件介紹

Rhino是美國Robert Mcneel公司開發的一款基于NURBS為主體結構的專業三維建模軟件，采用自由曲面的建模技術和特征實體的操作模式[24]。相較于POLYGON(多邊形)，NURBS可用極少的控制點創建具有高精度的復雜曲線并建立曲面模型，可精確的實現自由造型產品模型制作，同時支持ACIS、Parasolid、3ds、dwg以及點云資料等多種格式文件的輸入輸出。在Rhino建模之前為確保模型精度，在模型核心參數公差值的選擇上，采取Absolute tolerance默認的絕對公差值，默認Absolute tolerance值為0.01 mm。

1.2.3 復雜地形的FLAC 3D建模方法

雖然FLAC 3D專為巖土工程力學開發，內置有豐富的彈、塑性材料本構模型。但是在建立復雜計算模型時，仍是以命令驅動模式為主要輸入模式，這種方式對于模型的建立過于繁雜，因此本研究采用3D造型軟件Rhino結合FLAC 3D內置的fish語言在初始單元模型基礎上編寫了前處理程序，實現了對復雜多層地形建模的二次開發[25]。

首先使用Rhino讀入原始侵蝕坡面點云數據，采用MeshPatch指令進行網格補丁，建立三角網格，對網格進行修剪、調整網格坐標原點位置、調整模型單位為米，生成STL文件類型的地表模型[圖3(a)]。在FLAC 3D中通過fish語言程序提取STL文件的地表高程信息，然后使用generate zone brick生成六面體實體網格，采用全局坐標系定義坐標軸x、y、z方向網格單元數目及相鄰單元尺寸大小比率，然后調用generate topography geomset命令填充地表模型與新建六面體實體網格之間的間隙，建立接觸面。采用Rhino軟件較好地解決了在FLAC 3D中前期復雜侵蝕地形建模困難的問題，成功實現建模。仿真實例表明，所建立的三維模型能夠真實反映侵蝕坡面地形，模擬效果良好。

1.3 人工黃土坡面重力侵蝕的有限元數值模擬

1.3.1 人工黃土坡面概化模型及有限元計算模型

圖3(a)為采用Rhino軟件繪制的侵蝕坡面地形圖，圖3(b)為采用FLAC 3D建立的侵蝕坡面概化模型。該模型底部和四周設為固定約束邊界，坡面設為自由邊界[25]。該模型土層設置為單一土層，平均厚度為1 m，長度設為4 m，寬度設為1.2 m。坡面坡度為20°、25°，20°模型有限元網格共有節點445 511 個，單元40 萬個；25°模型有限元網格共有節點445 511 個，單元40 萬個。

圖3 20°及25°人工黃土坡面地形圖及概化模型Fig.3 Topographic maps and generalized models of 20° and 25° artificial loess slopes

1.3.2 黃土物理力學指標

根據邊坡工程經驗、現場資料分析、現場及室內巖土物理力學試驗和有限元計算的需要[25]，模型土層材料物理力學參數的具體特征取值見表1。

表1 黃土坡面性質參數Tab.1 Properties of loess slope

2 結果與分析

2.1 黃土坡面重力侵蝕發育情況

分別在坡面0.9、1.6、2、3 m處設置監測點，觀測坡面上中下3個部位處位移變化情況。模型運行到24 000 步時，計算最大不平衡力低于系統設置的默認值(10-5)，模型計算停止。從圖4中可以看出，各點的位移變化趨勢相同，隨著步長的增加，各點位移增大，然后降低至某一值后保持穩定。觀察圖4位移曲線變化規律可劃分為不同的侵蝕發育時期，將開始至位移最高點階段劃分為侵蝕發育期，從最高點降至穩定值階段劃分為侵蝕發育成熟期，從穩定值以后劃分為穩定期，由圖4可見，坡面下部侵蝕發育完成時間要快于中上部侵蝕發育，這與實驗過程觀察結果一致。且從圖4中可知，在坡面上部0.9 m處整體位移最大，坡面下部3 m處整體位移最小，當坡度增大時，中下部位監測點位移均有一定程度增大，其中1.6、2、3 m處增加幅度分別為10.7%、55.4%、33.4%。上部監測點位移減小，減小幅度為4.5%。分析原因為，坡度增加導致坡面徑流速度增大、坡面承雨面積減小，從而在一定程度上加劇了坡面中下部侵蝕發育程度，減緩了坡面上部的侵蝕發育。

圖4 監測點位移情況Fig.4 Displacement of the monitoring point

2.2 黃土坡面位移分布規律

在最大不平衡力低于系統默認值時，系統達到平衡狀態，得出模型3個方向的位移大小和坡面整體位移、豎直方向位移和水平方向位移分布圖(圖5、圖6、圖7、圖8)。從整體位移云圖來看，位移最大部分出現在坡頂及坡面上部；豎直方向位移云圖與總體位移云圖分布規律基本一致，且豎直方向最大位移數值也出現在坡頂及坡面上部，由此可知坡面位移是以向下侵蝕滑落為主。由圖7見，最大水平位移出現在坡面下部溝坡地帶，而坡面其他位置的水平位移基本為零，這表明相對于其他位置，下部出現水平位移的溝坡處易在水平方向發生侵蝕坍塌或滑坡，而朝溝底方向滑動。由圖8可知，縱向最大位移出現在坡面中下部，這與侵蝕實驗觀測結果一致，侵蝕實驗中坡面中下部侵蝕發育最為劇烈。

圖5 20°及25°坡面整體位移圖Fig.5 Overall slope displacement of 20° and 25°

圖6 20°及25°坡面豎直方向位移圖Fig.6 Vertical displacement of 20° and 25° slope

圖7 20°及25°坡面水平方向位移圖Fig.7 Horizontal displacements of 20°and 25°slopes

圖8 20°及25°坡面縱向位移圖Fig.8 Longitudinal displacements of 20°and 25°slopes

2.3 黃土坡面應力場分布規律

FLAC 3D分析計算出坡面模型在平衡狀態時所受到的應力大小及分布規律。如圖9、圖10中所示(FLAC 3D中定義為以拉應力為正，壓應力為負)，應力分布基本一致，從坡面應力圖分布來看，未出現拉應力，因此坡面受力基本以壓應力為主，即為坡面受到破壞時主要以壓剪形式產生的破壞模式為主。從圖9可以看出，邊坡內部，最大主應力(即第一主應力)方向與水平面夾角逐漸減小，且隨著坡度的增加土體內部受到的壓應力也隨之增大，表明深層土體主要受鉛直方向的壓應力，表現為受壓屈服。從圖10分析可知，由于坡度增大，導致應力集中，土體內部主應力及最大受壓屈服區體積增加，因此，在一定程度上坡度的增加，加大了土體的受力，從而導致加劇了重力侵蝕的發生。

圖9 坡面第一應力分布圖Fig.9 Distribution of the first stress on the slope

圖10 坡面第三應力分布圖Fig.10 Distribution of the third stress on the slope

2.4 黃土坡面塑性區分布規律

圖11為模型計算達到平衡狀態時坡面塑性狀態指示圖，圖中主要獲取剪切屈服區域(shear)和張拉屈服區(tension)[26]。每個屈服區函數賦有兩種狀態：now和past，其中now表示該單元在本次計算步數中正處于屈服面上，past表示曾處于屈服面上，但現處于彈性范圍，因此只有處于now狀態單元才對模型的破壞起作用[27]。由圖11可見，張拉屈服區主要分布在坡面頂部及坡面中上部位，說明在坡頂及坡面中上部位發生張拉破壞的可能性較大，張拉破壞較為嚴重；剪切屈服區域分布較為廣泛，其主要分布在坡面上下部溝坡區域及土體內部，說明這些區域易發生水平剪切變形，且當坡度增加時，坡面上部位剪切屈服區減少，張拉屈服區有所增加。

圖11 坡面塑性狀態分布Fig.11 Distribution of slope plastic state

從圖11中可知，以20°坡面為例，根據第一場次模型平衡狀態時的塑性區分布，對應第二場次連續性降雨結束后侵蝕地形，可以看出坡面中下部出現明顯侵蝕細溝且隨著降雨進行，細溝開始溯源侵蝕，侵蝕地形符合第一場次模型平衡狀態時的塑性區分布規律。同時對三次連續性降雨坡面進行受力分析，發現本場次侵蝕地形基本符合上一場模型平衡狀態時的塑性區分布情況，可以說明FLAC3D在一定程度上能夠對侵蝕地形發育起到預測的作用。

通過編程FISH語言，調用命令流對20°、25°兩個模型塑性屈服區體積進行計算，結果如表2所示。可以看出，在黃土坡面上塑性屈服區主要以剪切屈服區為主，其中在20°坡面上，剪切屈服區占總屈服區體積的87.9%，張拉屈服區占總屈服區體積的12.1%；在25°坡面上，剪切屈服區占總屈服區體積的85.4%，張拉屈服區占總屈服區體積的14.6%，說明坡面主要以剪切破壞模式為主。當坡度增加至25°時，剪切屈服區體積增加了0.271 m3，張拉屈服區體積增加了0.166 m3，總屈服區體積增加0.437 m3，相較于20°坡度下剪切屈服區、張拉屈服區和總屈服區體積分別增加了7.6%、34.08%、10.8%。說明隨坡度增加，張拉屈服區破壞體積占總屈服區破壞體積增大，但剪切塑性屈服區的體積占全部屈服區體積的大部分，因此重力侵蝕破壞仍是以剪切破壞模式為主。

表2 20°及25°塑性屈服區體積Tab.2 Plastic yield zone volumes of 20° and 25°

3 結 論

本文通過對FLAC 3D數字模擬軟件的二次開發建立人工黃土坡面侵蝕概化模型，并對黃土坡面重力侵蝕破壞區進行研究，得出以下結論：

(1)不同于其他學者在研究重力侵蝕時多采用的統計分析和GIS等模型方法，本研究通過可視化Rhino軟件和有限差分軟件FLAC3D結合，實現在復雜侵蝕坡面上的三維建模，解決了在FLAC 3D中復雜侵蝕地形建模困難的問題，能夠良好模擬侵蝕坡面地形。

(2)通過監測點位移規律，將侵蝕分為侵蝕發育階段、侵蝕成熟階段、侵蝕穩定階段，坡面下部侵蝕發育階段要早于坡面中上部侵蝕發育完成。

(3)坡面位移以向下侵蝕滑落為主，整體位移最大部分出現在坡頂及坡面上部，縱向最大位移出現在坡面下部；坡面土體應力分布主要以受壓屈服為主，且隨坡度增大導致應力集中，土體內部主應力及最大受壓屈服區體積增加，加劇重力侵蝕發生。

(4)隨坡度增加，張拉屈服區破壞體積占總屈服區破壞體積增大，但剪切塑性屈服區的體積占全部屈服區體積的大部分，因此重力侵蝕破壞仍是以剪切破壞模式為主。

(5)對坡面進行受力分析，得出每場次結束侵蝕地形基本符合上一場次模型平衡狀態時的塑性區分布，因此FLAC3D在一定程度上能夠對侵蝕地形發育起到預測的作用。

□