黃土沉積過程及微結構模型的非連續變形分析*

張 杰 李 萍 李同錄 張常亮 喬志甜 李 強 沈 偉

(①長安大學地質工程與測繪學院， 西安 710064， 中國) (②黃土高原水循環與地質環境教育部野外觀測研究站， 正寧 745399， 中國) (③博洛尼亞大學， 生物、地球與環境科學系， 博洛尼亞 67-40126， 意大利)

0 引 言

黃土是一種成分均一、多孔隙、水敏性強的風力搬運堆積物，具有顯著的結構性(朱興華等， 2017)，其變形破壞特征實際是微觀結構發生破壞的宏觀表現，探究黃土微觀結構的破壞行為有利于揭示黃土宏觀變形破壞的根源(袁志輝等， 2018)。

黃土的微觀結構是組成黃土的原始物質在沉積后，于特定的自然環境下，經過一系列物理化學作用形成的(方祥位等， 2013)。探究黃土微觀結構的變化，需首先建立黃土原始物質沉積后形成的初始結構。目前建立黃土微觀結構主要有兩種方法，一種是數值模擬法，另一種是圖像處理法。數值模擬法通過Monte Carlo等方法在一定空間內填充土體顆粒，從而形成黃土的微觀結構。如Rogers et al. (1994)通過在一定范圍內賦予顆粒隨機角度和位置的辦法建立了黃土的初始結構，并研究了顆粒之間的膠結作用，認為前期膠結作用是黃土后期產生濕陷的重要原因； Assallay et al. (1997)在假定顆粒水平、角度不變、膠結力與顆粒重力比為1，即顆粒在接觸瞬間便粘合的前提下，使用Monte Carlo法隨機充填并構建了黃土的二維初始結構； Dibben et al. (1998a，1998b)進一步拓展了Assallay et al. (1997)的方法，假定了一個臨界膠結力，當該膠結力大于重力時，兩顆粒粘合，否則顆粒移動，據此建立了黃土的初始結構，并且研究了顆粒的大小和形狀對結構的影響。圖像處理法則通過掃描電鏡以及斷面掃描(CT)等技術對黃土進行斷面掃描獲得黃土的微觀結構圖像，然后統計求解相關結構參數，基于參數控制建立黃土微觀結構。如陳開圣等(2009)基于數字圖像處理技術，定性和定量分析了不同壓實度下黃土的微結構特征； 谷天峰等(2011)基于掃描電鏡圖像分析了Q3黃土的微結構變化； 陳陽等(2015)基于掃描電鏡圖像和光學顯微鏡圖像分析了黃土濕陷前后微結構的變化特征； Liu et al. (2015)基于掃描電鏡圖像建立了一種4層結構的黃土概念微結構模型； Li et al. (2019)基于掃描電鏡圖像建立了不同壓力和含水率下的微結構模型。綜上所述，數值模擬法在構建黃土微觀結構時，假定條件過于理想，沒有考慮黃土在形成初始結構時顆粒之間力的作用； 圖像處理法建立的結構模型只能做統計分析，無法進行力學模擬，不利于進一步分析結構變化。

基于以上研究，本文提出了一種模擬黃土沉積過程并構建黃土初始結構的方法，即先利用Monte Carlo法生成沉積前的黃土顆粒群，然后引入非連續變形分析方法(DDA)模擬黃土的沉積過程，該方法能夠模擬顆粒下落過程中的相互碰撞及摩擦，由此建立黃土的初始結構模型。最后，利用建立的結構模型模擬了黃土的壓縮試驗，通過與相同條件下物理模型試驗結果的對比，初步表明該方法是可行的。

1 基于Monte Carlo法的黃土顆粒群生成

1.1 黃土顆粒形態的確定

黃土顆粒的形狀是其結構的要素之一，它是決定黃土結構排列并引起濕陷的重要因素(Smalley et al.，2017)。Krinsley et al. (1973)認為微小的沉積石英顆粒是一種扁平狀的顆粒。在此基礎上，Rogers et al. (1993)應用Monte Carlo抽樣統計分析了黃土顆粒形狀和大小的分布模式，表明黃土中的石英顆粒主要是一種粒徑為30μm的扁平狀顆粒，且其長、寬、高之比約為8︰5︰2，如圖 1a所示(Domokos， 2010； Howarth， 2010； Assadi Langroudi et al.， 2014， 2018)。Dibben et al. (1998)研究了二維平面上顆粒的形狀和大小對黃土結構的影響，結果表明在二維平面上用長寬比為3︰1的矩形顆粒代表石英顆粒較為合理。本文根據Dibben et al. (1998)的研究結論，在模擬二維平面上的黃土沉積過程時，采用3︰1的矩形顆粒代表以石英、長石為主的粗骨架顆粒(圖 1b)。

