廖優(yōu)斌， 蔣明鏡， 陳有亮， 奚邦祿, 譚亞飛鷗，3

(1.上海理工大學 環(huán)境與建筑學院，上海 200093； 2.同濟大學 土木工程學院，上海 200092；3.同濟大學 巖土及地下工程教育部重點實驗室， 上海 200092)

1 研究背景

20世紀50年代起，美國的Survey計劃、Apollo計劃及前蘇聯(lián)的Luna計劃開啟了人類探月計劃的第一輪高潮，揭開了月球神秘的面紗[1],發(fā)現(xiàn)月壤的形成機理、所處的重力場環(huán)境以及力學特性與地球土壤差別較大。21世紀，歐盟的曙光計劃、印度的飛船計劃、日本的月亮女神計劃和中國的嫦娥探月計劃掀起了第二輪探月高潮[2]，將人類對月球的認識推到了一個新的高度。通過更高分辨率的影像數(shù)據(jù)發(fā)現(xiàn)月球表面存在不少的高陡邊坡，其坡度一般都大于30°，月海區(qū)可達48.9°，月陸區(qū)甚至可達55.7°[3]。根據(jù)Apollo15探測器對月球表面雨海東部邊緣現(xiàn)場影像數(shù)據(jù)分析，部分地區(qū)發(fā)現(xiàn)了滑坡運動的證據(jù)[4]。人類活動常使自然邊坡的坡度變陡，一旦進入對月球資源的大規(guī)模開發(fā)，若工程設計不合理，就極易誘發(fā)坡體滑塌，甚至整體坍塌，這些地質(zhì)災害都會嚴重影響未來對月球資源開采的安全。我國的嫦娥探月計劃也明確指出當探測第3期工程基本完成不載人的月球探測任務后，隨后將實施載人登月和月球基地建設[2]。無論是資源開采還是基地建設，這些活動不可避免地將遇到不同坡角下的坡體工程。

圖1 月壤邊坡

坡角的大小會直接影響地面物質(zhì)的流動及能量轉換的規(guī)模和強度，月壤形成坡體的坡度比常規(guī)的略大，如圖1所示。另外現(xiàn)場的溝槽試驗、深陷足跡邊緣月壤不塌落及粉土結塊等現(xiàn)象[5]也表明了月壤具有較高的自穩(wěn)能力，這些現(xiàn)象引起了國內(nèi)外學者的廣泛研究[3-5]。坡角的大小和材料休止角直接相關，而休止角的大小是否會受到重力場的影響至今尚無定論。月面環(huán)境下重力場大約為1.62 m/s2,僅僅為地面環(huán)境下的1/6[1]，Horgan等[6]認為休止角隨著重力減小而增大，而Nakashima等[5]認為休止角不受重力大小的影響。另一方面，月面的特殊環(huán)境(低重力、高真空及晝夜溫差大)促成了顆粒之間的微距力，從而在一定程度上提高了月壤的自穩(wěn)能力[7-8]。削坡試驗作為研究坡體穩(wěn)定性中最直接最簡單的方法，研究不同切坡傾角的削坡試驗對理解月面環(huán)境下坡體滑塌機理具有一定的意義。然而，月壤的削坡試驗的研究仍處于起步階段，相關的研究成果較少，主要集中在試驗研究和數(shù)值模擬。在試驗研究方面，Zheng等[9]采用模擬月壤物研究了不同傾角下的月壤滑坡試驗，研究發(fā)現(xiàn)月壤滑坡過程中出現(xiàn)裂隙面和光滑面，然而這種模擬月壤只是級配及材料相似，對于低重力場及表面微距力無法考慮。傳統(tǒng)室內(nèi)試驗研究月壤具有成本高、還原度低、且對滑坡過程微觀信息難以全面觀察的缺點，而數(shù)值試驗作為科學研究的一種重要手段，能有效克服以上缺點，Bui等[10]采用離散單元法，將顆粒之間的自鎖能力引入到顆粒接觸模型中來研究月壤坡體的滑塌，從微觀機理上較好地解釋了月壤坡體具有較高的穩(wěn)定性。雖然國內(nèi)外目前對月壤滑坡研究取得了一定的成果，但能夠考慮顆粒間的抗轉動作用及范德華力的月壤滑坡破壞的研究報告尚少見。

