考慮尺寸的三維膠結顆粒法向接觸強度估算公式

劉 芳，周向楠，蔣明鏡

(1.同濟大學 土木工程防災國家重點實驗室，上海 200092；2.同濟大學 土木工程學院 地下建筑與工程系，上海 200092；3.同濟大學 巖土及地下工程教育部重點實驗室，上海 200092)

?

考慮尺寸的三維膠結顆粒法向接觸強度估算公式

劉 芳1,2,3，周向楠2,3，蔣明鏡1,2,3

(1.同濟大學 土木工程防災國家重點實驗室，上海 200092；2.同濟大學 土木工程學院 地下建筑與工程系，上海 200092；3.同濟大學 巖土及地下工程教育部重點實驗室，上海 200092)

粒間膠結的存在導致天然結構性土體具有獨特的宏觀變形和強度特性。為了研究法向加載條件下三維膠結顆粒的接觸破壞規(guī)律及尺寸效應，采用有限元法模擬了理想球形膠結顆粒的法向接觸力學試驗，并將有限元與離散元的模擬結果進行了對比，分析了不同直徑比和厚寬比情況下膠結顆?？估瓑簭姸鹊淖兓?guī)律，得到了膠結顆粒法向接觸強度與膠結尺寸參數的定量關系。結果表明：膠結物的直徑比和厚寬比顯著影響膠結顆粒法向抗壓強度，厚寬比占主導作用，抗壓強度隨著厚寬比增加顯著降低；但膠結顆?？估瓘姸鹊某叽缧伙@著。

巖土工程；膠結材料；接觸力學；有限元模擬；尺寸效應

粒間膠結廣泛存在于天然巖土材料中(如結構性砂土、軟巖和土石混合體等)[1]，顯著影響顆粒尺度上的粒間力傳遞方式。由于外部荷載和環(huán)境改變，粒間膠結可能發(fā)生破損引起顆粒重新排列，從而導致膠結巖土材料具有獨特的宏觀變形與強度特性[2-4]。合理描述膠結顆粒微觀接觸力學行為成為發(fā)展膠結巖土材料離散元模擬方法的關鍵[5]，而描述膠結顆粒承載能力的微觀膠結破壞準則是其重要組成部分。

目前，已有學者采用微觀模型試驗研究了膠結顆粒接觸強度的影響因素，并在二維試驗條件下證實膠結尺寸對膠結接觸強度的重要影響[6-13]。為了在更加可控的試驗條件下得到膠結接觸強度與粒間膠結幾何參數的定量關系，劉芳等[14-15]采用有限元或離散元數值試驗方法研究了二維理想圓盤膠結顆粒的法向接觸行為規(guī)律，并且獲得二維膠結法向接觸強度的定量估算方法。在更加接近真實情況的三維膠結顆粒模擬方面，Z.F.SHEN等[16]采用離散元法研究了三維條件下脆性膠結顆粒在不同復雜加載條件下的力學響應和破壞準則。相對于二維顆粒接觸行為研究而言，三維顆粒接觸行為的研究仍顯不足。

因此，筆者在前期二維膠結顆粒接觸數值試驗研究基礎上[14]，采用有限元方法模擬三維球形膠結顆粒的法向接觸力學試驗，進一步研究法向壓縮和拉伸加載條件下的三維膠結顆粒接觸破壞規(guī)律及尺寸效應，以期獲得三維條件下球形膠結顆粒法向接觸抗壓和抗拉強度的定量估算方法。

1 有限元模型

筆者采用的膠結顆粒模型如圖1。土顆粒近似為理想圓球，直徑為Ds，為了便于在有限元模擬中施加荷載，取半球考慮。顆粒間附著膠結物，粒間膠結寬度為Db，粒間膠結最小厚度為Hb，粒間膠結的幾何特征通過兩個無量綱物理量來描述，即直徑比ξ=Ds/Db和膠結厚寬比η=Hb/Db，其中ξ表征膠結物與顆粒接觸面的彎曲程度，sinθ=1/ξ，角度θ的定義見圖1，取值在0°到90°之間。在本研究中，為了便于與三維離散元模擬結果[16]進行對比，Db取定值6 mm，ξ在2～5范圍變化，η在0.2～0.5范圍變化。

