申志福，蔣明鏡，朱方園，胡海軍

(1.同濟(jì)大學(xué)地下建筑與工程系，上海 200092；2.同濟(jì)大學(xué)巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092）

0 引 言

離散單元法是上世紀(jì)70年代由Cundall[1]提出并逐步發(fā)展起來(lái)的一種數(shù)值計(jì)算方法，它將介質(zhì)看做一系列離散單元體的組合，單元體運(yùn)動(dòng)受牛頓經(jīng)典力學(xué)定律控制。離散單元法分析中假設(shè)顆粒單元為剛性體，接觸發(fā)生在很小范圍內(nèi)，接觸處允許有一定重疊量，重疊量大小與接觸力和接觸剛度有關(guān)。顆粒間的接觸和相互作用關(guān)系通過(guò)粒間接觸模型控制。該方法現(xiàn)已大量運(yùn)用于巖土體數(shù)值分析。為了研究不同特性的巖土體，各國(guó)學(xué)者提出了多種顆粒接觸模型，如理想線彈性接觸模型，Hertz-Mindlin非線性模型[2]，Iwashita K等[3]提出的抗轉(zhuǎn)動(dòng)模型，蔣明鏡等[4]提出的考慮毛細(xì)水的接觸模型等。顆粒接觸剛度和粒間摩擦系數(shù)是這些模型中兩個(gè)最基本的參數(shù)，也是離散元參數(shù)標(biāo)定過(guò)程中的主要調(diào)整參數(shù)。

從離散元計(jì)算過(guò)程來(lái)看，顆粒接觸剛度大小決定著計(jì)算時(shí)步。計(jì)算時(shí)步太大，離散元計(jì)算精度下降；時(shí)步太小結(jié)果波動(dòng)大且計(jì)算時(shí)間延長(zhǎng)[5]。從計(jì)算結(jié)果來(lái)看，顆粒接觸剛度和粒間摩擦系數(shù)對(duì)土體宏觀力學(xué)參數(shù)如初始切線模量、內(nèi)摩擦角有著很大影響，只有恰當(dāng)?shù)奈⒂^參數(shù)取值才可能得到與實(shí)際接近的數(shù)值計(jì)算結(jié)果。因此合理地選定微觀參數(shù)對(duì)離散元模擬的效率和結(jié)果可信度都至關(guān)重要。

目前，離散元微觀參數(shù)的選取有多種方法。武力等[6]將遺傳神經(jīng)網(wǎng)絡(luò)與三軸試驗(yàn)離散元數(shù)值模擬有機(jī)結(jié)合，用于改性砂土顆粒離散元接觸模型參數(shù)反演；邢紀(jì)波等[7]根據(jù)離散單元中的應(yīng)力波傳播條件確定粒間接觸剛度的理論公式；周健等[8]針對(duì)顆粒流仿真試樣的細(xì)觀力學(xué)特征與物理試樣宏觀力學(xué)響應(yīng)之間的關(guān)系進(jìn)行了參數(shù)研究,并揭示了一些規(guī)律；徐小敏等[9]利用顆粒材料單元體宏觀力學(xué)參數(shù)和顆粒細(xì)觀參數(shù)間的相關(guān)性，通過(guò)室內(nèi)三軸試驗(yàn)的模擬和結(jié)果的回歸分析，建立了宏微觀參數(shù)間的經(jīng)驗(yàn)公式。本文運(yùn)用二維離散元商業(yè)軟件Particle Flow Code（PFC2D）進(jìn)行66組雙軸壓縮實(shí)驗(yàn)以探討兩種不同孔隙比（分別代表密砂和松砂）試樣的粒間接觸剛度和摩擦系數(shù)對(duì)試樣初始切線模量（后文均稱(chēng)為彈性模量）和內(nèi)摩擦角的影響。

1 雙軸實(shí)驗(yàn)設(shè)計(jì)

1.1 接觸模型選擇

因本文只針對(duì)顆粒接觸剛度和粒間摩擦系數(shù)兩個(gè)微觀參數(shù)對(duì)宏觀參數(shù)的影響規(guī)律進(jìn)行分析，不針對(duì)特定土體進(jìn)行參數(shù)標(biāo)定，故接觸模型采用簡(jiǎn)單的線彈性接觸模型，不考慮粒間膠結(jié)作用。在 PFC中，顆粒間線性接觸模型如圖1所示：

圖1 粒間線性接觸模型Fig.1 Interparticle linear contact modle of PFC20.

