不同模型維數(shù)下常規(guī)三軸試驗的離散元研究

周 杰

(重慶交通大學 土木建筑學院，重慶 400074)

0 引 言

常規(guī)三軸試驗是實驗室研究土體力學性質(zhì)重要方法，但是實驗室試驗僅能獲得土體的宏觀力學性質(zhì)，并不能獲得土體的微觀結(jié)構(gòu)參數(shù)，故不能從本質(zhì)上解釋其宏觀力學響應。離散元因能同時獲得數(shù)值試樣的宏、微觀參數(shù)，而在巖土力學性質(zhì)的研究中顯示出很強的生命力[1-5]。羅勇，等[1]采用三維離散元數(shù)值模型對常規(guī)三軸試驗進行離散元數(shù)值研究，認為三維離散元數(shù)值試驗結(jié)果與實驗室試驗結(jié)果有良好的吻合性。雖然實驗室試驗試樣處于三維狀態(tài)，但二維離散元數(shù)值試驗的計算效率快、試樣微觀結(jié)構(gòu)的定量表達簡單，因而得到了廣泛的應用。周健、史旦達、蔣明鏡和Li，等[2-5]借助二維離散元進行三軸試驗或直剪試驗的數(shù)值研究，通過分析二維數(shù)值試樣的宏觀力學行為，結(jié)合接觸法向、顆粒旋轉(zhuǎn)等微觀參數(shù)的演化規(guī)律，深入研究土體宏觀力學性質(zhì)的微觀本質(zhì)。

但是，將實際的三維問題簡化為二維離散元數(shù)值模型，會涉及到以下問題：離散元模型維數(shù)的減少對模擬結(jié)果有多大的影響？通過對比二維離散元數(shù)值試驗結(jié)果與實驗室試驗結(jié)果而確定的模型參數(shù)，能否代表實驗室試樣本身的性質(zhì)？

筆者采用離散元程序——顆粒流PFC——分別進行了實驗室常規(guī)三軸試驗的三維和二維離散元模擬，保證兩種維數(shù)模型的材料參數(shù)和試樣密實度均相同，從試樣宏觀力學響應和微觀結(jié)構(gòu)特征的角度分析模型維數(shù)對數(shù)值結(jié)果的影響。

1 數(shù)值模型與參數(shù)

實驗室常規(guī)三軸試驗的三維和二維離散元數(shù)值分析，即三軸和雙軸離散元數(shù)值分析[6]。三軸數(shù)值試驗可按照試樣的實際尺寸建立圓柱形的數(shù)值試樣，而雙軸試驗的試樣為平面矩形，其寬和高分別為實驗室三維試樣的直徑和高。為了分析試樣維數(shù)對數(shù)值試驗結(jié)果的影響，需保證兩種維數(shù)試樣的尺寸、模型參數(shù)、試樣的孔隙比均相等。

1.1 三維、二維離散元數(shù)值模型的建立

為了保證本研究三維和二維數(shù)值試樣的孔隙比具有可比性，制備最密實的等粒徑數(shù)值試樣，即三維試樣顆粒在空間呈正三角錐形〔圖1(a)〕排列、二維試樣顆粒在平面呈正三角形〔圖1(b)〕排列。

圖1 數(shù)值試樣內(nèi)顆粒的排列結(jié)構(gòu)Fig.1 Particle packing structure in 2D/3D DEM model

1.2 數(shù)值模型參數(shù)的確定

研究所進行的離散元數(shù)值試驗，采用文獻[6]中推薦的模型參數(shù)，如表1。對數(shù)值試樣加載的水平“墻體”剛度kn=ks=1×108N/m，對試樣施加圍壓的豎直“墻體”剛度kn=ks=1×107N/m；顆粒的摩擦系數(shù)μb=0.5，墻體的摩擦系數(shù)μw=0.0。為了消除數(shù)值試樣中應力、應變的離散計算方法而產(chǎn)生的誤差，保證數(shù)值試樣內(nèi)顆粒的數(shù)目大于2 000個[7]。研究中設置數(shù)值試樣直徑方向35個顆粒、高度方向82個顆粒，顆粒半徑由程序根據(jù)試樣尺寸計算確定[8]。