離散單元法能夠很好地模擬顆粒材料大變形問題，深入土體的微觀領域去研究其本質(zhì)規(guī)律，故本文采用離散單元法進行模擬滑坡試驗。借助PFC2D軟件，通過植入二維月壤接觸模型，模擬不同切坡傾角下月壤坡體滑坡試驗，通過觀察其滑動過程的變形特性和分析其滑坡機理，建立初始坡角與工程災害指標的關系，為設計月面環(huán)境下放坡開挖及防范和治理月面環(huán)境下滑坡災害提供參考。

2 月壤的微觀模型介紹

本文所采用的微觀模型是本團隊最近提出來的，主要考慮月壤顆粒間抗轉動作用及范德華力，模型的有效性可以參考文獻[11-12],以下對模型進行簡要的介紹。

(1)顆粒間的抗轉動作用：月壤的形成與地壤的形成存在較大的差異，主要是隕石和微隕石的撞擊、宇宙射線和太陽風粒子的持續(xù)轟擊以及大幅度的晝夜溫差使巖石熱脹冷縮導致的破碎[1]，顆粒多屬多棱角多氣孔結構，表面凹凸不平，具有很強的嚙合作用且具有一定的抗轉動作用。

(2)顆粒間的范德華力：月壤的土性類似砂土，傳統(tǒng)土力學認為其不應該具有黏聚力，但由于月球表面大氣層極其稀薄，白天大氣中分子含量僅為104個/cm2，屬于超真空狀態(tài)[1]，月壤顆粒表面吸附氣體分子層厚度遠小于地面土壤顆粒，在地面環(huán)境下無需考慮的范德華力在月面環(huán)境下應該給予考慮。

根據(jù)文獻[11-12]，模型的力學響應如圖2所示，其中Fv為顆粒之間的范德華力。法向力學響應主要由彈性力和范德華力兩部分組成，其中線性函數(shù)的斜率表征顆粒之間的彈性接觸，而其截距表征范德華力。彈性力會隨粒間重疊量增加而增加，但范德華力不受重疊量改變的影響，恒為常數(shù)，當顆粒間沒有接觸時彈性力和范德華力均為零。切向力學響應采用增量法計算，最大值不超過滑動摩擦力μKnUn。轉動力學響應也采用增量法計算，最大值不超過Knunβr/6。其中Kn為法向接觸剛度，Ks為切向接觸剛度，Km為扭轉剛度，un為法向重疊量，us為剪切位移，θ為相對轉動角度，μ為摩擦系數(shù)。

圖2 月壤微觀接觸模型的力學響應[11-12]

3 月壤邊坡模型建立

在采用離散單元法模擬月壤邊坡滑坡過程試驗中，若直接采用與真實月壤顆粒相同的級配，則顆粒粒徑跨度過大，所需顆粒數(shù)目非常多，考慮到計算機能力及計算效率，對月壤級配做合理簡化[11-12]，如圖3所示。根據(jù)Mitchell等[13]測得真實月壤的黏聚力值范圍為0.1～1.0 kPa、內(nèi)摩擦角范圍為30°～50°，通過模擬雙軸試驗對月壤宏觀參數(shù)進行標定，微觀參數(shù)如表1所示，得峰值黏聚力為0.2 kPa，峰值內(nèi)摩擦角約為42.6°[11-12]，處于真實月壤的范圍之內(nèi)。

月壤滑坡數(shù)值模擬模型如圖4所示，地基附近及下部邊界采用3道無摩擦的剛性墻模擬，頂部采用自由邊界模擬，與常見的離心機模型試驗一致。

