孔隙空間分布對(duì)粗粒土強(qiáng)度變形特性影響研究

胡有方,袁俊平，盧 毅

(1. 河海大學(xué) 巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2. 河海大學(xué) 江蘇省巖土工程技術(shù)工程研究中心，江蘇 南京 210098)

粗粒土是一種典型的顆粒材料，由大量粒徑不等、形狀各異、排列隨機(jī)的土石顆粒組成，微觀上孔隙結(jié)構(gòu)的多樣性決定了其宏觀物理力學(xué)性質(zhì)的復(fù)雜性[1]。因此可以說(shuō)，粗粒土的變形過(guò)程實(shí)質(zhì)上就是其孔隙結(jié)構(gòu)變化的過(guò)程，如果能獲得粗粒土的孔隙結(jié)構(gòu)特征和變化規(guī)律，也就能掌握其變形規(guī)律[2]。近年來(lái)，許多國(guó)內(nèi)外學(xué)者對(duì)粗粒土的孔隙結(jié)構(gòu)特征進(jìn)行了研究，主要方法有實(shí)驗(yàn)法和數(shù)值重構(gòu)法[3-4]。實(shí)驗(yàn)方法最早起源于1982年，Petrovic等[5]首先在土壤的密度和含水率研究中引入CT掃描技術(shù)，之后，Warner[6]則在對(duì)CT掃描圖像的分析中發(fā)現(xiàn)，孔隙的位置、大小和數(shù)量信息都能由圖像精確地揭示。在國(guó)內(nèi)，學(xué)者呂菲[7]等在掃描圖像的基礎(chǔ)上建立了孔隙網(wǎng)絡(luò)模型，并成功預(yù)測(cè)了飽和土壤的水力學(xué)性質(zhì)。相方園[8]使用切片法進(jìn)行試驗(yàn)，將孔隙拓?fù)浣Y(jié)構(gòu)與形態(tài)學(xué)理論相結(jié)合，定量地獲得了孔隙結(jié)構(gòu)中的特征參數(shù)。與實(shí)驗(yàn)方法不同，數(shù)值重構(gòu)法的思路是基于二維孔隙信息和不同的數(shù)學(xué)模擬算法建立三維孔隙結(jié)構(gòu)模型，目前常用的方法有高斯隨機(jī)場(chǎng)法[9]、過(guò)程法[10]和模擬退火法[11]。比較成功的案例有Pilotti[12]基于過(guò)程法，用一種主要適用于球體、橢球體、圓柱體和平行六面體的算法建立了任意形狀顆粒的隨機(jī)堆積體。曾建邦[13]等利用改進(jìn)的模擬退火法重建了含水沉積物的三維孔隙結(jié)構(gòu)，并通過(guò)孔隙結(jié)構(gòu)的特征化分析揭示了分布模式對(duì)沉積物特性的影響。

綜上所述，目前的研究主要針對(duì)的是粗粒土孔隙結(jié)構(gòu)特征的描述方法，因此，如何將粗粒土的孔隙結(jié)構(gòu)特征與強(qiáng)度變形特性定量地聯(lián)系起來(lái)，仍是有待探討的問(wèn)題。本文在前人的研究基礎(chǔ)上，從孔隙空間分布的不均勻性出發(fā)進(jìn)行研究，定義了可以定量刻畫粗粒土孔隙空間變異性的孔隙空間分布系數(shù)FSD；基于PFC3D離散元軟件，分析了FSD與粗粒土強(qiáng)度變形參數(shù)之間的關(guān)系；同時(shí)，基于Abaqus有限元軟件，分析了FSD在有限元計(jì)算中的適用性。

1 土體孔隙空間分布不均勻性的定量刻畫

首先，粗粒土試樣中各巖土顆粒和孔隙的位置可以通過(guò)掃描的方法得到[14]。接下來(lái)，使用最大球算法[15]對(duì)孔隙進(jìn)行模擬，即對(duì)于每一部分孔隙，插入若干個(gè)以巖土顆粒為邊界的最大內(nèi)切球，使這些球充滿孔隙空間。以各球體的質(zhì)心作為孔隙的質(zhì)心，以各球體的體積作為孔隙的體積，從而達(dá)到以大量規(guī)則的球體模擬不規(guī)則的孔隙體的目的，最大球算法的示意圖如圖1所示。

圖1 最大球算法示意圖Fig.1 Schematic diagram of maximum sphere algorithm

下一步，將土體劃分為N個(gè)相同的區(qū)域?qū)紫肚蝮w進(jìn)行統(tǒng)計(jì)，由于粗粒土的不均勻性，每個(gè)區(qū)域內(nèi)的孔隙數(shù)量和孔隙體積分布都不相同，因此引入孔隙數(shù)量變異系數(shù)CVn和孔隙體積變異系數(shù)CVv，兩者的表達(dá)式分別為

