路德春， 張世浩， 田 雨， 杜修力

(北京工業(yè)大學巖土與地下工程研究所， 北京 100124)

土是由顆粒構成的離散介質，外力能夠引起細觀上土顆粒排列方式的改變，同時顆粒排列也能影響宏觀上土體的變形、強度等力學特性，因此宏細觀2種尺度之間存在著密切聯(lián)系.

離散單元法是研究宏細觀聯(lián)系的主要手段，該方法自Cundall等[1]提出后，在模型建立、算法實現(xiàn)、接觸本構等方面得到快速發(fā)展，可以實時監(jiān)測加載過程中的顆粒位置及顆粒間接觸情況. 利用該方法，國內外學者開展了大量數(shù)值模擬. 在簡單加載條件下，Garcia等[2]進行了大量的直剪試驗模擬，探究顆粒形狀對抗剪強度的影響規(guī)律. Nouguier-Lehon等[3]對不規(guī)則多邊形顆粒進行了一系列雙軸壓縮數(shù)值試驗，結果表明隨顆粒長徑比的增加，試樣達到臨界狀態(tài)所需的應變逐漸增加. 吳越等[4]采用幾種不同的制樣方式，研究了初始各向異性程度對剪脹性及臨界狀態(tài)的影響規(guī)律. 蔣明鏡等[5]分析了顆粒排列方向對試樣變形和強度的影響，發(fā)現(xiàn)隨著沉積角的增加，試樣的剪脹程度逐漸減弱，殘余內摩擦角逐漸減小. 在復雜加載條件下，劉洋等[6-7]使用顆粒流軟件PFC2D開展循環(huán)加載試驗，分析了松砂和密砂的顆粒配位數(shù)、粒間接觸力等細觀物理量在漸進破壞中的演化規(guī)律. 周健等[8]分析了細觀組構變化對宏觀力學響應的影響，探究了砂土液化的細觀機制. Sazzad等[9]發(fā)現(xiàn)，顆粒間摩擦角的取值越大，循環(huán)加載時應力- 應變滯回圈變得越來越窄. Phusing等[10]模擬了4種循環(huán)三軸壓縮試驗，采用配位數(shù)和滑動接觸率解釋了循環(huán)加載中的宏觀力學行為. 以上離散元模擬表明，顆粒形狀、排列等細觀結構特征與試樣變形、強度等宏觀力學性質之間有著密切聯(lián)系，但上述研究僅僅對這一聯(lián)系進行了定性分析，尚未給出宏細觀之間的定量關系.

受限于顆粒材料力學性質的復雜性，宏細觀定量關系的研究相對較少，一些學者主要是通過改變離散元模擬中顆粒細觀參數(shù)的取值，獲得細觀參數(shù)與宏觀力學參數(shù)之間的擬合關系式. 例如，Yang等[11]提出了平行黏結模型中顆粒剛度、顆粒尺寸等細觀參數(shù)和彈性模量、泊松比等宏觀力學指標之間的經驗公式. 徐小敏等[12]對離散元模擬和室內三軸試驗結果做回歸分析，建立了宏細觀彈性參數(shù)之間的經驗公式. 孔亮等[13]采用形狀系數(shù)量化顆粒的圓度和凹凸度，通過雙軸試驗離散元模擬提出了峰值強度和形狀系數(shù)的擬合公式. 周博等[14]研究了顆粒黏結強度、摩擦因數(shù)、剛度比等細觀參數(shù)和材料抗剪強度指標之間的相關性. 除了建立宏細觀參數(shù)之間的關系式以外，Zhao等[15]根據(jù)應力張量和組構張量的乘積定義了一個標量，發(fā)現(xiàn)在臨界狀態(tài)時，該標量的取值與中主應力系數(shù)無關，可以表示成平均有效正應力的單變量函數(shù). 從以上文獻綜述中可以發(fā)現(xiàn)，關于顆粒材料宏細觀定量關系的研究，還處在探索階段. 一些學者建立的細觀接觸模型參數(shù)與宏觀本構模型參數(shù)之間的關系式，與所選取的模型有關. 若能得到組構張量與應力張量之間的關系，則可以在現(xiàn)有的很多宏觀模型的基礎上，發(fā)展多尺度本構模型，既能描述顆粒材料的變形機理，又能在實際工程中得到應用.