為了制備相對密實均勻的地基，采用Jiang等[14]提出的分層欠壓法生成地基，生成的模型地基寬度為5 m，初始高度為2.5 m，如圖4(a)所示。顆?？倲?shù)為1.5×105，平面孔隙比為0.20。為了模擬月面環(huán)境真實重力場，試樣需要在(1/6)g重力場下完成固結平衡，固結完成后地基中水平應力及空隙比分布如圖5所示。由圖5可知固結平衡后的地基初始應力數(shù)值解和理論解基本相同，孔隙比分布均勻。而后需要對地基進行切坡來模擬不同初始坡角下的邊坡，以邊坡斜面中點為坐標原點，不斷改變傾角的同時保證坡后顆?？傮w積不變，刪除A-B-C-D區(qū)域內(nèi)的顆粒。圖4(b)所示為切坡完成后45°的坡體，在坡肩處、坡面處及坡趾處各放置1個監(jiān)測點。最后，在(1/6)g重力場下完成月壤滑坡系統(tǒng)的第二次平衡，觀察并記錄滑坡信息。

表1 離散元中的月壤微觀參數(shù)

4 結果分析

4.1 滑坡形態(tài)分析

圖6給出本文數(shù)值模擬及相關的室內(nèi)試驗[9]中坡度為90°時模擬月壤滑坡過程形態(tài)示意圖。不同坡角下滑坡變形規(guī)律基本一致，此處僅以90°為特例做進一步介紹。無論數(shù)值模擬還是室內(nèi)試驗的滑坡過程都經(jīng)歷3個階段：起滑階段、滑動階段及穩(wěn)定階段。起滑階段，數(shù)值模擬的坡肩處出現(xiàn)明顯的倒圓角，數(shù)值模擬中坡趾處輪廓線比室內(nèi)試驗大，這主要是因為室內(nèi)試驗采用的月壤黏聚力為1 kPa，比數(shù)值模擬0.2 kPa要大得多，此外，室內(nèi)試驗中模擬月壤僅干燥了12 h，可能殘留的水分增大了黏聚力。滑動階段，數(shù)值模擬和室內(nèi)試驗都出現(xiàn)葉片狀滑體，90°坡體在滑動階段出現(xiàn)前凸起上細下粗的葉片狀，這主要是因為在葉片中心線上下兩側含有一定黏聚力的土體顆粒運動速度不一致產(chǎn)生變形差異所致。

穩(wěn)定階段，此時滑動區(qū)域的顆粒速度較小，主體變形基本完成，數(shù)值模擬和室內(nèi)試驗的滑動區(qū)域都從葉片狀變?yōu)樗笞訝?，但兩者仍有區(qū)別，數(shù)值模擬滑坡穩(wěn)定后坡面較光滑且坡度較小，而室內(nèi)試驗坡面不光滑坡度也較大，這種差異主要有兩個原因：一是因為數(shù)值模擬采用的是簡化級配，粒徑跨度較小，形成坡面較光滑，而室內(nèi)試驗無需考慮計算效率問題，可以直接采用真實的月壤級配，其粒徑跨度大，形成的坡面較粗糙；二是因為數(shù)值模擬和室內(nèi)試驗模擬月壤的黏聚力相差較大，黏聚力小的容易光滑，黏聚力大的容易形成裂紋。根據(jù)Bui等[10]對模擬月壤90°滑坡試驗結果，出現(xiàn)黏性土坡的微裂縫并且最終坡度也較大，可歸結于其模擬的月壤含有較高的黏聚力2 kPa。而此處模擬的月壤黏聚力0.2 kPa，不足以形成高坡角微裂縫。以上這些差異都體現(xiàn)了月壤因月面獨特環(huán)境下而產(chǎn)生黏聚力，而這種黏聚力會直接影響破壞形態(tài)，決定著月壤偏向砂土還是偏向黏性土破壞。

圖4 數(shù)值模型邊界條件及切坡示意圖(單位：m)

圖5 1/6重力場下地基初始應力水平及孔隙比分布

圖7給出數(shù)值模擬和室內(nèi)試驗下5種傾角保持坡后體積不變時，滑坡前后形態(tài)示意圖。

從圖7(a)中可知盡管5種坡體滑坡前形態(tài)不一致，但滑坡完成后的形態(tài)基本相似；無論在室內(nèi)試驗還是數(shù)值模擬中，月壤滑坡的過程均為土體顆粒圍繞坡面中點旋轉再一次達到平衡的過程，如圖7(a)、7(b)所示。

圖6 滑坡過程形態(tài)對比分析

4.2 滑坡過程流滑距離和最終傾角