(1)

(2)

CVn和CVv可以較好地將試樣中孔隙數(shù)量的分布情況和孔隙體的分布情況反映出來(lái)，其值越小，則代表均勻性越好。理想情況下，孔隙分布完全均勻時(shí)，兩者值為0。但要注意的是，這兩個(gè)參數(shù)本身無(wú)法反映出孔隙質(zhì)心點(diǎn)在空間中的分布情況，因?yàn)榇嬖谄郜F(xiàn)象，如圖2所示，a試樣與b試樣的CVn和CVv相同，但顯然a試樣的偏聚程度小于b試樣。因此，還需要引入偏聚系數(shù)β來(lái)反映這種情況。

圖2 孔隙空間分布狀態(tài)示意圖Fig.2 Diagram of pore space distribution

(3)

式中：n表示樣本區(qū)域內(nèi)孔隙總數(shù)量。

接著，用x0表示試樣中孔隙完全規(guī)則地分布在試樣空間中時(shí)孔隙之間的距離，然而對(duì)于任意的n值，x0的值都難以計(jì)算，因此本文用下式代替：

(4)

式中：v表示樣本區(qū)域內(nèi)孔隙總體積。如此計(jì)算時(shí)，x0表示當(dāng)n個(gè)相同尺寸的球體孔隙總體積與樣本區(qū)域內(nèi)孔隙總體積相等時(shí)這些孔隙的直徑。

于是，偏聚系數(shù)β可以用下式計(jì)算：

(5)

在得到了孔隙數(shù)量變異系數(shù)CVn，孔隙體積變異系數(shù)CVv以及偏聚系數(shù)β之后，本文定義孔隙空間分布系數(shù)FSD作為綜合性指標(biāo)，從孔隙質(zhì)心點(diǎn)空間分布的均勻性、孔隙體積空間分布的均勻性以及孔隙質(zhì)心在空間上的偏聚程度三個(gè)方面綜合反映孔隙空間的分布情況，其表達(dá)式定義為

FSD=β·CVn·CVv

(6)

2 孔隙空間變異系數(shù)與土體強(qiáng)度變形參數(shù)關(guān)系研究

在定義了孔隙空間分布系數(shù)FSD之后，為了研究其對(duì)土體宏觀強(qiáng)度變形特性的具體影響，就需要定量分析FSD與土體各項(xiàng)強(qiáng)度參數(shù)之間的關(guān)系。因此本節(jié)中基于PFC3D軟件，建立了5組三軸壓縮試驗(yàn)試樣進(jìn)行數(shù)值模擬實(shí)驗(yàn)。試樣編號(hào)為PK1—PK5，高度H=200 mm，直徑D=101 mm。這些試樣的總孔隙率相同，但各試樣的內(nèi)部分層情況不同，詳細(xì)參數(shù)和示意圖如圖3所示。

試樣模型按照以下步驟建立：(1)將試樣空間按照試驗(yàn)方案均分為指定層數(shù)，從最底層開始生成顆粒；(2)使用半徑擴(kuò)大法生成試驗(yàn)方案中指定孔隙率的土顆粒，在每層頂部多預(yù)留5 cm的高度；(3)賦予土顆粒重力加速度，使土顆粒在重力作用下自由下落形成堆積體；(4)賦予該層頂部墻體一定速度使其下落直至獲得該層指定尺寸的試樣，運(yùn)行一定時(shí)步使試樣達(dá)到穩(wěn)定狀態(tài)，該層試樣生成完畢；(5)重復(fù)(2)—(4)步驟直至所有層的試樣均生成完畢，刪除中間墻體，運(yùn)行一定時(shí)步使試樣達(dá)到穩(wěn)定狀態(tài)，試樣生成完畢。模型生成過(guò)程如圖4所示。

完成試樣的生成后，為了防止偶然誤差給數(shù)值試驗(yàn)結(jié)果帶來(lái)影響，對(duì)各試樣進(jìn)行多次對(duì)照試驗(yàn)，結(jié)果取平均值。最終統(tǒng)計(jì)各試樣的孔隙空間分布系數(shù)FSD如表1所示。

注：數(shù)據(jù)為各層孔隙率。圖3 不同孔隙空間分布影響效應(yīng)試驗(yàn)方案示意圖Fig.3 Schematic diagram of test scheme for effect of different pore space distribution

圖4 模型生成過(guò)程示意圖Fig.4 Diagram of model generation process

