結構性黃土壓縮特性的微觀非連續(xù)變形分析

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

(1.長安大學 地質(zhì)工程與測繪學院，西安 710054；2.長安大學 黃土高原水循環(huán)與地質(zhì)環(huán)境教育部野外觀測研究站，甘肅 正寧 745399)

1 研究背景

黃土是一種成分均一、多孔隙、水敏性強的風力搬運堆積物，其在沉積過程中形成了特殊的大孔隙骨架結構，并經(jīng)過長期的物理化學作用逐漸形成了膠結，具有顯著的抗變形破壞能力，表現(xiàn)出較強的結構性[1-2]。

黃土的結構性對其力學特性具有很大影響，是其強度、變形的內(nèi)在決定因素[3]，被認為是“21世紀土力學的核心問題”[4]。土的微觀結構研究是土結構性研究的開端，充分認識土的微觀結構是揭示土結構性的重要途徑。因此，研究黃土的結構性需從微觀入手，對黃土在變形過程中的微觀結構變化進行分析。目前，研究黃土微觀結構變化主要有2種方法：一種是室內(nèi)試驗法，另一種是數(shù)值模擬法。

在室內(nèi)試驗方面，蒲毅彬等[5]對原狀黃土在側限壓縮、三軸壓縮、浸水及加荷-浸水過程中分別進行了CT掃描，分析了土樣在各個試驗過程中細觀結構的變化特征；雷勝友等[6]進行了原狀黃土三軸剪切、浸水試驗過程的CT掃描，分析了原狀黃土在各個試驗過程中的微觀結構變化；朱元青等[7]利用CT掃描試樣斷面獲得結構性黃土試樣的孔隙、孔洞變化，裂隙開展等微觀信息，提出了基于CT數(shù)均值的結構性參數(shù)，揭示了黃土在加載和濕陷過程中的結構變化；谷天峰等[8]利用掃描電鏡獲取了黃土試樣的微觀結構圖像，分析了黃土在動應力作用下微觀結構的變化；陳陽等[9]借助掃描電子顯微鏡和光學數(shù)碼顯微鏡獲取了黃土試樣的微結構圖像，分析了黃土濕陷前后土樣中大、中、小孔隙和微孔隙數(shù)量的變化，從微觀角度分析了黃土濕陷的成因機理。

在數(shù)值模擬方面，非連續(xù)變形分析方法(DDA)能夠獲得數(shù)值試驗過程中的顆粒、孔隙、接觸、應力等微觀信息，是一種很有發(fā)展前景的數(shù)值分析方法。郭培璽等[10-11]采用DDA研究粗粒料的力學特性，給出了反映顆粒相互作用的若干組構要素的分布情況；張國新等[12]將土顆粒簡化為相同大小的11邊形，采用DDA模擬了土的平面應變試驗，分析了土的彈塑性、剪脹性、應變軟化、卸載-再加載時的滯回圈、卸載體縮以及剪切帶形成的微觀機理；郭龍驍?shù)萚13]應用DDA模擬黃土的單向固結試驗，說明采用該方法分析黃土的微觀變形特性是可行的；Guo等[14]在DDA的基礎上，考慮了非飽和土中的毛細作用，計算了不同含水率下土體中的毛細水分布及基質(zhì)吸力，模擬了非飽和土的土水特征曲線。

綜上可知，前人在黃土微觀結構方面進行了大量研究，但仍存在一些問題。在室內(nèi)試驗方面，目前室內(nèi)測試技術無法對粒間膠結作用進行定量研究，難以測試顆粒和孔隙結構的變化，無法分析顆粒的運動過程。在數(shù)值模擬方面，研究結構性黃土的變形行為時，尚未有人將黃土在黃土化(loessification)過程中形成的黏粒膠結作用考慮進去?；谝陨蠁栴}，針對結構性黃土的大孔隙和膠結特性，本文在借鑒已有成果的基礎上對DDA程序進行擴展，考慮了黃土中的黏粒膠結作用，并利用Monte Carlo法和DDA模擬黃土的沉積過程，建立結構性黃土微觀結構模型。在此基礎上模擬黃土的壓縮試驗，分析了其在壓縮過程中微觀結構的變化規(guī)律。

2 黏粒膠結接觸模型簡介

