不同模型維數下常規三軸試驗的離散元研究

周 杰

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

0 引 言

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

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

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

1 數值模型與參數

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

1.1 三維、二維離散元數值模型的建立

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

圖1 數值試樣內顆粒的排列結構Fig.1 Particle packing structure in 2D/3D DEM model

1.2 數值模型參數的確定

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

表1 二維和三維離散元模型參數設置

2 離散元數值試驗結果分析

采用伺服控制程序對三維和二維數值試樣進行剪切[9]，實現對實驗室試驗過程的模擬。

2.1 數值試樣的宏觀力學性質

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

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

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

由以上的分析可知，不同維數數值試樣宏觀力學行為的差異，本質上反映了其微觀結構特征的不同。

2.2 數值試樣的微觀結構演化

圖3為剪切破壞后(圍壓0.1 MPa)三維和二維數值試樣的接觸力鏈和位移場，由圖3可知，三維數值試樣的剪切帶呈“X型”，剪切帶內顆粒之間的接觸力和位移均較小；試樣上、下邊界上或左、右邊界上的接觸力和位移基本相同，以“X型”剪切帶為分界線，試樣上、下邊界附近的顆粒發生指向試樣內側的位移，而左、右邊界的顆粒發生指向試樣外側的位移。二維試樣的剪切帶呈“XX型”，試樣核心部位顆粒的位移較小；同一邊界不同部位顆粒的接觸力和位移明顯不同，僅試樣上、下邊界中心部位的顆粒發生指向試樣內側的位移，其他部位的顆粒均發生指向試樣外側的位移。

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

(1)

式中：ΔN為接觸法向ni落在Δθ范圍內的數目；N為試樣內接觸的總數目。

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

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

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

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

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

當試驗圍壓由0.1 MPa增大至1 MPa時，三維試樣的剪切帶由“X型”變為“◇型”，二維試樣的剪切帶也由“XX型”變為“X型”。圍壓增大，試樣剪切帶的范圍趨于集中(圖6)。

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

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

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

3 結 論

針對實驗室常規三軸試驗，進行了最密實狀態的、參數相同的三維和二維的離散元數值分析，通過分析數值試樣的宏觀力學響應和微觀結構特征，獲得了以下結論：

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

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

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

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

[1] 羅勇,龔曉南,連峰.三維離散顆粒單元模擬無黏性土的工程力學性質[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] 周健,池永.砂土力學性質的細觀模擬[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] 蔣明鏡,王富周,朱合華.單粒組密砂剪切帶的直剪試驗離散元數值分析[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] 史旦達.單調與循環加荷條件下砂土力學性質細觀模擬[D].上海:同濟大學,2007.

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

[8] 周國慶,周杰,陸勇,等.顆粒流程序(PFC-2D)中阻尼參數的適用性研究[J].中國礦業大學學報,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].徐州:中國礦業大學,2011.

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