從表1可以看出，在其他顆粒細(xì)觀參數(shù)相同的條件下，各試樣的孔隙空間分布系數(shù)因分層的不同而發(fā)生了明顯的變化。試樣PK5相對(duì)于試樣PK1，空間分布系數(shù)FSD增大了46.8%，說(shuō)明該方法的合理性。

為了得到上述各試樣的應(yīng)力和應(yīng)變曲線，對(duì)表中的各組試樣進(jìn)行圍壓800 kPa的三軸固結(jié)排水試驗(yàn)數(shù)值模擬，結(jié)果如圖5和圖6所示。

從圖5和圖6中可以看出，在其他顆粒細(xì)觀參數(shù)相同的條件下，土體的彈性模量、峰值強(qiáng)度和泊松比都隨著孔隙空間分布系數(shù)FSD的增大而減小。試樣PK5相對(duì)于試樣PK1，彈性模量減小40.8%，峰值強(qiáng)度降低12.5%，泊松比減小39.5%。

表1 各試樣孔隙空間分布系數(shù)

圖5 偏應(yīng)力-軸向應(yīng)變曲線Fig.5 Relationship between deviatoric stress and axial strain

圖6 側(cè)向應(yīng)變-軸向應(yīng)變曲線Fig.6 Relationship between lateral strain and axial strain

由于加載剛進(jìn)行時(shí)，試樣的偏應(yīng)力-軸向應(yīng)變曲線近似為一條直線，因此將試樣軸向應(yīng)變達(dá)到1%時(shí)的割線模量作為試樣的彈性模量。相關(guān)性分析結(jié)果表明，彈性模量E與孔隙空間分布系數(shù)FSD之間呈指數(shù)函數(shù)關(guān)系，如式(7)。

E=A1+B1er1FSD

(7)

用指數(shù)函數(shù)對(duì)試樣彈性模量E和孔隙空間分布系數(shù)FSD進(jìn)行擬合，結(jié)果如圖7所示，其中A1的值為94.062 2，B1的值為-0.609 4，r1的值為-2 071.968 1。

將軸向應(yīng)變1%時(shí)的切線泊松比作為試樣的泊松比，相關(guān)性分析結(jié)果表明，泊松比ν與孔隙空間分布系數(shù)FSD呈指數(shù)函數(shù)關(guān)系，如式(8)。

v=A2+B2er2FSD

(8)

圖7 彈性模量隨空間分布系數(shù)變化曲線Fig.7 Relationship between elastic modulus and pore space distribution coefficient

圖8 泊松比隨空間分布系數(shù)變化曲線Fig.8 Relationship between Poisson′s ratio and pore space distribution coefficient

用指數(shù)函數(shù)對(duì)二者進(jìn)行擬合，結(jié)果如圖8所示，其中A2的值為0.034 3，B2的值為-21.302 8，r2的值為-3 967.103 5。

為了得到各試樣的內(nèi)摩擦角，在圍壓200、800和2 000 kPa的條件下，對(duì)各試樣進(jìn)行三軸試驗(yàn)數(shù)值模擬，得到了各試樣在不同圍壓下的峰值強(qiáng)度。根據(jù)莫爾-庫(kù)倫強(qiáng)度理論，各試樣的內(nèi)摩擦角如表2所示。

表2 各試樣內(nèi)摩擦角

為了研究試樣內(nèi)摩擦角與孔隙空間分布系數(shù)之間的定量關(guān)系，對(duì)內(nèi)摩擦角與孔隙空間分布系數(shù)進(jìn)行相關(guān)性分析，發(fā)現(xiàn)兩者之間呈指數(shù)函數(shù)關(guān)系，如式(9)。

φ=A3+B3er3FSD

(9)

用指數(shù)函數(shù)對(duì)試樣內(nèi)摩擦角和空間分布系數(shù)進(jìn)行擬合，結(jié)果如圖9所示，其中A3的值為38.469 2，B3的值為-0.156 7，r3的值為1 705.456 1。

3 孔隙空間分布系數(shù)的準(zhǔn)確性驗(yàn)證與分析

本節(jié)中使用Abaqus有限元軟件對(duì)孔隙空間分布系數(shù)FSD計(jì)算方法的準(zhǔn)確性進(jìn)行驗(yàn)證。相比于離散元方法，有限元方法更適合大規(guī)模的工程計(jì)算，這是由于有限元將粗粒土復(fù)雜的幾何結(jié)構(gòu)簡(jiǎn)化為了具有簡(jiǎn)單形狀的單元，單元內(nèi)的材料性質(zhì)和控制方程通過(guò)單元節(jié)點(diǎn)的未知量來(lái)進(jìn)行表達(dá)，從而使得計(jì)算的效率大大提高，然而這也使得有限元計(jì)算時(shí)忽略或低估了孔隙空間分布系數(shù)FSD的影響。那么孔隙空間分布系數(shù)FSD是否可以被忽略，是否可以使用Abaqus有限元軟件進(jìn)行孔隙空間變異性的模擬，這是本節(jié)中要討論的問(wèn)題。