黃土中含有一定量的黏粒，黏粒具有很大的表面能，在含水率<6.0%時被吸附的水分子為結合水，此時黏粒表面吸附著完全定向排列的水分子層和水化陽離子，構成黏粒膠團(如圖1所示)，從而具有黏粒膠結力，并且該膠結力的大小隨著含水量的變化而變化。黏粒膠結使黃土保持固體狀態(tài)，具有相當?shù)膹姸萚15]，這部分力主要以聚集體包膜(又稱為黏粒膠膜)的形式作用在黏粒和吸附結合水組成的集合體與粗顆粒的接觸部位，如圖2所示[16]。本文對原DDA算法進行了拓展，在粒間接觸處施加了黏粒膠結作用，黏粒膠結是一種長程膠結力[17]。先用DDA識別出顆粒接觸點，然后根據(jù)水的總量將聚集體包膜以水膜的形式分散到這些接觸點上，最后計算出黏粒膠結力的大小施加于顆粒之間。

圖1 黏粒膠團示意圖

圖2 黃土中黏粒膠結作用示意圖

圖3 圓心軌跡交匯算法流程

2.1 聚集體包膜的分配方式

聚集體包膜按照水膜的形式分配，采用圓心軌跡交匯算法[14]，其算法流程如圖3所示。下面以2個顆粒為例說明其分配方式。為了簡化計算，假設水-空氣界面的曲率半徑r在整個土體系統(tǒng)中是相同的。確定聚集體包膜分布的具體步驟如下：①以顆粒A為對象，在顆粒A邊上選取任意一點PA0為起始點，過該點作垂直于顆粒A邊的外法線，將外法線逆時針旋轉(zhuǎn)接觸角θ，在該線上找出點OA0，使PA0OA0長度為r；②使線段PA0OA0沿顆粒A表面逆時針移動一周，此時線段PA0OA0的端點OA0會在顆粒A外側形成一條軌跡線，該線平行于顆粒A表面，需要注意的是，由于顆粒的端點是顆粒邊界上不光滑的點，為了使顆粒A的圓心軌跡連續(xù)，將顆粒A的端點看成半徑很小的圓角;③在距顆粒A表面距離為2r的范圍內(nèi)尋找相鄰顆粒B，按照相同的方法作出相鄰顆粒B的圓心軌跡線，區(qū)別僅在于使PB0沿著顆粒的邊按順時針方向移動;④顆粒A與顆粒B的圓心軌跡線共有2個交點，分別為Oi和On，可根據(jù)給定的液面半徑r與交點Oi和On畫出液面，分別為C2和C1，C2與顆粒A的邊交于PAi，與顆粒B的邊交于PBi，C1與顆粒A的邊交于PAn，與顆粒B的邊交于PBn，其中凸液面C2為不合理液面，將其刪除，如圖4(a)所示；⑤互換顆粒A與顆粒B的位置，重復以上步驟，可得到另一側液面，如圖4(b)，將兩側合理液面組合，則確定了黏粒膠膜的分布區(qū)。

圖4 黃土顆粒圓心軌跡線及可能彎液面

2.2 黏粒膠結力的計算

黏粒膠結力的大小采用Kevin方程(式(1))計算，首先利用Kevin方程計算出黏粒之間的吸力ψ，再用ψ乘以膠結部分的斷面長度L獲得膠結力的大小，如圖5所示。

(1)

式中:ψ表示土中的吸力(kPa)；R表示通用氣體常數(shù)，值為8.314 32 J/(mol·K)；T為絕對溫度，T=273.16+t，t為攝氏溫度；ρw為水的密度；wv為水蒸氣的摩爾質(zhì)量，值為18.016 g/mol；RH為相對濕度。

圖5 黏粒膠結力的作用形式

開爾文方程中濕度RH可通過等溫吸附曲線來測定，選用陜西咸陽黃土塬北緣的剖面上L5的黃土樣的等溫吸附曲線，如圖6所示。

圖6 等溫吸附曲線

根據(jù)等溫吸附曲線，由Kevin方程可計算出吸力與含水率之間的關系，如圖7所示。

圖7 吸力與體積含水率關系曲線

用計算出的吸力乘以膠結部分的斷面長度L即得膠結力的大小。由于實驗室中對于天然黃土樣一般在105 ℃下進行烘干，此時黃土中仍存在部分弱結合水和強結合水[18]。為獲得在105 ℃下烘干后黃土的體積含水率，本文分別測定了陜西咸陽黃土塬北緣的剖面上L5層黃土樣天然含水率下的初始重量、105 ℃下烘干以及350 ℃下烘干的質(zhì)量，分別為98.520 5、96.469 6、94.383 3 g。然后通過式(2)計算出在105 ℃下黃土樣的質(zhì)量含水率ω為2.17%，通過式(3)換算成體積含水率為3.29%(干密度ρd=1.52 g/cm3)。

