離散元微觀參數對砂土宏觀參數的影響

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

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

0 引 言

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

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

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

1 雙軸實驗設計

1.1 接觸模型選擇

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

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

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

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

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

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

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

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

1.2 成樣參數及試驗方案

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

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

表1 試樣參數

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

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

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

表2 雙軸實驗方案

圖3 松、密沙的應力—應變關系曲線（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) 維持摩擦系數不變，在保證kn/ks=1.5下改變法向接觸剛度，在100 kPa，200 kPa和400 kPa圍壓下剪切；

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

2 模擬結果分析

2.1 應力應變關系

密樣和松樣的典型應力應變關系和體變關系如圖3、4所示。

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

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

2.2 接觸剛度變化對宏觀參數的影響

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

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

圖6 內摩擦角—法向接觸剛度關系曲線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中可以看到，剛度變化對密樣摩擦角的影響表現為隨剛度增大內摩擦角逐漸減小的特點，從23.6°減小到21.5°，減小了9%，但整體上減小不多，內摩擦角受剛度影響很小。

對于松樣，從圖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可以看到，松樣的內摩擦角隨剛度增大從15.3°增大到15.6°，增幅很小，內摩擦角基本不受剛度變化的影響。

由以上結果看到，顆粒接觸剛度對試樣彈性模量的影響在100 kPa和200 kPa低圍壓下很小，在400 kPa高圍壓下影響較明顯。這里需要指出的是，以上結論是在考慮了顆粒重疊量的情況下得到的。在離散元數值模擬中，顆粒允許發生一定的重疊量，但重疊量必須控制在允許的范圍內，因為實際砂土顆粒間的擠壓變形量非常小，故粒間接觸剛度不能取得太小。為分析重疊量與接觸剛度間的關系，在法向剛度分別為 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所示的顆粒最大重疊量與法向接觸剛度間的關系。圖中最大重疊量是指試樣中相鄰兩顆粒重疊區域寬度最大值與平均粒徑（7.5 mm）的比值。

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

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

2.3 粒間摩擦系數變化對宏觀參數的影響

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

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

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

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

圖9 內摩擦角—摩擦系數關系曲線Fig.9 Friction angle-frictional coefficientrelationship curves for both loose and dense sands.

2.4 離散元數值計算參數選取建議

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

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

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

3 結論

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

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

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

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

(4) 粒間摩擦系數對密樣和松樣的內摩擦角均有非常顯著影響。

[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]焦紅光，李靖如，趙繼芬，等. 關于離散元法計算參數的探討[J]. 河南理工大學學報(自然科學版), 2007, 26(1): 88-93.

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

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

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

[9]徐小敏，凌道盛，陳云敏，等.基于線性接觸模型的顆粒材料細–宏觀彈性常數相關關系研究[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.