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

(1.同濟大學地下建筑與工程系，上海 200092；2.同濟大學巖土及地下工程教育部重點實驗室，上海 200092）

0 引 言

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

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

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

1 雙軸實驗設(shè)計

1.1 接觸模型選擇

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

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

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

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

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

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

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

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

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

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

本實驗共生成兩種孔隙比試樣：0.21和0.25，分別代表密砂和松砂。密砂樣和松砂樣顆粒總數(shù)目均為6 000,直徑 6～9 mm,顆粒密度均為 2 600 kg·m-3。為能達到目標孔隙比，密樣成樣過程中粒間摩擦系數(shù)為0.1，松樣為1.0，顆粒與墻之間的摩擦系數(shù)均為 0。成樣過程中采用的最優(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。離散元雙軸壓縮試驗成樣過程參數(shù)見表1。生成試樣級配曲線如圖2所示。

表1 試樣參數(shù)

圖2 顆粒級配圖Fig.2 Distribution of grain size.

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

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

表2 雙軸實驗方案

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

從圖3（b）可見，松樣的應(yīng)力應(yīng)變關(guān)系為硬化型：各圍壓下應(yīng)力均隨軸向應(yīng)變增大而增大，增長趨勢漸趨平緩。由于顆粒較少，在400 kPa高圍壓下應(yīng)力波動較明顯。從圖4（b）可見，各圍壓下體變均表現(xiàn)為剪縮，剪縮體變隨軸向應(yīng)變增大而增大且增長趨勢漸緩。這些都是松砂的典型力學特性。

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

按照表2中第1組試驗方案，維持摩擦系數(shù)不變，在保證 kn/ks=1.5情況下改變法向接觸剛度進行雙軸試驗得到試樣彈性模量、內(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.

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

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

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

從圖7可以看出，在相同剛度下，圍壓越大，粒間最大重疊量越大；在相同圍壓下，隨法向接觸剛度增大，最大重疊量有顯著減小。本文以 45‰作為重疊量最大允許值，要保證在各圍壓下重疊量都不超過圖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ù)變化對宏觀參數(shù)的影響

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

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

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

圖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ù)值計算參數(shù)選取建議

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

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

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

3 結(jié)論

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

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

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

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

(4) 粒間摩擦系數(shù)對密樣和松樣的內(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)于離散元法計算參數(shù)的探討[J]. 河南理工大學學報(自然科學版), 2007, 26(1): 88-93.

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

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

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

[9]徐小敏，凌道盛，陳云敏，等.基于線性接觸模型的顆粒材料細–宏觀彈性常數(shù)相關(guān)關(guān)系研究[J].巖土工程學報, 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.