(2)

式中:mt為任意溫度下失去的水的質(zhì)量；M0為試驗結束后的樣品質(zhì)量。

(3)

式中：ρd為試樣的干密度；ρw為水的密度；ω為試樣的質(zhì)量含水率。

為與天然原狀黃土烘干后的試驗結果進行對比，本模型使用體積含水率為3.29%時的濕度0.78作為模型的參數(shù)，此時吸力值為3.36×104kPa。

3 結構性黃土微觀模型的建立

黃土具有典型的微觀架空結構，其架空結構是由于粉塵降落、粉粒搭接而產(chǎn)生的[19]。建立黃土的微觀結構模型，首先要確定黃土骨架顆粒的大小和形狀，Rogers等[20]運用Monte Carlo理論研究了黃土顆粒的形狀和大小，結果表明在三維空間中黃土以扁平狀石英顆粒為主，顆粒長、寬、高遵循8∶5∶2的比例，其粒徑多集中在20～30 μm，如圖8(a)所示。Dibben等[21]研究了二維平面上顆粒的形狀和大小對黃土結構的影響，結果表明在二維平面上用長寬比為3∶1的矩形顆粒代表石英顆粒較為合理。因此，本文沿用Dibben等的結論，在二維平面上采用3∶1(長∶寬)的矩形塊體單元模擬黃土顆粒形狀。為確定顆粒大小，統(tǒng)計了400組陜西自陜北到關中馬蘭黃土的粒度分析結果，作出其平均粒徑d50的頻率分布，最終確定黃土顆粒的最優(yōu)勢粒徑為28 μm。為簡化計算采用等效面積法將模型中顆粒的尺寸等效成直徑28 μm的圓，計算得矩形顆粒的長度為42.9 μm，寬度為14.3 μm，如圖8(b)。

圖8 黃土顆粒的形狀

實際上黃土顆粒的形狀和大小各異，并且黃土中存在團粒結構，若按實際情況模擬，計算將變得十分復雜，目前的計算能力基本難以實現(xiàn)。本文先注重于算法的探討，因此將黃土骨架顆粒按大小和形狀相同處理。

確定顆粒形狀和大小后，利用Monte Carlo法在模型盒的正上方生成沉積前的黃土顆粒群，然后利用非連續(xù)變形分析方法(DDA)模擬顆粒的自由下落，該方法能夠模擬顆粒下落過程中的相互碰撞及摩擦，待顆粒初步下落穩(wěn)定后對其施加5 kPa的正應力使其處于亞穩(wěn)態(tài)，由此構建黃土微觀結構模型，如圖9所示。微觀結構模型的建立過程詳見文獻[22]。圖9中模型的二維孔隙比為0.49，二維孔隙比與三維孔隙比的換算可用公式(4)確定(Hoomans 等, 1996)[23]，經(jīng)計算三維孔隙比為1.13。

(4)

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

圖9 結構性黃土微觀結構模型

本文按《土工試驗方法標準》(GB/T 237—1999)[24]測得陜西咸陽黃土塬北緣剖面上所取原狀黃土試樣孔隙比集中在1.050～1.120，少量試樣孔隙比達到1.30。由此可看出通過DDA數(shù)值模擬建立的黃土微觀結構模型的孔隙比和現(xiàn)實物理力學試驗中黃土的孔隙比接近。

4 結構性黃土壓縮試驗模擬

利用擴展的DDA的算法，在受荷板上施加若干集中荷載對微結構模型加載，通過指定測量點的位移變化確定變形參數(shù)并判定模型是否穩(wěn)定。當測量點的位移在2個時間步的差值<0.000 1 μm時，可認為模型在該級荷載下穩(wěn)定，開始施加下一級的荷載。對圖9中的黃土微觀結構模型設置了4個測量點和加載點，分別在無黏粒膠結和有黏粒膠結2種情況下對施加12.5、25、50、100、150、200、300、400、800、1 600、2 000 kPa的荷載進行模擬。由于試樣的初始孔隙比e0、高度H恒定不變，只需讀取試樣在壓縮過程中的變形量s，用式(5)計算試樣在不同應力下的孔隙比e，轉(zhuǎn)化為三維孔隙比后獲得的曲線如圖10所示。從圖10可以看出膠結試樣模型曲線變化量明顯小于不含膠結試樣模型，說明黏粒膠結力的存在顯著增強了黃土的強度，減小了黃土的壓縮性。