圖1 顆粒的模型示意Fig. 1 A schematic illustration of the particle model

筆者采用有限元程序ABAQUS對圖1所示的模型進行數值建模。由于土顆粒剛度一般遠大于粒間膠結材料剛度，為簡化分析，采用解析剛體模擬顆粒，不需要進行單元劃分，下部剛體參考點施加固端約束，上部剛體的參考點施加相應的位移和荷載條件。圖2為其中一個算例粒間膠結部分的單元網格圖(節(jié)點數1 325，單元數960)，單元類型為C3D8R。膠結物與顆粒之間的接觸設定為綁定約束，即假定接觸破壞產生在膠結物內部而不是顆粒和膠結物界面。

圖2 粒間膠結物的有限元網格Fig. 2 The finite element mesh of the inter-particle cementation

數值模擬中，通過位移控制進行法向荷載施加，即在上部剛體參考點施加法向的正/負位移模擬顆粒受拉/壓載荷方式，直至計算不收斂，從而得到膠結物的峰值拉伸荷載Rnt和峰值壓縮荷載Rnc，分別定義為膠結顆粒的抗拉和抗壓強度。

2 粒間膠結的材料模型及參數標定

粒間膠結物的材料模型采用ABAQUS程序自帶的混凝土塑性損傷模型(即CDP模型)[17]，以模擬土顆粒間的脆性膠結。模型參數根據Z.F.SHEN等的文獻中粒間膠結物材料試驗的離散元模擬曲線進行標定[16]。參數標定結果匯總見表1。

表1 粒間膠結物的材料模型參數

CDP模型中用戶自定義部分的屈服應力-非彈性應變、損傷-非彈性應變的取值由GB 50010—2011《混凝土結構設計規(guī)范》的經驗公式計算得到。圖3為粒間膠結物材料單元試驗的有限元模擬結果與離散元模擬結果的對比，有限元計算中所采用的本構模型參數見表1。單軸壓縮試驗有限元結果與離散元結果吻合較好，單軸拉伸試驗中，受本構模型所限，有限元模擬所得到的應力-應變關系呈現更為顯著的彈脆性，但峰值拉應力與離散元模擬結果在量級上比較接近。

圖3 材料參數標定Fig. 3 Material parameters calibration

3 壓縮試驗的數值模擬結果

3.1 法向荷載響應曲線

圖4和圖5為不同膠結尺寸下壓縮試驗數值模擬所得到的法向力與平均應變(即δ/Hb，δ為粒間法向位移)的關系。由圖4可見，隨著膠結直徑比ξ增大，法向力峰值略有增加；當ξ=2時，峰后法向荷載陡然下降到殘余值，呈現脆性破壞，當ξ逐漸增加時破壞模式由脆性轉為彈塑性破壞；有限元和離散元結果在峰值荷載、峰值荷載位移和破壞模式方面均比較接近。由圖5可見，隨著膠結厚寬比η的增加，法向力峰值明顯降低，有限元和離散元模擬結果基本一致。

圖4 不同膠結直徑比下法向壓縮荷載與位移關系Fig. 4 Relationship between normal compressive loading and displacement at different diameter ratios of the bond

圖5 不同厚寬比下的法向壓縮荷載與位移關系Fig. 5 Relationship between normal compressive loading and displacement at different thickness-to-width ratios of the bond

3.2 考慮尺寸的法向抗壓強度估算公式

膠結顆粒法向抗壓強度Rc定義為壓縮試驗中的峰值荷載，根據數值模擬結果，Rc與粒間膠結尺寸和顆粒大小有關，即具有尺寸效應。

若將Rc進行無量綱化，可根據F.LIU等[15]定義膠結顆粒法向抗壓強度因子如下：

(1)

根據式(1)可計算不同尺寸下有限元數值試驗所得到的膠結顆粒法向抗壓強度因子，繪制其與膠結尺寸參數的關系曲線。圖6為χc/(1+sinθ)與膠結厚寬比η的關系曲線。

圖6 抗壓強度因子與粒間膠結厚寬比關系Fig. 6 χc-η relationship obtained in FEM

由圖6可見，變量χc/(1+sinθ)與η具有較好歸一化關系，可通過冪函數進行擬合，即χc/(1+sinθ)=η-0.47。因此，膠結顆粒法向抗壓強度可用式(2)進行估算：