粒間法向接觸力采用全量法計(jì)算，由相對(duì)法向位移和法向接觸剛度決定：

式中 Fn為法向接觸力；kn為法向接觸剛度；Un為法向相對(duì)位移。

粒間切向接觸力采用增量法計(jì)算，由相對(duì)切向位移增量和切向接觸剛度決定：

式中ΔFs為切向接觸力增量；ks為切向接觸剛度；ΔUs為切向相對(duì)位移增量。

當(dāng)顆粒間無(wú)膠結(jié)時(shí)，切向力不能大于顆粒間摩擦系數(shù)與法向力的乘積：

式中Fs為切向接觸力；μ為粒間摩擦系數(shù)。

1.2 成樣參數(shù)及試驗(yàn)方案

雙軸試驗(yàn)的模擬過(guò)程分為成樣，固結(jié)和壓縮三部分。成樣方法采用蔣明鏡等[10]提出的分層欠壓法。為使試樣整體均勻，該方法將試樣分層生成，新的一層顆粒生成后試樣的平均孔隙比較前幾層平均孔隙比小，最后一層生成后試樣整體達(dá)到目標(biāo)孔隙比。該方法通過(guò)調(diào)節(jié)試樣孔隙比可形成松砂、中密砂及密砂試樣，同時(shí)能避免成樣過(guò)程中出現(xiàn)上松下密的現(xiàn)象。

本實(shí)驗(yàn)共生成兩種孔隙比試樣：0.21和0.25，分別代表密砂和松砂。密砂樣和松砂樣顆?？倲?shù)目均為6 000,直徑 6～9 mm,顆粒密度均為 2 600 kg·m-3。為能達(dá)到目標(biāo)孔隙比，密樣成樣過(guò)程中粒間摩擦系數(shù)為0.1，松樣為1.0，顆粒與墻之間的摩擦系數(shù)均為 0。成樣過(guò)程中采用的最優(yōu)孔隙比如下：密樣為ep(1)=0.233， ep(1+2)=0.228， ep(1+2+3)=0.225,ep(1+2+3+4)=0.219，ep(1+2+3+4+5)=0.21；松樣為ep(1)=0.27， ep(1+2)=0.269，ep(1+2+3)=0.265,ep(1+2+3+4)=0.259，ep(1+2+3+4+5)=0.25。離散元雙軸壓縮試驗(yàn)成樣過(guò)程參數(shù)見(jiàn)表1。生成試樣級(jí)配曲線如圖2所示。

表1 試樣參數(shù)

圖2 顆粒級(jí)配圖Fig.2 Distribution of grain size.

成樣后分別在100 kPa，200 kPa和400 kPa圍壓下等向固結(jié)；固結(jié)后豎向加壓；以10%/min的速率剪切；采用伺服方式保持圍壓穩(wěn)定；實(shí)驗(yàn)采用剛性邊界。在剪切過(guò)程中記錄試樣的軸向應(yīng)力、軸向應(yīng)變、體變和孔隙比。

為了探究微觀參數(shù)（顆粒接觸剛度和顆粒間摩擦系數(shù)）對(duì)宏觀參數(shù)（彈性模量和內(nèi)摩擦角）的影響，剪切過(guò)程按表2方案進(jìn)行：

表2 雙軸實(shí)驗(yàn)方案

圖3 松、密沙的應(yīng)力—應(yīng)變關(guān)系曲線（kn=1.5×108 N·m-1，ks=1.0×108 N·m-1）Fig.3 Stress-strain relationship curves of loose and dense sands.

(1) 維持摩擦系數(shù)不變，在保證kn/ks=1.5下改變法向接觸剛度，在100 kPa，200 kPa和400 kPa圍壓下剪切；

(2) 維持粒間接觸剛度不變，改變粒間摩擦系數(shù)，在100 kPa，200 kPa和400 kPa圍壓下剪切，其中密樣μ變化范圍為0.1～0.7，松樣μ變化范圍為0.1～0.5。

2 模擬結(jié)果分析

2.1 應(yīng)力應(yīng)變關(guān)系

密樣和松樣的典型應(yīng)力應(yīng)變關(guān)系和體變關(guān)系如圖3、4所示。

圖4 松、密砂的體變—軸向應(yīng)變關(guān)系曲線（kn=1.5×108 N·m-1，ks=1.0×108 N·m-1）Fig.4 Volumetric strain-axial strain relationship curves of loose and dense sands.