(5)

圖10 數(shù)值模型壓縮曲線

為了與實際情況對比，選取孔隙比大約1.130的3組原狀黃土試樣，將其置于105 ℃烘箱中烘干8 h。再分別在12.5、25、50、100、200、300、400、800、1 600、2 000 kPa的法向應力下壓縮，將3組平行試驗結果曲線與數(shù)值模擬曲線進行對比，如圖11所示?？梢钥闯鰯?shù)值模擬結果與室內(nèi)試驗結果展現(xiàn)了相同的趨勢，模擬結果可以較好地反映黃土的壓縮特性。需要指出的是，本文關于結構性黃土壓縮試驗，目前旨在從算法上探究其可行性，還未能實現(xiàn)對真實試樣的模擬，達到與試驗結果完全相吻合。

圖11 室內(nèi)試驗與數(shù)值模擬的壓縮曲線對比

5 結構性黃土顆粒運動分析

黃土的宏觀力學行為與其微觀結構變形密切相關，其在外力作用下發(fā)生的應變在微觀上是由于顆粒的運動引起的。為分析壓縮過程中顆粒的運動情況，提取各荷載下黃土的壓縮圖像以及對應的數(shù)據(jù)，其壓縮完成后圖像如圖12所示。從圖12可以看出，與膠結試樣模型相比，不含膠結試樣模型壓縮完成時表現(xiàn)出更大的壓縮量，更多顆?；氪罂紫吨?，大孔隙的數(shù)目明顯減少，整體結構看起來更加緊密。

圖12 壓縮完成后的模型

建立以樣盒左下角點為原點的平面直角坐標系，選擇長方形單元體質(zhì)心坐標作為顆粒的坐標，連接初始狀態(tài)與壓縮完成后顆粒的質(zhì)心坐標，標出方向，可得到顆粒的位移矢量，如圖13所示。

圖13 顆粒位移矢量示意圖

圖14 顆粒位移矢量

圖14分別為膠結試樣和不含膠結試樣顆粒的位移矢量。從圖14可以明顯看出膠結試樣顆粒位移比不含膠結試樣顆粒位移小，但整體規(guī)律相似。每個顆粒位移的大小和方向不盡相同，但主要以豎向位移為主。在壓縮過程中，不同區(qū)域的顆粒位移存在較大差異，上部顆粒位移大，越往下部顆粒的位移越小，且兩側及中部顆粒較密集的部位也出現(xiàn)了較大的位移，這是因為該部位密集，使得顆粒被擠向了孔隙。值得注意的是，在壓縮過程中粒間的相對運動，使得少量顆粒出現(xiàn)了上浮現(xiàn)象，但上浮幅度較小，模型總體高度仍在下降。

圖15 顆粒豎向位移與水平位移示意圖

5.1 平動結果分析

對膠結試樣進行平動

分析，計算顆粒的豎向位移Δy與水平位移Δx，如圖15所示。豎向位移Δy和水平位移Δx的公式為：

式中:(x1,y1)為顆粒運動起點質(zhì)心坐標;(x2,y2)為顆粒運動終點的質(zhì)心坐標。

圖16為壓縮完成時不同位置(x,y)顆粒的平動(Δx,Δy)成果。從圖16(a)可以看出顆粒的豎向位移(Δy)隨著高度(y)的增大而增大，但同一高度的不同顆粒豎向位移是不同的，這說明了同一高度處左右相鄰的顆粒有錯動，顆粒位移矢量圖也反映了這點。從圖16(b)可以看出不同高度處顆粒的水平位移呈現(xiàn)出對稱式分布的特征，300 μm高度以上的部位其水平位移相對較大，最大可達20 μm，中下部顆粒的水平位移很小且水平位移大致相同，主要集中于-5～5 μm。從圖16(c)可以看出不同水平位置處顆粒的水平位移也呈現(xiàn)出對稱式分布的特征，但同一水平位置處顆粒的水平位移差異較大。圖16(b)、圖16(c)總體說明黃土顆粒在水平方向上的運動是隨機的，即在顆粒左、右兩側平均分配。

圖16 顆粒位移散點圖