為了建立宏細觀之間的定量聯(lián)系，本文開展了一系列應力比保持恒定的雙軸壓縮離散元模擬，分析顆粒之間的接觸力(細觀量)與試樣受到的外力(宏觀量)之間的規(guī)律，得到了材料各向異性程度與應力比的函數(shù)關系式，進而構建了組構張量與應力張量的隱式定量聯(lián)系.

1 等應力比離散元模擬

土體的變形、強度等宏觀力學特性與土顆粒的大小、形狀、排列方式等細觀結構特征之間有著密切聯(lián)系. 在加載過程中，隨著應力狀態(tài)的改變，顆粒的排列方式不斷演化，因而直接建立復雜加載條件下的宏細觀定量關系是十分困難的. 為簡化起見，本文采用離散單元法，開展等應力比雙軸壓縮試驗，此時宏觀物理量(即應力比)是明確且恒定的，只需要觀察細觀物理量(即顆粒組構)是如何演化的，從而減少了宏細觀關系式的變量個數(shù). 之所以采用等應力比這種加載方式，也是考慮到顆粒材料具有摩擦特性，決定其屈服和破壞的是應力比，而并非某一應力的絕對值，通過對比不同應力比下的試驗結果，可以更直觀地研究顆粒材料的變形和強度規(guī)律.

1.1 加載方案

借助二維顆粒流程序PFC2D開展一系列雙軸壓縮模擬. 采用ball distribute功能，生成一個初始尺寸為70 mm×70 mm、孔隙率為0.15的試樣. 其中顆粒均為圓形，直徑為0.1～0.6 mm，密度為2 500 kg/m3，總數(shù)超過9 000個. 然后，為了保證顆粒的初始排列為各向同性，給試樣施加100 kPa的等向壓力. 將按上述方法制得的試樣保存，記錄每個顆粒的位置以及顆粒間接觸力的大小和方向，后續(xù)加載均在此顆粒排列的基礎上進行.

在等應力比加載過程中，豎向應力σy從100 kPa開始，每個增量步增加5 kPa；同時，等比例地施加水平應力σx，使得σy與σx的比值r始終不變.共進行10組數(shù)值試驗，r的取值分別為1.1，1.2，…，2.0.以r=2為例，加載開始時，σy=100 kPa，σx=50 kPa，第一個增量步完成后，σy增加至105 kPa，而σx變?yōu)?2.5 kPa.等應力比加載相當于對試樣進行各向異性壓縮，試樣內部的顆粒排列將隨之發(fā)生改變，以適應外荷載的變化.在加載過程中，實時監(jiān)測顆粒之間接觸的變化，直至接觸法向的空間分布達到基本穩(wěn)定.

1.2 模擬參數(shù)

在PFC2D中，顆粒本身被視為不可變形的剛體，但允許相鄰兩顆粒存在一定程度的重疊，并根據(jù)重疊量的大小，通過接觸模型計算顆粒間的接觸力. 本文采用最簡單的線性接觸模型，來探究宏細觀之間的定量聯(lián)系. 模型參數(shù)的取值參考蔣明鏡等[16]的推薦值，如表1所示，該組參數(shù)已被證明能夠使得離散元模擬結果與室內試驗結果基本一致.

表1 數(shù)值模擬參數(shù)

2 模擬結果及定性分析

首先，通過觀察雙軸壓縮過程中顆粒間接觸力的傳遞和接觸法向的空間分布，定性地分析二者與應力比之間的規(guī)律，為下文建立宏細觀之間的定量關系提供依據(jù).

2.1 接觸力的傳遞