圖 1 黃土顆粒的形狀Fig. 1 The shape of loess particles

為了確定黃土顆粒的最優勢粒徑，統計了400組陜西省自陜北到關中馬蘭黃土的粒度分析結果，統計出其平均粒徑d50的頻率分布，總共獲得400個平粒粒徑，作出平均粒徑頻率分布直方圖(圖 2)。可以看出黃土的粒徑主要集中于26～30μm，因此選取中位數28μm作為黃土的優勢粒徑。采用等效面積粒徑法確定二維平面上顆粒的粒徑，使矩形顆粒的面積與直徑為28μm的圓面積相等，可得矩形石英顆粒的等效長度為42.9μm，等效寬度為14.3μm。

圖 2 黃土顆粒平均粒徑頻率分布Fig. 2 Frequency distribution of average particle size of loess

實際上黃土顆粒的形狀和大小各異，若按實際形狀和大小模擬，計算將變得十分復雜，目前的計算能力基本難以實現。本文先注重于算法的探討，因此基于上述統計結果，將黃土顆粒按大小和形狀相同處理。

1.2 黃土顆粒群的生成

利用選定大小和形狀的二維顆粒，使用Monte Carlo法在預先設置的區域內，隨機生成沉積前的黃土顆粒群，本文給定一1.0mm×0.7mm的區域。具體步驟為：

(1)在已確定的垂直立面上的一個方形區域中，建立直角坐標系，水平為x坐標，垂直為y坐標。

(2)使用Monte Carlo法，抽取3個隨機數，將其分別換算為顆粒在指定區域的(x，y)值，及顆粒長軸與x軸的夾角α， 0≤α≤360°。

(3)只要新顆粒與已有顆粒發生一定的重疊，認為該顆粒生成失敗，同時賦予該顆粒一個新坐標。

(4)由于顆粒的位置和方向是隨機的，隨著顆粒增多，新顆粒和已有顆粒可能發生重疊，當發生重疊時，放棄該顆粒，重新抽樣。

(5)經反復測試，當顆粒連續15～20次與已有顆粒重疊時，顆粒已經均勻地填滿指定空間，則結束運算。

利用該方法生成的黃土顆粒群如圖 3所示。

圖 3 利用Monte Carlo法生成的黃土顆粒群Fig. 3 Loess particle group generated by Monte Carlo method

圖 4 黃土沉積DDA初始數值模型Fig. 4 Initial numerical model of DDA in loess deposition

2 基于DDA的顆粒群沉積過程模擬

首先建立用于存儲顆粒的模型盒，模型盒內部尺寸為0.7mm×0.5mm。利用Monte Carlo法隨機生成顆粒群置于模型盒上方，如圖 4所示。

圖 5 初始黃土微觀結構模型Fig. 5 Initial loess microstructure model a. 臺階結構， b. 堆疊結構， c. 點接觸結構， d. T型結構； ①. 大孔隙， ②. 中孔隙， ③. 小孔隙

由于黃土顆粒微小，在實際下落過程中，受到重力與空氣阻力的共同作用，并非加速下落，而是時快時慢，運動軌跡也非垂線，但其速度和軌跡對模型影響不大，因此在建模時不考慮速度和軌跡的影響，將顆粒的初速度設置為0。由于黃土顆粒的成分主要是長石和石英，類似于花崗巖(謝星等， 2013)，因此顆粒的彈性模量、泊松比參考花崗巖的參數取值，密度可通過室內試驗獲得，松散狀態下黃土的內摩擦角與其休止角相等，可取休止角。黃土顆粒及土樣盒的物理力學參數如表 1。相對于土顆粒，模型盒為剛體，因此模擬時使其彈性模量遠大于土體。

表 1 模型物理力學參數Table 1 Physical and mechanical parameters of model

初始數值模型建立后，同時釋放所有顆粒，顆粒在下落中會相互碰撞和摩擦，顆粒群在重力作用下下降并在模型盒中堆積。最終堆積穩定，形成的黃土初始微觀結構模型如圖 5所示。可以看出，生成的模型中，底部較密實，向上變疏松，中上部有架空的大孔隙，和實際沉積的黃土結構有相似性。