表1 二維和三維離散元模型參數(shù)設置

2 離散元數(shù)值試驗結(jié)果分析

采用伺服控制程序?qū)θS和二維數(shù)值試樣進行剪切[9]，實現(xiàn)對實驗室試驗過程的模擬。

2.1 數(shù)值試樣的宏觀力學性質(zhì)

三維試樣和二維試樣的剪切強度-剪切應變的關系曲線如圖2，由于三維和二維數(shù)值試樣均處于最密實的狀態(tài)，故試樣在剪切過程中發(fā)生明顯的軟化特征；圍壓越高，試樣抗剪強度越大、其所對應的應變也越大。數(shù)值試驗結(jié)果與實驗室試驗結(jié)果的規(guī)律一致。

圖2 數(shù)值試樣的剪切應力-應變關系曲線Fig.2 Shear stress-strain curves of numerical specimens

在相同圍壓的作用下，三維試樣達到強度峰值時的應變明顯小于二維試樣的峰值應變，兩者之比約為1/4。當圍壓較低(<0.5 MPa)時，二維試樣峰值強度大于三維試樣的峰值強度；隨著圍壓的增大，三維試樣的峰值強度逐漸大于二維試樣的峰值強度。對于文中所研究最密實的數(shù)值試樣，不同維數(shù)數(shù)值試樣的峰值強度在數(shù)值上相差不大，其間的差距會隨試樣的松散而變得明顯。三維和二維數(shù)值試樣的黏聚力分別為0.28和0.46 MPa、內(nèi)摩擦角分別為48.7°和43.9°，由于數(shù)值試樣處于理想的最密實狀態(tài)，故數(shù)值試驗所得試樣抗剪強度指標偏大。試樣黏聚力在微觀上與顆粒位置調(diào)整的能力有關，而其內(nèi)摩擦角在微觀上與顆粒間的摩擦作用點的數(shù)目有關。三維數(shù)值試樣內(nèi)顆粒在空間位置調(diào)整的可能性更大，故宏觀上表現(xiàn)為三維試樣的內(nèi)黏聚力較??；而三維試樣內(nèi)顆粒間的摩擦作用點較多，在圍壓作用下顆粒間的摩擦力使宏觀上表現(xiàn)為其內(nèi)摩擦角較大。三維和二維數(shù)值試樣在不同圍壓條件下宏觀力學響應的差異，本質(zhì)上都是由于試樣內(nèi)部結(jié)構(gòu)的不同、結(jié)構(gòu)對外力的抵抗能力所決定的：試樣內(nèi)顆粒之間的約束越多、位置調(diào)整越難，則其宏觀的抗剪強度越大。

由以上的分析可知，不同維數(shù)數(shù)值試樣宏觀力學行為的差異，本質(zhì)上反映了其微觀結(jié)構(gòu)特征的不同。

2.2 數(shù)值試樣的微觀結(jié)構(gòu)演化

圖3為剪切破壞后(圍壓0.1 MPa)三維和二維數(shù)值試樣的接觸力鏈和位移場，由圖3可知，三維數(shù)值試樣的剪切帶呈“X型”，剪切帶內(nèi)顆粒之間的接觸力和位移均較?。辉嚇由?、下邊界上或左、右邊界上的接觸力和位移基本相同，以“X型”剪切帶為分界線，試樣上、下邊界附近的顆粒發(fā)生指向試樣內(nèi)側(cè)的位移，而左、右邊界的顆粒發(fā)生指向試樣外側(cè)的位移。二維試樣的剪切帶呈“XX型”，試樣核心部位顆粒的位移較?。煌贿吔绮煌课活w粒的接觸力和位移明顯不同，僅試樣上、下邊界中心部位的顆粒發(fā)生指向試樣內(nèi)側(cè)的位移，其他部位的顆粒均發(fā)生指向試樣外側(cè)的位移。

為了便于定量比較三維和二維數(shù)值試樣內(nèi)接觸法向ni的分布，筆者僅統(tǒng)計三維試樣在x-z平面內(nèi)的接觸信息。定義Δθ角度范圍內(nèi)的接觸法向密度EC(ni)為：