為了減少滑坡帶來的危害，對滑坡過程常見的指標進行預測是非常有必要的。根據(jù)Scheidgger[15]采用的滑坡體積、垂直滑落高度和水平滑移距離三者存在一定的聯(lián)系來預測滑動距離的研究思路，本文通過控制坡后體積保持一定，簡化分析滑坡距離與邊坡傾角之間的關系，如圖8(a)所示。流滑距離為坡趾處到滑坡舌的距離，為了便于理解及使用，將流滑距離和坡高的比值定義為歸一化流滑距離Smax。隨著邊坡傾角的增大，歸一化流滑距離也不斷增大，這種增大趨勢可以通過一條開口向下的拋物線進行擬合：

Smax=-0.01395α2+3.0455α-70.088

(1)

通過對坡高進行歸一化處理，并且代入臨界值為50%時，得到一條實用的結論：針對此類月壤，當坡角高于51.7°時，沖擊到最遠的土體的歸一化值大于50%，即超過了邊坡高度的一半；對于最危險的90°坡體其流滑距離也小于坡高。圖8(b)給出了初始切坡傾角與坡體最終傾角的關系，由圖8(b)可看出，初始坡角對最終穩(wěn)定的坡度幾乎沒有影響，最終坡度在29°附近，以45°傾角時呈現(xiàn)微弱的峰值。根據(jù)文獻[12]此類月壤在低圍壓的峰值內(nèi)摩擦角為42.6°、殘余內(nèi)摩擦角為20.2°。5種傾角下滑坡失穩(wěn)后的穩(wěn)定邊坡的角度都小于峰值內(nèi)摩擦角，而大于殘余內(nèi)摩擦角，且基本穩(wěn)定在休止角30.6°[16]。

圖7 滑坡前后形態(tài)對比分析

4.3 滑坡過程最大速度分析

滑坡的速度直接決定滑坡時長，亦是人類對抗滑坡災害反應的時間，如快速滑坡比緩慢滑坡更易造成對生命財產(chǎn)的破壞，因此研究滑坡過程的速度變化規(guī)律具有重要的工程意義。圖9(a)給出了5種切坡傾角下最大速度Vmax隨時間變化的規(guī)律。隨著時間的增長，最大速度逐漸增長并且在0.2～0.8 s內(nèi)相繼達到峰值，之后隨著時間的增長最大速度逐漸接近零。這與Laouafa等[17]從速度的角度上描述滑坡過程一致，認為滑坡的本質(zhì)其實就是受到外力的觸發(fā)由“固態(tài)”變?yōu)椤耙簯B(tài)”，最后受到外部阻力又從“液態(tài)”變?yōu)椤肮虘B(tài)”的過程。根據(jù)最大速度變化規(guī)律可以將不同坡度的滑坡分為類崩塌型和流滑型，類崩塌型的以60°、75°和90°為典型，出現(xiàn)較大時間跨度下的高速度場區(qū)，類似于突然間的崩塌；流滑型以30°和45°為典型，整個滑坡過程中不會有明顯的特大速度，而是始終以穩(wěn)定的低速度進行流動，更似顆粒的流動滑移。若將5種傾角在滑坡過程中出現(xiàn)的最大速度單獨提取出來，會發(fā)現(xiàn)最大速度基本上隨著切坡傾角呈現(xiàn)拋物線式的增大，如圖9(b)所示，這條曲線也可以通過二次多項式進行擬合：

Vmax=0.0018α2-0.0338α+1.2175

(2)

但這種增大趨勢與歸一化流滑距離增長趨勢不一樣，最大速度增長速率隨著坡度增大不斷加劇，流滑距離增大速率隨著坡度增大逐漸放緩。國際地質(zhì)科學聯(lián)合會[18]規(guī)定對于3 m/min～5 m/s之間的滑坡可以定位為快速滑坡,而大于5 m/s的可以認為是極速滑坡，而極速滑坡易造成極大的生命財產(chǎn)損失；將極速滑坡的下限值5 m/s代入擬合公式(2)得到α=56°，這表明當邊坡坡角高度超過56°時，有很大的可能出現(xiàn)月壤的極速滑坡。所以從滑坡速度方面設計月壤坡體時，盡管月壤有較強的自穩(wěn)能力，仍應該減少高陡坡體。