由圖 5生成的初始黃土微觀結構模型可以看出，顆粒間主要有4種接觸結構類型，分別為：a. 臺階結構(staircase)(Smalley et al.，1969)、b. 堆疊結構(stack)、c. 點接觸結構(point contact)(Assallay et al.，1997)、d. T型結構。其中，臺階結構(staircase)是黃土中最基本的一種接觸結構，在黃土中大量存在。堆疊結構(stack)中兩個顆粒之間的接觸面積比臺階結構(staircase)大，當兩顆粒之間的接觸面積大于50%時為堆疊結構，小于50%時為臺階結構，因此堆疊結構(stack)比臺階結構(staircase)穩定。點接觸結構非常不穩定，但是在初始結構中這種點接觸數量居多，在上覆土層壓力作用下，這種接觸結構首先被破壞，這種不穩定的點接觸是評價黃土塌陷的關鍵。T型結構也是一種不穩定的結構，在壓縮過程中，上部顆粒和下部顆粒均有可能倒塌，形成其他更穩定的接觸結構，在以往學者構建的模型中并沒有發現這種結構，實際黃土中可見到這種結構，如圖 6所示(Milodowski et al., 2015)。以上顆粒接觸結構往往不是單獨存在的，會與其他結構相互組合，組成更復雜的接觸結構。

圖 6 黃土中的T型接觸結構Fig. 6 T-type contact structure in loess

在圖 5中，可以明顯看到大、中、小3種類型的孔隙。其中，大孔隙半徑>0.016mm，主要由骨架顆粒相互支撐形成，孔徑大于骨架顆粒粒徑，這類孔隙較多地存在于黃土的初始結構中，黃土濕陷時，這類孔隙首先破壞，顆粒可以掉落其中，變為中孔隙或小孔隙。中孔隙半徑為0.004～0.016mm，主要由黃土顆粒相互穿插嵌套形成，孔徑比骨架顆粒粒徑小，在黃土受荷時，支架顆粒可以移動和轉動，使得這類孔隙變小。小孔隙的孔徑遠遠小于骨架顆粒的粒徑，為0.001～0.004mm，體積小，數量多，呈縫隙狀或者似三角形，這類孔隙十分穩定，不足以變形。

在圖 5所示模型的上方畫一條水平線，計算出初始結構的孔隙比為0.60。該模型在二維平面上建立，所生成的結構為二維結構，計算所得孔隙比與三維空間的孔隙比不同，二維要小于三維(黃磊等， 2016)。利用Hoomans et al. (1996)提出的等粒徑顆粒二維孔隙與三維孔隙的轉化式(1)，可實現二維與三維的轉化。經計算，所生成的初始結構三維孔隙比為1.36。

(1)

式中：n2D、n3D分別為土體的二維、三維孔隙率。

3 壓縮試驗模擬

利用DDA的算法，通過在受荷位置施加點荷載的方式給生成的微結構模型結構加載，通過指定測量點的位移變化確定變形參數并判定是否穩定。通常當測量點的位移在兩個時間步的差值為0.0001時，可認為在該級荷載下達到穩定，施加下一級的荷載(郭龍驍等， 2017)。對圖 5中的黃土初始結構先施加5kPa的正應力壓縮，變形穩定后黃土的亞穩態(metastable)結構二維孔隙比為0.49(圖 7)，對應的三維孔隙比為1.13。在此基礎上，設置了4個測量點和荷載點，對其施加10kPa、25kPa、50kPa、100kPa、200kPa、400kPa、800kPa的荷載進行模擬。

圖 7 黃土的亞穩態結構Fig. 7 Metastable structure of loess

當荷載小于50kPa時，由于荷載較小，結構變化不太明顯。故在這里只展示荷載大于50kPa的壓縮結果，如圖 8所示。從圖中可以明顯地看出，隨著荷載的增加，架空結構和大孔隙首先破壞，部分土顆粒進入大孔隙中，孔隙變小，當加載到800kPa時，大部分大孔隙基本閉合。

圖 8 各級荷載穩定后的模型Fig. 8 Model after load stabilization at each level a. 加載100kPa； b. 加載200kPa； c. 加載400kPa； d. 加載800kPa

圖9為黃土模型孔隙比轉化為三維后的e-lgp曲線，可以看出，該曲線和物理試驗獲得的e-lgp曲線非常相似。壓力為100kPa之前曲線平緩，其后曲線陡降，也與圖 5中大孔隙和架空結構的破壞相對應。

圖 9 數值模型壓縮曲線Fig. 9 Compression curve of the numerical model

圖 10 選定顆粒的徑向分布函數Fig. 10 The radial distribution function of selected particle