從圖3（a）可見(jiàn)，密樣的應(yīng)力應(yīng)變關(guān)系為軟化型：不同圍壓下的應(yīng)力均先隨軸向應(yīng)變的增大而增大，在軸向應(yīng)變?yōu)?%～2%時(shí)達(dá)到強(qiáng)度峰值后迅速減小，最終趨于穩(wěn)定。由于顆粒較少，在400 kPa高圍壓下峰值后應(yīng)力波動(dòng)較明顯。從圖4（a）可見(jiàn)，各圍壓下密樣的體變隨軸向應(yīng)變的變化均表現(xiàn)出先剪縮后剪脹的特性，即在軸向應(yīng)變?yōu)?%～2%時(shí)達(dá)到剪縮最大值后剪縮量逐漸減少進(jìn)而表現(xiàn)出剪脹，剪脹體變隨軸向應(yīng)變?cè)龃蠖鴿u趨穩(wěn)定。這些都是密砂的典型力學(xué)特性。

從圖3（b）可見(jiàn)，松樣的應(yīng)力應(yīng)變關(guān)系為硬化型：各圍壓下應(yīng)力均隨軸向應(yīng)變?cè)龃蠖龃螅鲩L(zhǎng)趨勢(shì)漸趨平緩。由于顆粒較少，在400 kPa高圍壓下應(yīng)力波動(dòng)較明顯。從圖4（b）可見(jiàn)，各圍壓下體變均表現(xiàn)為剪縮，剪縮體變隨軸向應(yīng)變?cè)龃蠖龃笄以鲩L(zhǎng)趨勢(shì)漸緩。這些都是松砂的典型力學(xué)特性。

2.2 接觸剛度變化對(duì)宏觀參數(shù)的影響

按照表2中第1組試驗(yàn)方案，維持摩擦系數(shù)不變，在保證 kn/ks=1.5情況下改變法向接觸剛度進(jìn)行雙軸試驗(yàn)得到試樣彈性模量、內(nèi)摩擦角與法向接觸剛度之間的關(guān)系，如圖5、6所示。

圖5 彈性模量—法向接觸剛度關(guān)系曲線Fig.5 Young’s modulus-normal stiffness relationship curves of loose and dense sands.

圖6 內(nèi)摩擦角—法向接觸剛度關(guān)系曲線Fig.6 Friction angle-normal stiffness relationship curve.

對(duì)于密樣，從圖5(a)中可以看到，在相同剛度下，圍壓越大彈性模量越大；在100 kPa和200 kPa較低圍壓下接觸剛度變化對(duì)試樣彈性模量的影響并不明顯；在400 kPa高圍壓下，當(dāng)接觸剛度從1.5×108N·m-1變化到1.5×109N·m-1時(shí)，彈性模量從42.9 MPa增長(zhǎng)到84 MPa,增長(zhǎng)了近一倍。從圖6中可以看到，剛度變化對(duì)密樣摩擦角的影響表現(xiàn)為隨剛度增大內(nèi)摩擦角逐漸減小的特點(diǎn)，從23.6°減小到21.5°，減小了9%，但整體上減小不多，內(nèi)摩擦角受剛度影響很小。

對(duì)于松樣，從圖5(b)中可以看到，在相同剛度下，圍壓越大彈性模量越大；與密樣相似，在100 kPa和200 kPa較低圍壓下彈性模量基本不受剛度變化的影響；在400 kPa高圍壓下試樣彈性模量隨接觸剛度增大而有明顯增大，當(dāng)接觸剛度從1.5×108N·m-1變化到1.5×109N·m-1時(shí)，彈性模量從23.8 MPa增長(zhǎng)到33.2 MPa，增大了40%。從圖6可以看到，松樣的內(nèi)摩擦角隨剛度增大從15.3°增大到15.6°，增幅很小，內(nèi)摩擦角基本不受剛度變化的影響。

由以上結(jié)果看到，顆粒接觸剛度對(duì)試樣彈性模量的影響在100 kPa和200 kPa低圍壓下很小，在400 kPa高圍壓下影響較明顯。這里需要指出的是，以上結(jié)論是在考慮了顆粒重疊量的情況下得到的。在離散元數(shù)值模擬中，顆粒允許發(fā)生一定的重疊量，但重疊量必須控制在允許的范圍內(nèi)，因?yàn)閷?shí)際砂土顆粒間的擠壓變形量非常小，故粒間接觸剛度不能取得太小。為分析重疊量與接觸剛度間的關(guān)系，在法向剛度分別為 1.5×107N·m-1、7.5×107N·m-1、1.5×108N·m-1、4.5×108N·m-1和7.5×108N·m-1五種情況下進(jìn)行密砂雙軸試驗(yàn)得出了如圖7所示的顆粒最大重疊量與法向接觸剛度間的關(guān)系。圖中最大重疊量是指試樣中相鄰兩顆粒重疊區(qū)域?qū)挾茸畲笾蹬c平均粒徑（7.5 mm）的比值。