4.4 滑坡過程S-t曲線分析

斜坡從變形的產(chǎn)生到最終失穩(wěn)破壞，累計位移-時間曲線(以下簡稱S-t曲線)，一般都會經(jīng)歷起滑階段、滑動階段及穩(wěn)定階段。前人根據(jù)S-t曲線斜率的特點將斜坡變形演化分成不同的階段，并據(jù)此建立滑坡預警判據(jù)[19-20]，所以研究不同傾角下的斜坡S-t曲線具有重要的意義。本次模擬中，選取具有代表性的3個特征點動態(tài)追蹤其累計位移，3個特征點分別布置在5種邊坡的坡趾處、坡面中點處及坡肩處，共計15個，主要研究滑坡過程邊坡傾角及監(jiān)測點布置對S-t曲線的影響，如圖10所示。

由圖10可知所有的特征值點的S-t曲線都分為了起滑階段、滑動階段及穩(wěn)定階段。邊坡傾角對S-t曲線的影響在不同布置區(qū)域呈現(xiàn)不同的規(guī)律：

圖8 月壤滑坡的流滑距離和最終坡角的預測

圖9 滑坡過程最大速度隨時間和邊坡傾角的變化曲線

圖10 特征點累計位移曲線

(1)坡肩處(圖10(a))，邊坡傾角越大，則同一時間監(jiān)測點的累計位移越大，表明坡肩處土體運動的距離和切坡傾角是正相關的。與此同時，隨著切坡傾角的增大，坡體再次達到穩(wěn)定階段的時間也在逐步延長。

(2)坡面中點處(圖10(b))，坡面中點累計位移基本上仍隨著切坡傾角的增大而增加，但90°例外。對于這種反?，F(xiàn)象，通過輸出全過程的變形圖，發(fā)現(xiàn)坡面中部特征點的位移受到坡頂土體垂直下落的影響，阻礙了這部分顆粒向前運動，因此累計位移減小。

(3)破趾處，中等坡度(如45°)的位移增長得最快，坡度過陡(如90°)或者過緩(如30°)的坡體其S-t曲線都會出現(xiàn)累計位移增長較慢。這主要是因為30°坡趾處由于坡面較緩，滑坡呈現(xiàn)輕微流動，對坡趾處監(jiān)測點的土體提供的推力較?。?0°坡趾處在滑坡過程中，坡肩處顆粒垂直落在坡趾，使得坡趾處的顆粒壓實變密，坡后的土體不但沒能提供坡趾處土體足夠的推力，反而成了其阻力，這部分與坡面中點特征點極為相似，只是阻力作用更為顯著。而介于這兩個極端的坡體即中間坡度的坡體能夠由坡后土體向前推動，使得其位移最大。綜上所述，采用S-t曲線做滑坡預警時，需要注意監(jiān)測點布置位置差異帶來對監(jiān)測數(shù)據(jù)的影響。若需要根據(jù)坡角信息特征及簡單時間段估算滑坡過程處于何種階段，最好將監(jiān)測點布置在坡肩頂部，此時所得信息較為可靠。

5 結 論

本文通過引入考慮范德華力及抗轉動作用的月壤微觀接觸模型，采用離散單元法模擬了月壤切坡試驗，從多個角度分析了不同切坡傾角對滑坡災害的影響，并且建立了滑坡過程中流滑距離、最終傾角及最大速度與切坡傾角的大致關系。得到主要結論如下：

(1)通過與試驗對比，在考慮月壤之間的微距力的情況下，坡體具有一定的黏性，容易形成滑坡凸體及可能的裂縫，危害較大。初始切坡傾角對最終傾角幾乎沒有影響；但發(fā)現(xiàn)最大流滑距離隨著初始坡角的增大而增大，當初始坡角超過51.7°時，最大流滑距離會超過坡高的一半。

(2)從工程危害角度分析了5種傾角下滑坡過程最大速度的規(guī)律，將其分為類崩塌型和流滑型，土體的最大速度隨初始坡角的增大呈現(xiàn)拋物線式增長，通過拋物線擬合發(fā)現(xiàn)，當坡角低于56°時，可有效避免極速滑坡。