為了從微觀角度研究黃土初始結構模型在壓縮過程中的變化，引入徑向分布函數。徑向分布函數反映了以任意顆粒為中心，其周圍顆粒的對稱徑向分布狀態。在所生成的黃土模型中選擇某一顆粒，以該顆粒為中心畫出3個同心矩形區域，分別命名為內側(I)、中部(M)、外側(O)，如圖 5所示。每個區域的徑向分布函數的計算方法如式(2)。

g(r)=ρ(r)/ρa

(2)

式中：ρ(r)為粒子徑向分布的數密度；ρa為整個系統的平均粒子數密度。

在所建立初始模型中選擇3個不同部位的顆粒，顆粒1、2、3分別為大、中、小3類孔隙附近的顆粒，在不同荷載下，其徑向分布函數如圖 10所示。

從圖 10可以看出圖 5的黃土初始結構極為不穩定，在5kPa正應力作用下，顆粒1內側區域(I)的徑向分布密度顯著增加，中間區域(M)值顯著減少，說明中間區域(M)的顆粒滑向了內側區域(I)，使得最內側區域變得擁擠，在200kPa時顆粒1的內側區域(I)值達到最大，此后顆粒內側區域(I)值不再變大，說明此時顆粒1內側達到最擁擠狀態，隨后中間區域(M)值開始變大，并且內側區域(I)值出現減少，說明其內側區域(I)某些顆粒發生了側移，移向了中間區域(M)。顆粒2在5kPa的正應力作用下，內側區域(I)與中間區域(M)值均增加，但沒有顆粒1顯著，外側區域(O)值減小，說明外側區域(O)顆粒滑向了內側區域(I)和中間區域(M)。顆粒3在壓縮過程中徑向分布密度基本不變。綜上所述，黃土在壓縮過程中大孔隙和架空結構會首先被破壞，且破壞程度最大，而小孔隙基本不會變化。

圖 11 黃土沉積過程試驗方案示意圖Fig. 11 Schematic diagram of test scheme for loess deposition process

4 與物理模型試驗的對比

4.1 試驗設計

數值模擬只有在與天然物質的行為相匹配的情況下才有用，因此為了驗證數值模型的可行性，本文設計了一個物理模型試驗，如圖 11所示。土樣取自陜西咸陽黃土塬北緣剖面上的L5層黃土。先將土樣碾碎、烘干，測得其比重為2.71。將試樣置于75μm的篩中，在固結儀應力盒正上方一定高度輕輕晃動篩子，將土粒篩落在壓力盒中，待裝滿壓力盒時停止。本試驗將篩子水平置于壓力盒上方大約0.5m的位置。

待樣品裝滿壓力盒后，仔細平整表面，盡量避免擾動，此時測得樣品的孔隙比為1.83。接著在樣品上輕輕地放置一塊透水石并蓋好蓋子。固結儀的蓋子和透水石可提供1kPa的正應力，將正應力增加到5kPa，此時樣品形成了亞穩態(metastable)結構，孔隙比為1.44。隨后調整樣品，往壓力盒中篩入更多的樣品，并輕微地壓縮，制取兩組孔隙比約為1.14(與數值模擬形成的孔隙比大致相等)的樣品，將樣品放在一維固結儀中依次用5kPa、10kPa、25kPa、100kPa、200kPa、400kPa、800kPa的正應力壓縮。

4.2 結果對比與分析

繪制物理試驗與數值模擬的壓縮曲線e-lgp如圖 12所示，可以看出，曲線的趨勢相同，但兩者之間仍然存在一定的偏差，并且數值模擬形成的初始結構的孔隙比與模型試驗土樣的孔隙比存在一定差距，考慮主要是由于實際黃土顆粒的大小和形狀各異，而本文為了簡化模型用一種理想化的顆粒代替。后續將在考慮該因素的基礎上開展進一步研究。

圖 12 數值模擬與實際試驗對比Fig. 12 Comparison between numerical simulation and actual experiment

5 結 論

本文提出模擬黃土沉積過程和構建黃土初始結構的方法，并對黃土的壓縮試驗進行模擬，得出以下結論：

(1)生成黃土的初始結構，對該結構中的孔隙及接觸結構進行分析，發現所生成結構中主要存在大、中、小3種孔隙； 存在臺階(staircase)、堆疊(stack)、點接觸(point contact)、T型4類主要的接觸結構。

(2)模擬黃土的壓縮試驗，得到了壓縮曲線，并選取特定顆粒作出了徑向分布函數，從微觀角度說明黃土在壓縮過程中架空結構和大孔隙首先會被破壞，而小孔隙基本不會變化。

(3)設計物理模型試驗，與數值模擬結果進行對比，兩種壓縮曲線趨勢大致相同，表明本文方法是可行性，為黃土微觀力學的進一步研究提供了基礎。