(1)

式中：ΔN為接觸法向ni落在Δθ范圍內(nèi)的數(shù)目；N為試樣內(nèi)接觸的總數(shù)目。

將研究平面內(nèi)數(shù)值試樣顆粒接觸法向的分布表示在圖4中，由于最密實三維和二維試樣內(nèi)顆粒均呈正三角形排列，故剪切后試樣的接觸法向主方向為0°，60°，120°；三維數(shù)值試樣內(nèi)顆粒的排列在受剪切作用的影響，剪切后少量顆粒接觸法向調(diào)整為50°，90°，130°方向，三維數(shù)值試樣內(nèi)接觸法向的分布較二維數(shù)值試樣內(nèi)接觸法向的分布均勻。

圖4 破壞后數(shù)值試樣的接觸法向密度分布規(guī)律(圍壓0.1 MPa)

剪切后三維和二維數(shù)值試樣的切向接觸力分布規(guī)律如圖5，由于三維數(shù)值試樣采用球形顆粒，而二維數(shù)值試樣采用單位厚度的圓盤顆粒，數(shù)值試樣的接觸面積不同，要使試樣達到相同的圍壓水平，二維數(shù)值試樣需要更大的外力。因此筆者關于試樣內(nèi)接觸力的分析，僅關注兩種維數(shù)數(shù)值試樣接觸力的分布規(guī)律。

圖5 破壞后數(shù)值試樣的切向接觸力分布規(guī)律(圍壓0.1 MPa)

三維試樣內(nèi)顆粒間的切向接觸力集中在25～75°和115～145°之間，這與三維數(shù)值試樣宏觀上所表現(xiàn)出來的“X型”剪切帶相一致。二維數(shù)值試樣切向接觸力主要集中在4個角度區(qū)間：15～25°、75～85°、105～115°和135～155°，故宏觀上二維試樣的剪切帶為“XX 型”。試樣內(nèi)切向接觸力的微觀分布規(guī)律與宏觀剪切帶的一致性，進一步證明了試樣微觀結(jié)構(gòu)對宏觀響應的決定性作用。

當試驗圍壓由0.1 MPa增大至1 MPa時，三維試樣的剪切帶由“X型”變?yōu)椤啊笮汀?，二維試樣的剪切帶也由“XX型”變?yōu)椤癤型”。圍壓增大，試樣剪切帶的范圍趨于集中(圖6)。

圖6 破壞后數(shù)值試樣的接觸力鏈和位移場(圍壓1 MPa)

分析圖7所示試樣內(nèi)切向接觸力的分布規(guī)律可知，圍壓增大至1 MPa后，三維試樣內(nèi)顆粒間的切向接觸力集中在45～65°和115～145°范圍內(nèi)，二維數(shù)值試樣內(nèi)顆粒間的切向接觸力集中在35～75°、105～125°和155～165°范圍內(nèi)。圍壓增大使二維和三維試樣的切向接觸力的分布更加集中，其方向接近接觸法向的分布方向(0°，60°，120°)上，與宏觀上試樣剪切帶的變化規(guī)律一致。

圖7 破壞后數(shù)值試樣的切向接觸力分布規(guī)律(圍壓1 MPa)

3 結(jié) 論

針對實驗室常規(guī)三軸試驗，進行了最密實狀態(tài)的、參數(shù)相同的三維和二維的離散元數(shù)值分析，通過分析數(shù)值試樣的宏觀力學響應和微觀結(jié)構(gòu)特征，獲得了以下結(jié)論：

1)最密實三維和二維數(shù)值試樣在剪切過程中均表現(xiàn)出軟化的力學性質(zhì)，隨圍壓的增大，不同維數(shù)試樣峰值強度增大，三維試樣的峰值強度逐漸大于二維試樣的峰值強度；不同圍壓下，三維試樣峰值強度所對應的應變一直約為二維試樣的1/4倍。