從圖7可以看出，在相同剛度下，圍壓越大，粒間最大重疊量越大；在相同圍壓下，隨法向接觸剛度增大，最大重疊量有顯著減小。本文以 45‰作為重疊量最大允許值，要保證在各圍壓下重疊量都不超過(guò)圖7中的控制線，接觸剛度至少為1.5×108N·m-1。考慮了重疊量影響后，各剛度下的重疊量在100 kPa和200 kPa低圍壓下變化不大，在400 kPa高圍壓下變化較為明顯，這也是上述彈性模量隨接觸剛度變化規(guī)律的微觀解釋。

圖7 重疊量—法向接觸剛度關(guān)系曲線Fig.7 Overlap-normal stiffness relationship curves.

2.3 粒間摩擦系數(shù)變化對(duì)宏觀參數(shù)的影響

按照表2中第2組試驗(yàn)方案，維持粒間接觸剛度不變，改變粒間摩擦系數(shù)進(jìn)行雙軸試驗(yàn)，得到試樣彈性模量、內(nèi)摩擦角與粒間摩擦系數(shù)之間的關(guān)系，如圖8、9所示。

對(duì)于密樣，從圖8(a)中可以看到，各圍壓下試樣彈性模量隨摩擦系數(shù)的增大而增大，增長(zhǎng)趨勢(shì)漸緩；在100 kPa、200 kPa和400 kPa下，當(dāng)摩擦系數(shù)由0.1增加到0.7時(shí)，彈性模量增幅分別為51%、56%和65%，這說(shuō)明粒間摩擦系數(shù)對(duì)彈性模量有較明顯影響。從圖9中可以看到，試樣內(nèi)摩擦角隨粒間摩擦系數(shù)的增大幾乎呈線性增大，從12.1°增加到26°，反應(yīng)了顆粒間摩擦作用對(duì)密砂強(qiáng)度的顯著影響。

對(duì)于松樣，從圖8(b)中可以看到，彈性模量隨摩擦系數(shù)增大而增大，彈性模量的增幅在較大摩擦系數(shù)時(shí)有明顯增大且圍壓越大這種增大的趨勢(shì)越明顯。在100 kPa、200 kPa和400 kPa下，當(dāng)摩擦系數(shù)由0.1增加到0.5時(shí)，彈性模量分別增加了1.24倍、1.66倍和2.55倍，說(shuō)明摩擦系數(shù)對(duì)松樣彈性模量的影響比密樣大得多。從圖9中可以看到，松樣的內(nèi)摩擦角隨摩擦系數(shù)的增大而增大，當(dāng)摩擦系數(shù)較大時(shí)增大的幅度有所減緩；當(dāng)摩擦系數(shù)由0.1增加到0.5時(shí)，內(nèi)摩擦角從10.9°增加到15.3°，說(shuō)明摩擦系數(shù)對(duì)松砂強(qiáng)度有明顯貢獻(xiàn)。

圖8 彈性模量—摩擦系數(shù)關(guān)系曲線Fig.8 Young’s modulus-frictional coefficient relationship curves of loose and dense sands.

圖9 內(nèi)摩擦角—摩擦系數(shù)關(guān)系曲線Fig.9 Friction angle-frictional coefficientrelationship curves for both loose and dense sands.

2.4 離散元數(shù)值計(jì)算參數(shù)選取建議

由以上結(jié)果可以看出，無(wú)論密樣還是松樣，兩個(gè)離散元微觀參數(shù)各自對(duì)宏觀力學(xué)參數(shù)的影響程度是不同的，可形象地用圖10說(shuō)明。微觀參數(shù)對(duì)宏觀參數(shù)的影響程度用箭頭表示，箭頭越大顏色愈深，表示影響程度越大。從圖中可以看到，對(duì)于微觀參數(shù)，粒間摩擦系數(shù)對(duì)宏觀參數(shù)的影響較顆粒接觸剛度產(chǎn)生的影響大；對(duì)于宏觀參數(shù)，初始切線彈性模量受到粒間摩擦系數(shù)影響較受到接觸剛度的影響大，內(nèi)摩擦角則主要由粒間摩擦系數(shù)控制而基本不受接觸剛度影響。