Rc=αηcαθc(Aσc)

(2)

式中：αηc和αθc分別為粒間膠結物的厚寬比系數和界面形狀系數，根據本次有限元數值試驗的擬合結果，αηc=η-0.47，αθc=1/(1+sinθ)。

圖7為粒間膠結物法向應變分別達到1.2εp,2εp和εmax時膠結物軸向垂直截面等效塑性應變的發(fā)展情況，其中εp為法向力達到峰值時的應變值，εmax為加載結束時的應變值。

圖7 膠結物軸對稱面上等效塑性應變發(fā)展過程Fig. 7 Development of plastic strain at the axial-symmetric cross-section of the inter-particle cementation

由圖7可見，膠結物外邊緣首先發(fā)生塑性變形，然后向內部擴展，這與試驗所觀測到的破壞首先以外緣破碎形式發(fā)生[13]是一致的。對比圖7(a)～(c)和圖7(g)～(i)，隨著η增加，塑性區(qū)更容易貫通，因此膠結顆粒法向抗壓強度有所降低；當應變達到εmax時，η=0.5情況下，膠結物塑性區(qū)已基本貫通形成剪切帶，而η=0.3情況下，受接觸面端部約束比較明顯，膠結物塑性區(qū)仍未貫通。對比圖7(a)～(c)和圖7(d)～(f)，隨著ξ增加，塑性區(qū)擴展速度減緩，表明顆粒與膠結物接觸界面越接近平面，接觸面端部約束越明顯，越有利于提高膠結顆粒法向抗壓強度。

圖8為加載過程中膠結物水平中心截面法向應力的變化情況，膠結中心點為橫坐標原點。由圖8可見，當法向荷載達到峰值之前，截面各點法向應力均有所增加，中心處應力值高于邊緣處；隨后，邊緣處產生塑性變形破壞，邊緣處應力陡然跌落，而中心處仍然處于彈性區(qū)，應力驟然上升。對于η=0.5的情況〔見圖8(c)〕，當應變達到εmax時，膠結物塑性區(qū)已基本貫通，中心處也產生塑性變形破壞，截面各點均發(fā)生應力降低；而η=0.3的情況〔見圖8(a)和(b)〕，中心處塑性區(qū)未貫通，中心處仍有較大應力值。

4 拉伸試驗的數值模擬結果

4.1 法向荷載響應曲線

圖9和圖10為不同膠結尺寸下拉伸試驗數值模擬所得到的法向拉伸荷載與平均拉伸應變(即δ/Hb)的關系。由圖9，圖10可見，區(qū)別于壓縮試驗，粒間膠結在拉伸情況下發(fā)生脆性破壞；與壓縮試驗類似，拉伸荷載峰值隨著膠結直徑比ξ增大和膠結厚寬比η減小而增加，但是尺寸影響不如壓縮試驗顯著。有限元和離散元結果在峰值荷載、峰值荷載位移和破壞模式方面均比較接近。

圖9 不同膠結直徑比下法向拉伸荷載與位移關系Fig. 9 Relationship between normal tensile loading and displacement at different diameter ratios of the bond

圖10 不同膠結厚寬比下法向拉伸荷載與位移關系Fig. 10 Relationship between normal tensile loading and displacement at different thickness-to-width ratios of the bond

4.2 法向抗拉強度及尺寸效應

與前述類似，膠結顆粒法向抗拉強度Rt可定義為拉伸試驗中的峰值拉伸荷載，相應地，可定義無量綱的膠結顆粒法向抗拉強度因子如下：

(3)

式中：σt為膠結材料單軸抗拉峰值應力，在本分析中σt=13.8 MPa(見圖3的單軸拉伸試驗模擬結果)。

圖11為χt與膠結厚寬比η的關系曲線，χt隨著η增加略有增加，趨近于1，若采用指數函數進行數據擬合，可得到最佳擬合方程為χt=e-0.036/η。因此，膠結顆粒法向抗拉強度可用以式(4)進行估算：

Rt=αηtαθt(Aσt)

(4)

式中：αηt和αθt分別為粒間膠結物的厚寬比系數和界面形狀系數，αθt近似為1.0，根據本次有限元數值試驗的擬合結果：αηt=e-0.036/η。