為了解決此問(wèn)題，建立編號(hào)分別為A1—A5的5種三軸固結(jié)排水實(shí)驗(yàn)試樣，試樣的高度、直徑均與PK1—PK5相同。同時(shí)，通過(guò)上節(jié)中的結(jié)果推算出試樣的彈性模量、泊松比、內(nèi)摩擦角，使得兩組試樣的等效強(qiáng)度參數(shù)全部相同，具體數(shù)值如表3所示。最后，設(shè)置試樣的孔隙率均為35%，使其為均勻試樣。

表3 A1—A5試樣強(qiáng)度參數(shù)

對(duì)試樣施加800 kPa的圍壓時(shí)，各試樣的偏應(yīng)力-軸向應(yīng)變曲線如圖10所示。

圖10 試樣A1—A5偏應(yīng)力-軸向應(yīng)變曲線Fig.10 Relationship between deviatoric stress and axial strain of sample A1—A5

從圖10中可以看出，由于試樣A1—A5忽略了孔隙空間分布系數(shù)FSD，消除了空間分布差異對(duì)強(qiáng)度的不利影響，使得雖然試樣的各項(xiàng)強(qiáng)度參數(shù)均與離散元計(jì)算時(shí)相同，但試樣破壞時(shí)的偏應(yīng)力均偏大，并且隨著離散元試樣中FSD的增加，這種差異更加明顯。試樣A1相比試樣PK1，破壞時(shí)的偏應(yīng)力增大了1.6%，試樣A5相比試樣PK5，破壞時(shí)的偏應(yīng)力增大了12.0%。

為了降低有限元單元的均勻性，增加孔隙空間分布系數(shù)的影響，采用上節(jié)圖3中的分層方法建立5種三軸固結(jié)排水實(shí)驗(yàn)試樣B1—B5，在前者的基礎(chǔ)上額外考慮不同土層孔隙率的差異。同樣對(duì)試樣施加800 kPa的圍壓時(shí)，各試樣的偏應(yīng)力-軸向應(yīng)變曲線如圖11所示。以試樣5為例，將三次試驗(yàn)的結(jié)果進(jìn)行比較，如圖12所示。

圖11 試樣B1—B5偏應(yīng)力-軸向應(yīng)變曲線Fig.11 Relationship between deviatoric stress and axial strain of sample B1—B5

圖12 試樣A5,B5,PK5偏應(yīng)力-軸向應(yīng)變曲線Fig.12 Relationship between deviatoric stress and axial strain of sample A5,B5,PK5

從圖11中可以看出，對(duì)于考慮分層的有限元試樣，其峰值強(qiáng)度相比均勻的有限元試樣有所降低，但仍然比離散元試樣大。這是因?yàn)?，雖然分層增加了層與層之間的孔隙空間分布差異，但單獨(dú)每個(gè)層內(nèi)的孔隙仍然是均勻的，其FSD雖然大于0，但仍小于離散元試樣。

從圖12中可以看出，對(duì)試樣5來(lái)說(shuō)，雖然三次試驗(yàn)的等效強(qiáng)度完全相同，但三次試驗(yàn)的峰值偏應(yīng)力分別為2.74，3.07，2.99 MPa。再次驗(yàn)證了試樣的峰值強(qiáng)度隨著孔隙空間分布系數(shù)FSD的增大而減小的規(guī)律，并進(jìn)一步論證了孔隙空間分布系數(shù)FSD計(jì)算方法的合理性。同時(shí)也說(shuō)明，在試樣離散型較大時(shí)，不應(yīng)當(dāng)忽略孔隙空間分布變異性的影響。

4 結(jié)論

1)孔隙空間分布系數(shù)FSD可以較好地模擬粗粒土中孔隙空間分布的不均勻性。這種不均勻性包括孔隙質(zhì)心點(diǎn)空間分布的不均勻性、孔隙體積空間分布的不均勻性以及孔隙質(zhì)心在空間上不同的偏聚程度，因此FSD是一項(xiàng)比較系統(tǒng)的綜合性指標(biāo)。

2)孔隙空間分布系數(shù)FSD在離散元和有限元分析中均能得到較好的運(yùn)用。在其他顆粒細(xì)觀參數(shù)相同的條件下，土體的彈性模量、峰值強(qiáng)度和泊松比都隨著孔隙空間分布系數(shù)FSD的增大呈指數(shù)函數(shù)形式減小。