由于采用圓形顆粒，將相鄰兩顆粒的圓心相連，即可得到法向接觸力的作用軌跡，如圖1所示. 該圖能夠顯示外荷載是如何在試樣內部傳遞的，圖中連線的粗細反映接觸力的大小. 圖1分別對應初始狀態(tài)，以及r=1.5、2.0的等應力比加載的穩(wěn)定時刻. 可知，在等應力比加載前，各方向接觸力的傳遞較為均衡，連線的粗細也相差不大，表明接觸力的大小在不同方向上基本相同，試樣為各向同性. 經過等應力比加載后，水平方向的接觸力減弱、甚至部分消失，豎直方向成為接觸力傳遞的主導方向. 另一方面，通過連線粗細的變化可以看出，隨著應力比的增大，豎直方向的接觸力較水平方向越來越大. 特別是當r=2.0時，可以很明顯地看到幾條自上而下的接觸力鏈(見圖1(c)的放大區(qū)域)，這是試樣為承受外荷載而形成的各向異性細觀結構.

圖1 不同應力比下的接觸力鏈Fig.1 Contact force chains under different stress ratios

由以上分析可知，顆粒間接觸力(細觀物理量)的方向、大小受外荷載(宏觀物理量)的方向、大小的顯著影響. 當某一方向上的應力較大時，該方向成為接觸力傳遞的主導方向，其接觸力的大小明顯高于其他方向. 因此，宏觀應力與細觀接觸力之間存在密切聯(lián)系.

2.2 接觸法向的空間分布

為了更加直觀地分析顆粒排列的各向異性，統(tǒng)計了顆粒之間所有接觸法向的空間分布情況，繪出玫瑰圖，如圖2所示. 采用極坐標系，每個扇形對應的轉角，等于接觸法向與水平方向的夾角；扇形的半徑，表示該方向上的接觸數(shù)占總接觸數(shù)的百分比. 根據(jù)圖2(a)可知，初始狀態(tài)時，各方向接觸法向的數(shù)目基本相同，表明顆粒排列接近各向同性. 而圖2(b)～(d)顯示，在等應力比加載完成后，接觸法向的空間分布發(fā)生顯著變化，豎直方向(即轉角等于90°)的接觸法向明顯增多，且當r=2.0時，將近12%的接觸法向沿著豎直方向，而沿水平方向(即轉角等于0°或180°)的接觸法向數(shù)僅占總數(shù)的7%～8%.

圖2 不同應力比下接觸法向的空間分布Fig.2 Spatial distribution of contact normal under different stress ratios

3 組構張量及其演化

在定性分析的基礎上引入組構張量的概念，定量地描述接觸法向的空間分布及其演化規(guī)律，進而尋找宏細觀之間的定量關系.

3.1 組構張量的定義

組構張量是描述顆粒材料各向異性的物理量，可以根據(jù)顆粒長軸方向[17]、孔隙長軸方向[18]、接觸法向[19]等細觀幾何要素進行定義. 由于本文的離散元模擬采用圓形顆粒，不存在顆粒長軸，但接觸法向能夠非常方便地確定；而且有研究表明[20]，相比于顆粒長軸方向和孔隙長軸方向，接觸法向對應力的敏感程度更高，更容易隨著外荷載的改變而實時調整其空間分布，因此本文利用接觸法向來定義組構張量. 如圖3所示，接觸法向垂直于相鄰兩顆粒的公切線，與顆粒圓心連線的方向相同，沿該方向的單位向量為n=(nx,ny)(x軸沿水平方向，y軸沿豎直方向)，參考Oda[21]的定義式，計算組構張量

圖3 組構張量的定義Fig.3 Definition of fabric tensor

(1)

3.2 組構演化規(guī)律

直接分析一個二階張量的變化規(guī)律較為困難，因此，首先引入Fij的第二偏不變量

(2)

式中：δij為Kronecker符號(當i=j時，δij=1；當i≠j時，δij=0).因此，當且僅當Fij=δij/2，即試樣為各向同性時，J才等于0；否則J≠0，代表試樣為各向異性，且J值越大，各向異性程度越高.

