描述砂性土宏細(xì)觀力學(xué)關(guān)系的數(shù)學(xué)模型

閆洪超， 魯杰， 饒振興， 翟洪濤， 關(guān)鵬

(1. 河南省地質(zhì)礦產(chǎn)勘查開發(fā)局第四地質(zhì)勘查院， 鄭州 450000； 2. 中國(guó)地質(zhì)大學(xué)工程學(xué)院， 武漢 430074)

由地質(zhì)作用或人工碎裂形成的土既可作為建筑材料，如磚、瓦、混凝土骨料等，又可作為地基或路基承受荷載[1-2]。其力學(xué)特性是土力學(xué)最關(guān)注的問題之一。土體總是被當(dāng)成連續(xù)介質(zhì)，這是因?yàn)?，?dāng)土被視為連續(xù)介質(zhì)時(shí)，可采用連續(xù)力學(xué)理論來(lái)研究土體，進(jìn)而構(gòu)建基于連續(xù)介質(zhì)理論的應(yīng)力-應(yīng)變本構(gòu)關(guān)系，用于服務(wù)工程建設(shè)。實(shí)踐證明，采用連續(xù)介質(zhì)力學(xué)理論研究土體是可行的，解決了很多實(shí)際工程問題。然而土體畢竟是由土粒構(gòu)成的，把它視為連續(xù)介質(zhì)理論上是行不通的，土體力學(xué)行為的研究最終應(yīng)回歸土粒或土粒群[3]。自顆粒流離散元方法問世以來(lái)，砂性土細(xì)觀顆粒層面的力學(xué)行為得到了深入的研究，為構(gòu)建土體宏觀力學(xué)響應(yīng)與土粒相關(guān)的細(xì)觀力學(xué)行為之間的聯(lián)系提供了途徑，突破性的研究成果不斷出現(xiàn)，Christoffersen等[4]基于準(zhǔn)靜態(tài)平衡理論，提出了采用粒間接觸力表述砂性土等顆粒材料宏觀應(yīng)力張量的方法。Oda[5]提出采用二階組構(gòu)張量描述顆粒材料顆粒定向、接觸力大小和方向的各向異性。砂性土在荷載作用下宏觀應(yīng)力張量發(fā)生變化時(shí)，必然會(huì)有粒間接觸力各向異性的變化。Rothenburg等[6]提出采用傅里葉級(jí)數(shù)描述粒間接觸力的各向異性，并建立起了細(xì)觀接觸力各向異性指標(biāo)與宏觀應(yīng)力比的關(guān)系，同樣的，砂性土在荷載作用下宏觀應(yīng)變張量發(fā)生變化時(shí)，必然會(huì)對(duì)應(yīng)著與體積相關(guān)的某一細(xì)觀指標(biāo)各向異性的變化。Kruyt等[7]將相互接觸的顆粒形心連接成多邊形，且定義了多邊形的方向，試樣在變形過程中多邊形的方向表現(xiàn)出各向異性，基于此，建立了宏觀應(yīng)變與多邊形方向各向異性指標(biāo)間的關(guān)系。

雖然宏觀應(yīng)力和應(yīng)變已經(jīng)與細(xì)觀接觸力和孔隙多邊形各向異性建立起了聯(lián)系，但這些結(jié)論均基于對(duì)某種試驗(yàn)數(shù)據(jù)的唯象總結(jié)，缺乏理論依據(jù)。Rothenburg等[6]提出的宏觀應(yīng)力比與接觸力各向異性指標(biāo)的聯(lián)系僅適用于三軸壓縮試驗(yàn)，而對(duì)于直剪試驗(yàn)則不適用。針對(duì)這個(gè)問題，基于已有研究，現(xiàn)采用理論推導(dǎo)和數(shù)值試驗(yàn)驗(yàn)證相結(jié)合的方法建立通用的粒狀砂性土宏觀應(yīng)力與細(xì)觀接觸力各向異性指標(biāo)的關(guān)系。以期為構(gòu)建砂土宏細(xì)觀本構(gòu)模型提供新途徑。

1 顆粒流離散元方法

顆粒流離散元法以粒間接觸為計(jì)算單元，基于力-位移法則和牛頓第二運(yùn)動(dòng)定律，采用顯式積分算法模擬顆粒在粒間接觸力作用下的運(yùn)動(dòng)。采用Itasca[8]公司的顆粒流離散元軟件PFC3D進(jìn)行計(jì)算，計(jì)算流程如圖1所示。

圖1 PFC計(jì)算流程簡(jiǎn)圖Fig.1 PFC calculation flow diagram