統(tǒng)計各個加載階段位移L，如表1所示。L的最小值為7.5×10-7μm，最大值為11.14 μm，最大值與最小值相差非常大，說明黃土在壓縮過程中顆粒運動過程復雜, 對顆粒的約束力或限制力大小不同, 使得黃土顆粒的位移沒有明顯規(guī)律，表1中的標準差數(shù)據(jù)也能反映這一點。

表1 各加載階段位移結果

圖17 顆粒的長軸方向

5.2 轉(zhuǎn)動結果分析

顆粒在荷載作用下會發(fā)生轉(zhuǎn)動，對膠結試樣進行顆粒的轉(zhuǎn)動結果分析，用顆粒的長軸方向作為指標來描述顆粒的轉(zhuǎn)動行為，顆粒的長軸方向是指顆粒長軸與水平方向的夾角，如圖17所示，用α表示(0°<α<180°)。為了更好地反映加載完成后顆粒的轉(zhuǎn)動規(guī)律，繪制出壓縮完成時不同位置(x,y)顆粒的轉(zhuǎn)動(ΔФ)成果,如圖18所示。從圖18(a)、圖18(b)可以看出顆粒的轉(zhuǎn)動角度(ΔФ) 呈現(xiàn)對稱式分布，主要集中于-30°～30°，最大轉(zhuǎn)動角度可達到±180°，整個試樣發(fā)生逆時針(正值)和順時針(負值) 轉(zhuǎn)動的顆粒數(shù)及轉(zhuǎn)動量的分布范圍大致相同，不同高度以及不同水平位置處顆粒的轉(zhuǎn)動角度差異不大，說明顆粒的轉(zhuǎn)動方向僅與顆粒的長軸方向的隨機分布有關，與顆粒的平動關系不大。

圖18 顆粒轉(zhuǎn)動角度散點圖

為了反映在不同加載階段顆粒的轉(zhuǎn)動規(guī)律，將顆粒長軸與水平方向的夾角α在[90°,180°]區(qū)間內(nèi)的角度換算為[-90°,0°]，得到新的轉(zhuǎn)動角度β，引入γ和Δ兩個參數(shù)來描述顆粒主軸方向的變化情況[25]。其中:

(8)

(9)

式中:β為α調(diào)整后顆粒的轉(zhuǎn)動角度(-90°≤α≤90°);N為顆粒個數(shù)。γ反映主軸方向的大小，即被“壓倒”的程度，而Δ反映全部顆粒主軸方向的統(tǒng)一性。Δ介于0～100，當Δ=0時，說明顆粒主軸方向完全趨于隨機分布;Δ=100時，說明顆粒主軸方向完全一致。Δ與γ的計算結果如圖19所示。

圖19 γ和Δ與荷載的關系

從圖19可以看出，隨著荷載的增加，γ逐漸減小，Δ逐漸增加。γ減小說明了隨著荷載的不斷增大，更多的顆粒趨向于水平排列來維持自身的穩(wěn)定；Δ增加說明隨著荷載的不斷增大，顆粒的主軸方向趨于一致,顆粒排列趨于規(guī)則。此外，從圖19曲線中可明顯看出一個拐點，該拐點是由于試樣結構在上一級荷載(300 kPa)作用下已經(jīng)達到相對穩(wěn)定的狀態(tài)，且荷載梯度(100 kPa)較小，不足以打破原有的平衡狀態(tài)，使得試樣的γ和Δ變化不大，當對試樣施加下一級荷載時，由于荷載梯度(400 kPa)較大，打破了原有的平衡狀態(tài)，使γ和Δ急劇變化。

6 結 論

本文將黏粒膠結作用嵌入已有的DDA算法中，模擬結構性黃土的壓縮試驗，并對黃土在壓縮過程中的顆粒運動情況進行了分析，得出以下結論：

(1)提出了黏粒膠結接觸模型，首先用DDA識別出顆粒接觸點，然后根據(jù)水的總量將聚集體包膜以水膜的形式分散到這些接觸點上，最后計算出黏粒膠結力的大小施加于顆粒之間。

(2)建立結構性黃土的數(shù)值模型，模擬黃土的一維壓縮試驗，并與室內(nèi)試驗結果進行對比，論證了所建數(shù)值模型的合理性。

(3)對結構性黃土在一維壓縮過程中的顆粒的平動結果和轉(zhuǎn)動結果進行分析，從微觀上說明了黃土在壓縮過程中顆粒的運動規(guī)律。