(3)通過對5種坡角的坡肩處、坡面中部及坡趾處共15個特征點的S-t曲線分析，發(fā)現(xiàn)S-t曲線都經(jīng)歷了起滑階段、滑動階段及穩(wěn)定階段，但只有當監(jiān)測點布置在坡肩處時累計位移才會隨著角度的增大而增大。若需要根據(jù)坡角信息特征及簡單時間段估算滑坡過程處于何種階段，最好將監(jiān)測點布置在坡肩頂部。

參考文獻：

[1] 歐陽自遠. 月球科學概論[M]. 北京: 中國宇航出版社, 2005.

[2] 歐陽自遠. 我國月球探測的總體科學目標與發(fā)展戰(zhàn)略[J]. 地球科學進展, 2004,19(3):351-358.

[3] 周增坡, 程維明, 萬 叢,等.月球正面撞擊坑的空間分布特征分析[J]. 地球信息科學學報, 2012,14(5):618-626.

[4] 喬 樂, 劉小倩, 趙健楠,等. 月球雨海地區(qū)三個著陸點的地質(zhì)特征對比研究[J]. 中國科學(物理學 力學 天文學),2016,46(2):31-45.

[5] NAKASHIMA H, SHIOJI Y, KOBAYASHI T, et al. Determining the angle of repose of sand under low-gravity conditions using discrete element method [J]. Journal of Terramechanics, 2011,48(1):17-26.

[6] HORGAN B H N, BELL J F. Seasonally active slipface avalanches in the north polar sand sea of Mars: Evidence for a wind-related origin [J]. Geophysical Research Letters, 2015,39(9):9201-9207.

[7] CHANG C S, HICHER P Y. A constitutive model for granular materials with surface energy forces[J]. Journal of Aerospace Engineering, 2014,22(1):43-52.

[8] PERKO H A, NELSON J D, SADEH W Z. Surface cleanliness effect on lunar soil shear strength[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2001,127(4):371-383.

[9] ZHENG Hu, HUANG Yu. Model tests on flow slide of lunar regolith simulant[J]. Environmental Earth Sciences, 2014,73(8):1-7.

[10] BUI H H, KOBAYASHI T, FUKAGAWA R, et al. Numerical and experimental studies of gravity effect on the mechanism of lunar excavations[J]. Journal of Terramechanics, 2009,46(3):115-124.

[11] 鄭 敏, 蔣明鏡, 申志福.簡化接觸模型的月壤離散元數(shù)值分析[J].巖土力學, 2011,31(S1):766-771.

[12] JIANG Mingjing, SHEN Zhifu, THORNTON C. Microscopic contact model of lunar regolith for high efficiency discrete element analyses[J]. Computers and Geotechnics, 2013,54(10):104-116.

[13] MITCHELL J K, BROMWELL L G, DAVID C W, et al. Soil mechanical properties at the Apollo 14 site[J]. Journal of Geophysical Research, 1972,77(29):5641-5664.

[14] JIANG Mingjing, KONRAD J M, LEROUEIL S. An efficient technique for generating homogeneous specimens for DEM studies [J]. Computers and Geotechnics, 2003,30(7):579-597.

[15] SCHEIDEGGER A E. On the prediction of the reach and velocity of catastrophic landslides[J]. Rock Mechanics, 1973,5(4):231-236.

[16] 蔣明鏡, 奚邦祿, 申志福,等. 不同重力下月壤水平推剪阻力離散元數(shù)值分析[J]. 巖土工程學報, 2015，37(7):1300-1306.

[17]LAOUAFA F, DARVE F. Modelling of slope failure by a material instability mechanism[J]. Computers and Geotechnics, 2002,29(4):301-325.

[18] International Union of Geological Sciences Working Group on Landslides. A suggested method for describing the rate of movement of a landslide[J]. Bulletin of the International Association of Engineering Geology, 1995,52(1):75-78.

[19] 許 強, 曾裕平, 錢江澎,等.一種改進的切線角及對應的滑坡預警判據(jù)[J]. 地質(zhì)通報,2009,28(4):501-505.

[20] 許 強,湯明高,徐開祥,等.滑坡時空演化規(guī)律及預警預報研究[J].巖石力學與工程學報,2008,27(6):1104-1112.