(1)

式(1)中：d1和d2為相互接觸的兩個(gè)顆粒的半徑，m；G和ν為接觸的剪切模量和泊松比，參考Mitchell等[3]總結(jié)的硅質(zhì)砂的力學(xué)指標(biāo)范圍，取90 GPa和0.2；U為接觸的重疊量，m；Fn為法向接觸力，N。由式(1)可知，接觸的剛度kn和ks是與接觸重疊量相關(guān)的，重疊量越大，剛度越大。另外，接觸法向力與

I和J為兩個(gè)相互接觸顆粒的代號(hào)；l(IJ)為相互接觸的顆粒形心距；δ(IJ)為兩顆粒在接觸處的重疊量；及為X軸和Y軸方向上 的單位向量圖2 Hertz接觸模型簡(jiǎn)圖Fig.2 Schematic diagram of Hertz contact model

接觸切向力滿足Mohr-Coulomb準(zhǔn)則：

Fs≤μFn

(2)

式(2)中：μ為顆粒的摩擦系數(shù)。

2 數(shù)值模型的建立

采用顆粒流離散元軟件PFC3D進(jìn)行數(shù)值建模，在邊長(zhǎng)為1 m的立方體計(jì)算空間中生成5 000個(gè)粒徑為0.01 m的球形顆粒，顆粒摩擦系數(shù)設(shè)為零，由于顆粒粒徑很小，懸浮在空間中，互不接觸。計(jì)算空間的邊界為周期性邊界，當(dāng)顆粒形心越過邊界后，顆粒不會(huì)消失，會(huì)以相同的速度出現(xiàn)在空間中中心對(duì)稱的位置處，這種設(shè)置消除了剛性邊界條件的影響，用于模擬土粒真實(shí)的受力環(huán)境。采用自定義的功能函數(shù)，持續(xù)放大顆粒的粒徑，直到空間中顆粒的體積分?jǐn)?shù)超過52.36%時(shí)停止。等粒徑的顆粒在三維空間中能堆積出的最松散結(jié)構(gòu)為正六面體結(jié)構(gòu)(六面體的每個(gè)節(jié)點(diǎn)位置上有一個(gè)顆粒)，對(duì)應(yīng)的最大的孔隙率是47.64%，而在三維空間中能堆積出的最密實(shí)結(jié)構(gòu)為正四面體結(jié)構(gòu)(四面體的每個(gè)節(jié)點(diǎn)位置上有一個(gè)顆粒)，對(duì)應(yīng)的最小的孔隙率是25.95%。所以當(dāng)顆粒的體積分?jǐn)?shù)達(dá)到52.36%時(shí)，理論上對(duì)應(yīng)空間中最松散狀態(tài)時(shí)的孔隙率，此時(shí)改變顆粒的摩擦系數(shù)為0.5，模型如圖3所示。在最松散的孔隙率下，由于顆粒未能形成正六面體結(jié)構(gòu)，因此，顆粒系統(tǒng)處于流動(dòng)狀態(tài)(unjammed state)。繼續(xù)放大顆粒粒徑或減小計(jì)算空間，會(huì)增加平均顆粒配位數(shù)，使顆粒系統(tǒng)由流動(dòng)狀態(tài)變?yōu)樽枞麪顟B(tài)(jammed state)。自然界中，土體在未破壞之前常處于阻塞狀態(tài)。

圖3 模型圖Fig.3 Model diagram

模型建好后，移動(dòng)計(jì)算空間的底面以模擬單軸壓縮試驗(yàn)，具體步驟為：使計(jì)算空間的底面按照恒定的速度v向上(Z方向)運(yùn)動(dòng)，而保持其他面不動(dòng)，壓縮過程中輸出粒間接觸力信息，計(jì)算試樣的宏觀應(yīng)力，當(dāng)Z方向上的總應(yīng)變量達(dá)到3%時(shí)停止。為了使顆粒體在壓縮過程中保持準(zhǔn)靜態(tài)，v不能太大，取值0.005 m/s。

3 宏觀應(yīng)力張量的表示方法

(3)

在剪切或壓縮過程中，接觸切向力相對(duì)于接觸法向力來(lái)說(shuō)，可以忽略不計(jì)(數(shù)量級(jí)倍的差別)。因此，簡(jiǎn)化平均宏觀應(yīng)力張量，式(3)可轉(zhuǎn)變?yōu)?/p>

(4)