2)由于數(shù)值試樣處于理想的最密實狀態(tài)，故數(shù)值試驗所得試樣抗剪強度指標偏大：三維和二維數(shù)值試樣的黏聚力分別為0.28和0.46 MPa、內(nèi)摩擦角分別為48.7°和43.9°。試樣黏聚力在微觀上與顆粒位置調(diào)整的能力有關，內(nèi)摩擦角在微觀上與顆粒間的摩擦點的作用有關。

3)在相同圍壓下，三維和二維試樣宏觀剪切帶的形狀不同，試樣切向接觸力的集中方向決定了其宏觀剪切帶的形狀。圍壓的提高使二維和三維試樣內(nèi)切向接觸力的分布均變得更加集中在接觸法向主方向范圍內(nèi)，宏觀上表現(xiàn)為試樣剪切帶變得集中。

三維和二維數(shù)值試驗結(jié)果的差異表明，在利用簡化的二維離散元模型分析實際的三維問題時，不能一味追求數(shù)值試驗與物理試驗在數(shù)量上的一致，而應發(fā)揮二維試樣在微觀結(jié)構(gòu)定量分析方面的優(yōu)勢，從機理上解釋試樣的宏觀力學性質(zhì)，努力建立能描述試樣宏觀力學響應與其微觀結(jié)構(gòu)演化之間的定量規(guī)律。

[1] 羅勇,龔曉南,連峰.三維離散顆粒單元模擬無黏性土的工程力學性質(zhì)[J].巖土工程學報,2008,30(2):292-297.

Luo Yong,Gong Xiaonan,Lian Feng.Simulation of mechanical behaviors of granular materials by three-dimensional discrete element method based on particle flow code [J].Chinese Journal of Geotechnical Engineering,2008,30(2):292-297.

[2] 周健,池永.砂土力學性質(zhì)的細觀模擬[J].巖土力學,2003,24(6):901-906.

Zhou Jian,Chi Yong.Mesomechanical simulation of sand mechanical properties [J].Rock and Soil Mechanics,2003,24(6):901-906.

[3] 史旦達,周健,劉文白,等.砂土直剪力學性狀的非圓顆粒模擬與宏細觀機理研究[J].巖土工程學報,2010,32(10):1557-1565.

Shi Danda,Zhou Jian,Liu Wenbai,et al.Exploring macro-and micro-scale responses of sand in direct shear tests by numerical simulations using non-circular particles [J].Chinese Journal of Geotechnical Engineering,2010,32(10):1557-1565.

[4] 蔣明鏡,王富周,朱合華.單粒組密砂剪切帶的直剪試驗離散元數(shù)值分析[J].巖土力學,2010,31(1):253-258.

Jiang Mingjing,Wang Fuzhou,Zhu Hehua.Shear band formation in ideal dense sand in direct shear test by discrete element analysis [J].Rock and Soil Mechanics,2010,31(1):253-258.

[5] Li Xia,Yu Haisui.Influence of loading direction on the behavior of anisotropic granular materials [J].International Journal of Engineering Science,2009,47(11/12):1284-1296.

[6] Itasca Consulting Group,Inc.PFC2D (Particle Flow Code in 2 Dimensions),Version 3.10 [M].Minneapoli:ICG,2002.

[7] 史旦達.單調(diào)與循環(huán)加荷條件下砂土力學性質(zhì)細觀模擬[D].上海:同濟大學,2007.

Shi Danda.Micromechanical Simulations of Sand Behavior under Monotonic and Cyclic Loading [D].Shanghai:Tongji University,2007.

[8] 周國慶,周杰,陸勇,等.顆粒流程序(PFC-2D)中阻尼參數(shù)的適用性研究[J].中國礦業(yè)大學學報,2011,40(5):662-667.

Zhou Guoqing,Zhou Jie,Lu Yong,et al.Selection of damping parameters used in particle flow code (PFC-2D) [J].Journal of China University of Mining & Technology,2011,40(5):662-667.

[9] 周杰.砂土制樣過程的仿真試驗及宏-微觀分析[D].徐州:中國礦業(yè)大學,2011.

Zhou Jie.Multi-Scale Investigation on the Effect of the Specimen Preparation on Sand Behavior [D].Xuzhou:China University of Mining & Technology,2011.