在等應力比加載過程中，通過對接觸法向的實時監(jiān)測，得到了J隨豎直應力σy的演化規(guī)律，如圖4所示.以應力比r=1.5的曲線為例進行分析.在加載初期，J接近0，說明此時試樣基本處于各向同性狀態(tài).隨著σy的提高，J逐漸增大，試樣的各向異性程度變大，這是因為顆粒之間的接觸法向在不斷調整，趨向于豎直方向.此后，J的增長速率逐漸變緩，最終收斂到一個穩(wěn)定值，外荷載的繼續(xù)提高不會改變試樣內部的細觀結構，顆粒排列已達到穩(wěn)定狀態(tài).對比不同應力比加載下的J演化曲線，可以看出：r越大，J增長速度越快，表明接觸法向的調整變得更加迅速；另一方面，最終J的穩(wěn)定值也越來越大，且對應的σy越來越高，這說明試樣需要施加更大的外荷載才能達到穩(wěn)定，所形成的細觀結構具有更強的各向異性程度.

圖4 不同應力比下第二偏不變量的演化規(guī)律Fig.4 Evolution of the second partial invariant under different stress ratios

3.3 穩(wěn)定時的組構值

由于J值在σy較大時仍不斷波動，僅憑肉眼觀察確定J的穩(wěn)定值具有很強的主觀性，因此引入統(tǒng)計學中的變異系數(shù)δ作為判斷J是否達到穩(wěn)定的標準，其表達式為

(3)

圖5的橫坐標為應力比r，縱坐標為相應的Js.由圖可知，當r∈[1.1～2.0]時，隨著r的增大，Js逐漸增長(說明顆粒排列達到穩(wěn)定時的各向異性程度更高)，且增長速度越來越快.這一規(guī)律可用二次函數(shù)描述，即

圖5 第二偏不變量的穩(wěn)定值與應力比的關系Fig.5 Relation between the stable value of the second partial invariant and the stress ratio

Js=a(r-1)2

(4)

式中a為反映Js增長速度的系數(shù).與數(shù)值模擬結果的對比表明，式(4)可以很好地擬合Js隨r的變化規(guī)律，相關系數(shù)R2達到0.991.而且，當r=1時，根據(jù)式(4)可得Js=0，這符合等向壓縮條件下顆粒排列最終演化成各向同性的試驗規(guī)律.

3.4 適用性分析

式(4)是在特定情況下得到的擬合結果，為了檢驗其對其他情況的適用性，調整了顆粒的初始排列方式和摩擦因數(shù)的取值，采用相同的等應力比加載和數(shù)據(jù)處理方法，分析了上述2種因素對模擬結果的影響.