為了證明此結(jié)論的可靠性，分別采用式(3)和式(4)獲取單軸壓縮試驗(yàn)過程中的主應(yīng)力，將式(4)和式(3)的差值與式(3)作比定義為相對(duì)誤差，用Rerror表示。同時(shí)，設(shè)置多種粒間摩擦系數(shù)，μ= 0.2～0.5，以驗(yàn)證粒間摩擦系數(shù)是否對(duì)該結(jié)論產(chǎn)生影響。

由圖4可知，忽略接觸切向力以后，簡(jiǎn)化的平均宏觀應(yīng)力張量與簡(jiǎn)化前的差別甚小，在±6%以內(nèi)，且摩擦系數(shù)越大，相對(duì)誤差就越大。

圖4 簡(jiǎn)化前和簡(jiǎn)化后平均宏觀應(yīng)力張量的相對(duì)值Fig.4 The relative value of the average macroscopic stress tensor before and after simplification

對(duì)于單分散顆粒體(即構(gòu)成顆粒體的顆粒具有相同的直徑)來(lái)說(shuō)，相互接觸的兩顆粒形心的距離l(I,J)為

l(I,J)=Dp-δ(I,J)

(5)

式(5)中：Dp為顆粒的直徑，m。由于接觸的剛度參數(shù)非常大，因此顆粒間的接觸重疊量相對(duì)于顆粒直徑來(lái)說(shuō)，可以忽略不計(jì)，因此可以得

l(I,J)=Dp

(6)

將式(6)代入式(4)可得

(7)

顆粒材料的孔隙率和配位數(shù)之間存在聯(lián)系[9-10]，式(7)中V可采用顆粒配位數(shù)Z、顆粒體積分?jǐn)?shù)ns以及顆粒接觸數(shù)量Nc進(jìn)行換算：

(8)

將式(8)代入式(7)可得

(9)

而在Hertz接觸模型中，接觸法向力的計(jì)算公式為

(10)

式(10)中：n為接觸法向單位向量；E*和R*分別為接觸的彈性模量和相對(duì)半徑，計(jì)算公式為

(11)

式(11)中：EI、EJ、νI、νJ分別為相互接觸的兩顆粒的彈性模量和泊松比；RI和RJ分別為相互接觸的兩顆粒的半徑。

(12)

4 細(xì)觀組構(gòu)張量的表述方法

在荷載作用下，顆粒材料的粒間接觸力分布具有不均勻性，即各向異性，是顆粒材料的特點(diǎn)之一。對(duì)于天然土體，土粒的形狀是隨機(jī)不規(guī)則的，在這種情況下，顆粒的擇優(yōu)定向、分支向量(連接兩個(gè)相互接觸的顆粒形心的向量)、孔隙定向和接觸力的分布均是不均勻分布的，然而采用統(tǒng)計(jì)學(xué)方法可以描述它們的各向異性特點(diǎn)[4]。組構(gòu)一詞常表述的是顆粒的空間排列特點(diǎn)，在這里定義為接觸力分布、分支向量和顆粒定向等在顆粒體中的排列特征。組構(gòu)張量是組構(gòu)的量化形式，常采用二階張量形式進(jìn)行描述[11-12]，即

(13)

(14)

圖5展示了在軸向壓縮試驗(yàn)中獲取的接觸力法向組構(gòu)張量，分別采用特征值之差與特征值之比描述接觸法向各向異性在壓縮過程中的變化，可見，兩種定義都可以很好地描述接觸法向各向異性。出于簡(jiǎn)化的目的，采用式(14)表征組構(gòu)張量的各向異性程度。

圖5 Φ1與Φ3差值和比值在單軸壓縮過程中的變化Fig.5 The difference and ratio of Φ1 to Φ3 during uniaxial compression

5 假定

為建立宏觀應(yīng)力與粒間接觸力各向異性的聯(lián)系，需作適當(dāng)?shù)募僭O(shè)。假設(shè)在單軸壓縮試驗(yàn)中，與壓應(yīng)力作用面平行的方向上(x-y水平面)，即橫向方向上，接觸力向量的大小及數(shù)量分布相同，則Φ2=Φ3，且有Φ1+Φ2+Φ3=1，如圖6所示。

圖6 單軸壓縮試驗(yàn)過程中接觸力法向組構(gòu) 張量的特征值Fig.6 The eigenvalues of the normal fabric tensor of the contact force during the uniaxial compression test

6 平均宏觀應(yīng)力與細(xì)觀組構(gòu)張量的聯(lián)系