圖11 抗拉強度因子與粒間膠結厚寬比關系Fig. 11 χt-η relationship obtained in FEM

圖12為加載過程中膠結物水平中心截面法向拉應力分布情況，膠結中心點為橫坐標原點。由圖12可見，隨著法向拉伸荷載增加，截面各點法向應力均有所增加，中心處應力值高于邊緣處，但是中心與邊緣處應力區(qū)別不如壓縮試驗顯著，說明在拉伸荷載條件下，顆粒與膠結物接觸界面約束效應并不顯著。對比圖12(a)和圖12(c),隨著η增加，應力分布趨于均勻，表明當膠結厚度超過一定范圍內，端部約束影響不明顯。對比圖12(a)和圖12(b), 隨著ξ增加，應力分布區(qū)別不大，表明在拉伸荷載條件下，界面形狀影響可忽略不計。

圖12 膠結物水平中心截面上法向拉應力分布Fig. 12 Distribution of normal tensile stress at the horizontal central cross-section of the inter-particle cementation

5 結 論

為了研究膠結顆粒接觸破壞規(guī)律及尺寸效應的影響，筆者對壓縮和拉伸加載條件下的三維球形膠結顆粒進行了有限元數值試驗。研究了不同直徑比和厚寬比條件下膠結物抗拉壓強度的變化規(guī)律，并得到考慮膠結物尺寸的膠結顆?？箟汉涂估瓘姸鹊墓浪愎剑玫揭韵轮饕Y論：

1)膠結球形顆粒法向抗壓強度與粒間膠結物的幾何特征顯著相關，該特征可以通過兩個無量綱物理量進行描述，即膠結厚寬比η和直徑比ξ。膠結球形顆粒法向抗壓強度隨著η減小和ξ增加而增加，η的影響占主導，可以采用式(2)進行估算。

2)膠結球形顆粒法向抗拉強度的尺寸效應不是特別顯著，隨著η減小和ξ增加而略有增加，可采用式(4)進行估算，為簡化起見可忽略ξ的影響。

需要指出的是，估算公式式(2)和式(4)所涉及的擬合公式與膠結物的材料參數有關，是否適合其它膠結材料還需要進一步的驗證工作。顆粒膠結法向接觸的尺寸效應不僅限于顆粒與膠結物的直徑比和膠結物的厚寬比兩個方面，還包括膠結顆粒之間的尺寸比，筆者只考慮了簡單的等同雙球的情況，在后續(xù)的研究中需進一步的完善。

[1] CUCCOVILLO T, COOP M R. On the mechanics of structured sands [J].Géotechnique，1999, 49(6):741-760.

[2] COOP M, ATKINSON J H, HUANG J T, et al. The mechanics of cemented carbonate sands[J].Géotechnique,1993,43(1):53-67.

[3] SCHNAID F, PRIETTO P, CONSOLI N. Characterization of cemented sand in triaxial compression[J].JournalofGeotechnicalandGeoenvironmentalEngineering, 2001, 127(10): 857-868.

[4] ISMAIL M, JOER H A, RANDOLPH M F, et al. Cementation of porous materials using calcite[J].Géotechnique, 2002, 52(5): 313-324.

[5] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J].Géotechnique,1979, 29(1): 47-65.

[6] DELENNE J Y, YOUSSOUFI M S E, CHERBLANC F, et al. Mechanical behavior and failure of cohesive granular materials[J].InternationalJournalforNumerical&AnalyticalMethodsinGeomechanics, 2004, 28(15):1577-1594.

[7] JIANG M J, SUN Y G, XIAO Y. An experimental investigation on the mechanical behavior between cemented granules[J].GeotechnicalTestingJournal, 2012, 35(5):104-108.

[8] JIANG M J, SUN Y G, LI L Q, et al. Contact behavior of idealized granules bonded in two different inter-particle distances: an experimental investigation[J].MechanicsofMaterials, 2012, 55(14):1-15.

