人工黃土坡面重力侵蝕數(shù)值模擬

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

(1.桂林理工大學(xué)廣西環(huán)境污染控制理論與技術(shù)重點實驗室，廣西 桂林 541004；2.桂林理工大學(xué)巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心，廣西 桂林 541004

0 引 言

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

1 實驗材料與方法

1.1 實驗概況

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

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

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

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

1.2.1 FLAC 3D數(shù)值模擬基本原理

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

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

(1)

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

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

(2)

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

圖1 四面體Fig.1 Tetrahedron

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

(3)

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

(4)

對于靜態(tài)問題，在不平衡力中加入非黏性阻尼，使得系統(tǒng)振動逐漸衰減至達到平衡狀態(tài)。此時式3變?yōu)椋?/p>

(5)

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

(6)

式中：∝為阻尼系數(shù)，默認(rèn)值為0.8。

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

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

1.2.2 Rhino三維建模軟件介紹

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

1.2.3 復(fù)雜地形的FLAC 3D建模方法

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

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

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

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

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

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

1.3.2 黃土物理力學(xué)指標(biāo)

根據(jù)邊坡工程經(jīng)驗、現(xiàn)場資料分析、現(xiàn)場及室內(nèi)巖土物理力學(xué)試驗和有限元計算的需要[25]，模型土層材料物理力學(xué)參數(shù)的具體特征取值見表1。

表1 黃土坡面性質(zhì)參數(shù)Tab.1 Properties of loess slope

2 結(jié)果與分析

2.1 黃土坡面重力侵蝕發(fā)育情況

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

圖4 監(jiān)測點位移情況Fig.4 Displacement of the monitoring point

2.2 黃土坡面位移分布規(guī)律

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

圖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 黃土坡面應(yīng)力場分布規(guī)律

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

圖9 坡面第一應(yīng)力分布圖Fig.9 Distribution of the first stress on the slope

圖10 坡面第三應(yīng)力分布圖Fig.10 Distribution of the third stress on the slope

2.4 黃土坡面塑性區(qū)分布規(guī)律

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

圖11 坡面塑性狀態(tài)分布Fig.11 Distribution of slope plastic state

從圖11中可知，以20°坡面為例，根據(jù)第一場次模型平衡狀態(tài)時的塑性區(qū)分布，對應(yīng)第二場次連續(xù)性降雨結(jié)束后侵蝕地形，可以看出坡面中下部出現(xiàn)明顯侵蝕細(xì)溝且隨著降雨進行，細(xì)溝開始溯源侵蝕，侵蝕地形符合第一場次模型平衡狀態(tài)時的塑性區(qū)分布規(guī)律。同時對三次連續(xù)性降雨坡面進行受力分析，發(fā)現(xiàn)本場次侵蝕地形基本符合上一場模型平衡狀態(tài)時的塑性區(qū)分布情況，可以說明FLAC3D在一定程度上能夠?qū)η治g地形發(fā)育起到預(yù)測的作用。

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

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

3 結(jié) 論

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

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

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

(3)坡面位移以向下侵蝕滑落為主，整體位移最大部分出現(xiàn)在坡頂及坡面上部，縱向最大位移出現(xiàn)在坡面下部；坡面土體應(yīng)力分布主要以受壓屈服為主，且隨坡度增大導(dǎo)致應(yīng)力集中，土體內(nèi)部主應(yīng)力及最大受壓屈服區(qū)體積增加，加劇重力侵蝕發(fā)生。

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

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

□