通過對隨機生成的試樣施加100 kPa的等向預壓力，獲得初始各向同性試樣；為得到各向異性試樣，將某方向的預壓力提高，進行非等向壓縮，使得顆粒的初始排列產生一定的方向性，再進行r=1.8下的等應力比加載.得到的J值隨σy的演化規(guī)律如圖6所示.對于試樣I，水平方向的預壓力大于豎直方向，因此初始時顆粒間的接觸力以水平方向為主(Fxx>Fyy).在等應力比加載階段，隨著σy的提高，原本占主導的水平接觸力逐漸減弱，試樣變?yōu)楦飨蛲?Fxx=Fyy)，因此J先減小到0.然后，隨著豎向接觸力的生成和不斷發(fā)展，豎直方向成為接觸力傳遞的主導方向(Fxx

圖6 初始顆粒排列對組構演化的影響Fig.6 Effect of the initial particle arrangement on fabric evolution

圖7為顆粒間摩擦因數(shù)μ分別取0.3、0.4、0.5時，Js隨r的變化規(guī)律.可知，在其他條件相同的情況下，μ的取值對Js有一定影響，因為顆粒間的摩擦力(即切向接觸力)也能夠分擔一部分外荷載，從而改變最終穩(wěn)定時的顆粒排列方式.然而，μ只能影響擬合系數(shù)a的大小，而在本文的模擬條件下，Js與r的二次關系仍然成立，且相關系數(shù)R2均大于0.980.

圖7 顆粒間摩擦系數(shù)對第二偏不變量的穩(wěn)定值的影響Fig.7 Effect of the particle friction coefficient on the stable value of the second partial invariant

4 組構張量與應力張量的定量聯(lián)系

得到了組構張量的第二偏不變量J與應力比r的公式，初步建立了宏細觀之間的定量聯(lián)系，但J和r都是標量，只能反映組構和應力的各向異性程度的大小，而無法表示二者的方向.若要全面考慮大小和方向，需要建立組構張量Fij與應力張量σij的關系式，即Fij=f(σij)或σij=f(Fij).然而，直接得到2個張量之間的顯式方程較為困難.本文借鑒Zhao等[15]的方法，先將Fij與σij以不同形式組合，得到一個新的二階張量；若該組合張量與加載方式無關，即f(Fij,σij)=C(C為常數(shù))，說明Fij與σij始終滿足一個隱式函數(shù)關系.該隱式函數(shù)也能夠描述宏細觀定量關系.

4.1Fij與σij的組合形式

首先，細觀物理量，即組構張量Fij的定義式見式(1)；而對于宏觀物理量，根據(jù)雙軸壓縮模擬中豎直和水平加載板上的平均應力，得到應力張量

(5)

為了消除應力大小的影響，將σij各元素都除以平均正應力p(p=(σy+σx)/2)，得到量綱一的應力比張量

(6)

在等應力比加載過程中，ηij不發(fā)生變化.

然后，將Fij與ηij以不同形式做點積，構造了如下幾種組合張量：

(7)

(8)

(9)

4.2 組合張量隨應力比的變化規(guī)律

仿照式(2)，求出顆粒排列達到穩(wěn)定狀態(tài)時Xij對應的第二偏不變量，仍記為Js，進而繪出Js與r的關系，如圖8所示.因為對于不同形式的Xij，Js的變化范圍有所不同，所以利用r=1.1時的Js值進行了歸一化，以便比較Js隨r的相對變化情況.

圖8 各種組合張量的第二偏不變量隨應力比的變化規(guī)律Fig.8 Variation of the second partial invariants of different joint tensors with stress ratio

需要說明的是，根據(jù)Zhao等[15]的工作，在臨界狀態(tài)，即應力比等于臨界狀態(tài)應力比時，F(xiàn)ij與σij的某種聯(lián)合不變量，是一個與中主應力系數(shù)無關的常量.本文的工作可視為對上述規(guī)律的發(fā)展，即在應力比達到臨界狀態(tài)應力比之前，試樣壓縮穩(wěn)定(顆粒排列不再隨應力的增加而改變)時Fij與σij的某種組合，也是一個常量.

為了得到顆粒材料宏細觀之間的定量關系，本文在離散元模擬中進行了大幅簡化：顆粒形狀為圓形，所有顆粒之間的摩擦因數(shù)取均一值，并且采用二維模擬，加載條件為等應力比壓縮，直至顆粒排列達到穩(wěn)定狀態(tài).而實際中，土顆粒的細觀特性十分復雜，土體受到的應力狀態(tài)多種多樣，因此本文所建立的宏細觀定量關系式存在一定的局限性，相關問題還應進行深入研究，以更加貼近真實情況.

5 結論

1) 在等應力比加載過程中，大主應力方向逐漸成為顆粒間接觸力傳遞的主導方向，且應力比越大，這一趨勢越明顯，顆粒排列的各向異性程度越高.

2) 隨著外荷載的增加，組構張量的第二偏不變量逐漸趨于一個穩(wěn)定值，該穩(wěn)定值與應力比呈二次函數(shù)關系，且這種關系不隨顆粒的初始排列方式、摩擦因數(shù)取值的改變而變化.

3) 將顆粒排列達到穩(wěn)定時的組構張量與應力張量按某種形式進行組合，該組合張量與應力比無關，從而建立了組構張量和應力張量的隱式關系.