黃土沉積過程及微結(jié)構(gòu)模型的非連續(xù)變形分析*

張 杰 李 萍 李同錄 張常亮 喬志甜 李 強(qiáng) 沈 偉

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

0 引 言

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

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

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

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

1.1 黃土顆粒形態(tài)的確定

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

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

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

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

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

1.2 黃土顆粒群的生成

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

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

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

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

(4)由于顆粒的位置和方向是隨機(jī)的，隨著顆粒增多，新顆粒和已有顆?？赡馨l(fā)生重疊，當(dāng)發(fā)生重疊時，放棄該顆粒，重新抽樣。

(5)經(jīng)反復(fù)測試，當(dāng)顆粒連續(xù)15～20次與已有顆粒重疊時，顆粒已經(jīng)均勻地填滿指定空間，則結(jié)束運(yùn)算。

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

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

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

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

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

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

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

表 1 模型物理力學(xué)參數(shù)Table 1 Physical and mechanical parameters of model

初始數(shù)值模型建立后，同時釋放所有顆粒，顆粒在下落中會相互碰撞和摩擦，顆粒群在重力作用下下降并在模型盒中堆積。最終堆積穩(wěn)定，形成的黃土初始微觀結(jié)構(gòu)模型如圖 5所示?？梢钥闯?，生成的模型中，底部較密實，向上變疏松，中上部有架空的大孔隙，和實際沉積的黃土結(jié)構(gòu)有相似性。

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

圖 6 黃土中的T型接觸結(jié)構(gòu)Fig. 6 T-type contact structure in loess

在圖 5中，可以明顯看到大、中、小3種類型的孔隙。其中，大孔隙半徑>0.016mm，主要由骨架顆粒相互支撐形成，孔徑大于骨架顆粒粒徑，這類孔隙較多地存在于黃土的初始結(jié)構(gòu)中，黃土濕陷時，這類孔隙首先破壞，顆?？梢缘袈淦渲校?yōu)橹锌紫痘蛐】紫?。中孔隙半徑?.004～0.016mm，主要由黃土顆粒相互穿插嵌套形成，孔徑比骨架顆粒粒徑小，在黃土受荷時，支架顆?？梢砸苿雍娃D(zhuǎn)動，使得這類孔隙變小。小孔隙的孔徑遠(yuǎn)遠(yuǎn)小于骨架顆粒的粒徑，為0.001～0.004mm，體積小，數(shù)量多，呈縫隙狀或者似三角形，這類孔隙十分穩(wěn)定，不足以變形。

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

(1)

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

3 壓縮試驗?zāi)M

利用DDA的算法，通過在受荷位置施加點(diǎn)荷載的方式給生成的微結(jié)構(gòu)模型結(jié)構(gòu)加載，通過指定測量點(diǎn)的位移變化確定變形參數(shù)并判定是否穩(wěn)定。通常當(dāng)測量點(diǎn)的位移在兩個時間步的差值為0.0001時，可認(rèn)為在該級荷載下達(dá)到穩(wěn)定，施加下一級的荷載(郭龍驍?shù)龋?2017)。對圖 5中的黃土初始結(jié)構(gòu)先施加5kPa的正應(yīng)力壓縮，變形穩(wěn)定后黃土的亞穩(wěn)態(tài)(metastable)結(jié)構(gòu)二維孔隙比為0.49(圖 7)，對應(yīng)的三維孔隙比為1.13。在此基礎(chǔ)上，設(shè)置了4個測量點(diǎn)和荷載點(diǎn)，對其施加10kPa、25kPa、50kPa、100kPa、200kPa、400kPa、800kPa的荷載進(jìn)行模擬。

圖 7 黃土的亞穩(wěn)態(tài)結(jié)構(gòu)Fig. 7 Metastable structure of loess

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

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

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

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

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

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

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

(2)

式中：ρ(r)為粒子徑向分布的數(shù)密度；ρa(bǔ)為整個系統(tǒng)的平均粒子數(shù)密度。

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

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

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

4 與物理模型試驗的對比

4.1 試驗設(shè)計

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

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

4.2 結(jié)果對比與分析

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

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

5 結(jié) 論

本文提出模擬黃土沉積過程和構(gòu)建黃土初始結(jié)構(gòu)的方法，并對黃土的壓縮試驗進(jìn)行模擬，得出以下結(jié)論：

(1)生成黃土的初始結(jié)構(gòu)，對該結(jié)構(gòu)中的孔隙及接觸結(jié)構(gòu)進(jìn)行分析，發(fā)現(xiàn)所生成結(jié)構(gòu)中主要存在大、中、小3種孔隙； 存在臺階(staircase)、堆疊(stack)、點(diǎn)接觸(point contact)、T型4類主要的接觸結(jié)構(gòu)。

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

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