圖10 微觀參數(shù)對(duì)宏觀參數(shù)影響Fig.10 Sketch of influence of micro-parameters on macro-parameters.

在進(jìn)行離散元數(shù)值模擬時(shí)，可通過(guò)圖10的規(guī)律合理選取、調(diào)整微觀參數(shù)。由于接觸剛度對(duì)宏觀參數(shù)影響較粒間摩擦系數(shù)產(chǎn)生的影響小得多，故可考慮調(diào)整接觸剛度的取值，如選取較小的接觸剛度以增大計(jì)算時(shí)步提高效率而不致引起宏觀參數(shù)的太大變化；但剛度不能太小，否則重疊量過(guò)大，影響計(jì)算可信度，建議法向接觸剛度不宜小于1.5×108N·m-1。粒間摩擦系數(shù)因?qū)暧^參數(shù)影響顯著，故當(dāng)模擬結(jié)果與目標(biāo)結(jié)果出入較大時(shí)調(diào)整摩擦系數(shù)將使參數(shù)調(diào)整效果更為明顯有效。

3 結(jié)論

本文運(yùn)用PFC2D通過(guò)66組雙軸壓縮試驗(yàn)對(duì)密砂樣和松砂樣進(jìn)行了離散元模擬，探究了顆粒接觸剛度和粒間摩擦系數(shù)兩個(gè)微觀參數(shù)對(duì)土體初始切線彈性模量和內(nèi)摩擦角的影響，得出以下結(jié)論：

(1) 當(dāng)考慮了顆粒重疊量限制后，在400 kPa較高圍壓下，粒間接觸剛度對(duì)密樣和松樣彈性模量均有一定程度影響且對(duì)密樣影響大于松樣；在100 kPa和200 kPa較低圍壓下密樣和松樣彈性模量均基本不受接觸剛度影響。

(2) 粒間接觸剛度對(duì)密樣和松樣的內(nèi)摩擦角均影響很小。

(3) 粒間摩擦系數(shù)對(duì)密樣和松樣彈性模量都有明顯影響且對(duì)松樣的影響大于密樣。

(4) 粒間摩擦系數(shù)對(duì)密樣和松樣的內(nèi)摩擦角均有非常顯著影響。

[1]Cundall P A. A Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock System[C]//Proceedings of the Symposium of the International Society of Rock Mechanics. Nancy, France: [s.n.],1971.

[2]Mindlin R D, H Deresiewicz. Elastic Spheres in Contact under Varying Oblique Forces[J]. Journal of Applied Mechanics, 1953, 20(3): 327-344.

[3]Iwashita K, Oda M. Rolling Resistance at Contacts in Simulation of Shear Band Development by DEM[J]. Journal of Engineering Mechanics,1998,124(3): 285-292.

[4]Jiang M J, Leroueil S., Konrad J M. Insight into Shear Strength Functions of Unsaturated Granulates by DEM Analyses[J]. Computers and Geotechnics,2004, 31(6):473-489.

[5]焦紅光，李靖如，趙繼芬，等. 關(guān)于離散元法計(jì)算參數(shù)的探討[J]. 河南理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 26(1): 88-93.

[6]武力，屈福政，李守巨. 土壓平衡盾構(gòu)改性砂土離散元模型參數(shù)反演方法研究[J]. 大連理工大學(xué)學(xué)報(bào),2010, 50(6): 860-866.

[7]邢紀(jì)波，俞良群，張瑞豐，等. 離散單元法的計(jì)算參數(shù)和求解方法選擇[J]. 計(jì)算力學(xué)學(xué)報(bào), 1999, 16(1): 47-51, 99.

[8]周健，廖雄華，池永，等.土的室內(nèi)平面應(yīng)變?cè)囼?yàn)的顆粒流模擬[J].同濟(jì)大學(xué)學(xué)報(bào), 2002, 30(9):1044-1050.

[9]徐小敏，凌道盛，陳云敏，等.基于線性接觸模型的顆粒材料細(xì)–宏觀彈性常數(shù)相關(guān)關(guān)系研究[J].巖土工程學(xué)報(bào), 2010,32(7):991-998.

[10]Jiang M J, Konrad J M, Leroueil S. An Efficient Technique for Generating Homogeneous Specimens for DEM Studies[J]. Computers and Geotechnics, 2003, 30(7): 579-597.