宏觀應(yīng)力比σr和細(xì)觀組構(gòu)比Φr分別采用式(15)和式(16)進(jìn)行計(jì)算。式(15)由式(12)推導(dǎo)得出，式(16)由式(13)進(jìn)行推導(dǎo)。

(15)

式(15)中：σ11、σ22和σ33均為宏觀主應(yīng)力；∑表示基于統(tǒng)計(jì)樣本中總接觸數(shù)量的求和。

(16)

式(16)中：Φ11、Φ22和Φ33為接觸法向在三維空間中3個(gè)主方向上的各向異性。

為求出宏觀應(yīng)力比與細(xì)觀組構(gòu)比的關(guān)系。首先定義兩個(gè)無(wú)量綱參數(shù), 令

(17)

圖7為在改變粒間摩擦系數(shù)μ和試樣體積分?jǐn)?shù)ns的情況下，單軸壓縮試樣過程中Rσ和RΦ的變化。其中μ取值為0.0、0.01、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9，ns取值范圍為0.594～0.64。圖7結(jié)果表明，Rσ和RΦ具有很好的線性關(guān)系(線性擬合的相關(guān)性系數(shù)R2超過0.96)，且這種線性關(guān)系與試樣中顆粒的體積分?jǐn)?shù)ns有關(guān)，而與粒間摩擦系數(shù)無(wú)關(guān)μ。采用下式表征Rσ和RΦ的線性關(guān)系：

RΦ=c1Rσ+c2

(18)

式(18)中：c1和c2為與顆粒體積分?jǐn)?shù)相關(guān)的參數(shù)?？梢?，顆粒的體積分?jǐn)?shù)決定著顆粒體的相對(duì)密度，而相對(duì)密度是粒狀土宏細(xì)觀力學(xué)響應(yīng)的主要影響因素之一。

圖7 不同粒間摩擦系數(shù)和顆粒體積分?jǐn)?shù)情況下 Rσ和RΦ在壓縮過程中的變化Fig.7 Changes of Rσ and RΦ during compression under different inter-particle friction coefficients and particle volume fractions

同樣，基于不同體積分?jǐn)?shù)的單軸壓縮試驗(yàn)結(jié)果，擬合得到c1和c2與體積分?jǐn)?shù)的關(guān)系式：

(19)

需要注意的是，Rσ和RΦ在單軸壓縮試驗(yàn)中是隨著顆粒體積分?jǐn)?shù)在變化的，而式(19)中c1和c2是與試樣初始體積分?jǐn)?shù)相關(guān)的。

根據(jù)式(14)，細(xì)觀組構(gòu)各向異性指標(biāo)計(jì)算公式為

(20)

由式(20)可得，細(xì)觀接觸力各向異性與宏觀平均應(yīng)力的關(guān)系為

(21)

7 結(jié)論

土體等顆粒材料的宏觀力學(xué)響應(yīng)與顆粒層面細(xì)觀力學(xué)行為之間的聯(lián)系是現(xiàn)代土力學(xué)中的難點(diǎn)和熱點(diǎn)之一。針對(duì)已有的宏細(xì)觀本構(gòu)關(guān)系多缺乏普適性的問題，根據(jù)已有理論研究成果，采用赫茲接觸理論，基于顆粒流離散元程序PFC3D開展顆粒體宏觀應(yīng)力張量與細(xì)觀接觸力組構(gòu)張量關(guān)系研究，提出了更具普適性的宏細(xì)觀本構(gòu)模型。首先對(duì)宏觀應(yīng)力張量進(jìn)行簡(jiǎn)化，簡(jiǎn)化后的宏觀應(yīng)力張量與實(shí)際值的誤差小于6%。然后定義了細(xì)觀接觸力組構(gòu)各向異性指標(biāo)?；诖耍诩僭O(shè)單軸壓縮橫向方向上組構(gòu)各向異性相同的情況下，定義了兩個(gè)無(wú)量綱的參數(shù)Rσ和RΦ。多組單軸壓縮試驗(yàn)結(jié)果表明，Rσ和RΦ具有明顯的線性關(guān)系，由它們的線性關(guān)系得到了宏觀應(yīng)力張量與細(xì)觀接觸力組構(gòu)張量的關(guān)系。研究成果為聯(lián)系土體宏觀應(yīng)力與細(xì)觀粒間接觸力的關(guān)系提供了參考。

需要注意的是，研究結(jié)果是基于單分散顆粒體，在周期性邊界條件下，單軸壓縮試驗(yàn)中論證的，對(duì)于多分散顆粒體和更復(fù)雜的邊界情況，則需要進(jìn)一步研究。