[9] JIANG M J, LIU F, ZHOU Y P. A bond failure criterion for DEM simulations of cemented geomaterials considering variable bond thickness[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 2014, 38(18): 1871-1897.

[10] JIANG M J, ZHANG N, CUI L, et al. A size-dependent bond failure criterion for cemented granules based on experimental studies[J].ComputersandGeotechnics, 2015, 69: 182-198.

[11] 蔣明鏡, 孫渝剛, 李立青. 復雜應力下兩種膠結顆粒微觀力學模型的試驗研究[J]. 巖土工程學報, 2011, 33(3): 354-360. JIANG Mingjing, SUN Yugang, LI Liqing. Experimental study on micro-mechanical model for two different bonded granules under complex stress conditions[J].ChineseJournalofGeotechnicalEngineering, 2011, 33(3): 354-360.

[12] 蔣明鏡, 肖俞, 孫渝剛, 等. 水泥膠結顆粒的微觀力學模型試驗[J]. 巖土力學, 2012, 33(5): 1293-1300. JIANG Mingjing, XIAO Yu, SUN Yugang, et al. Experimental investigation on micromechanical model of cement-bonded particles [J].RockandSoilMechanics, 2012, 33(5): 1293-1300.

[13] JIANG M J, JIN S L, SHEN Z F, et al. Preliminary experimental study on three-dimensional contact behavior of bonded granules[C]//IOPConferenceSeries:EarthandEnvironmentalScience. Warwick, UK: IOP Science, 2015.

[14] 劉芳, 李震, 蔣明鏡. 粒間膠結尺寸對圓形顆粒法向接觸響應的影響[J]. 地下空間與工程學報, 2015, 11(6):1450-1456. LIU Fang, LI Zhen, JIANG Mingjing. Effects of bond size on the normal contact behavior of two bonded disks[J].ChineseJournalofUndergroundSpaceandEngineering, 2015, 11(6): 1450-1456.

[15] LIU F, JIANG M J, SUN C, et al. On the size-dependent compressive resistance of bonded granules [J].GéotechniqueLetters, 2015, 5: 104-111.

[16] SHEN Z F, JIANG M, WAN R. Numerical study of inter-particle bond failure by 3D discrete element method [J].InternationalJournalforNumerical&AnalyticalMethodsinGeomechanics, 2015, 40(4):523-545.

[17] LEE J，FENVES G L. Plastic-damage model for cyclic loading of concrete structures [J].JournalofEngineeringMechanics, 1998, 124(8) : 892 -900．

(責任編輯：譚緒凱)

Formulas for Estimating Size-Dependent Normal Contact Resistance of 3D Bonded Spheres

LIU Fang1, 2, 3, ZHOU Xiangnan2, 3, JIANG Mingjing1, 2, 3

(1. State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, P. R. China; 2. Department of Underground Architecture & Engineering, School of Civil Engineering, Tongji University, Shanghai 200092, P. R. China; 3. Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai 200092, P. R. China)

The presence of inter-particle bonds leads to the unique macroscopic deformation and strength characteristics of natural structural soils. In order to study 3D cement particle contact failure regularity and size effect under the normal loading conditions, the contact mechanics test of ideal spherical particle cementation was simulated by finite element method. The results of the finite element simulation and those of discrete element simulation were compared. The finite element method was used to simulate normal contact response of two spheres sandwiching cementation in order to reveal the size effect on normal resistance of bonded spheres. The simulation results of finite element were compared with those of discrete element. The rule of compression/tension resistance of the bonded spheres was analyzed under different values of two dimensionless geometric parameters of the bond, i.e., diameter ratio and thickness-to-width ratio of the bond. Quantitative relationships between contact strength and cement size parameters were obtained. The results show that the compressive resistance of the bonded spheres is highly affected by the diameter ratio and the thickness-to-width ratio of the cementaneous material, and the thickness-to-width ratio plays a dominant role. The compressive resistance of the bonded spheres significantly decreases as the thickness-to-width ratio of the bond increases, while the size effect of the cement particle tensile strength is not significant.

geotechnical engineering; cementaneous material; contact mechanics; finite element simulation; size effect

10.3969/j.issn.1674-0696.2017.07.11

2016-02-05；

2016-12-25

國家自然科學基金項目(415722675；51239010；51579178)；國家科技支撐計劃課題(2013BAB06B02)；流域水循環(huán)模擬與調控國家重點實驗室自主研究課題(2015ZY05)

劉 芳(1978—)，女，廣東河源人，副教授，博士，主要從事巖土力學與工程方面的研究。E-mail：liufang@tongji.edu.cn。

TU43

A

1674-0696(2017)07-066-07