負泊松比巖體顆粒流模型及模擬分析

金愛兵，王凱，楊振偉，孟新秋

?

負泊松比巖體顆粒流模型及模擬分析

金愛兵，王凱，楊振偉，孟新秋

(北京科技大學土木與環(huán)境工程學院，北京，100083)

為了從微觀結構角度研究巖體特殊的負泊松比效應，基于顆粒流理論和顆粒流程序(PFC)，提出采用顆粒流數(shù)值模擬方法。從已知具有負泊松比效應的無機晶體材料出發(fā)，研究其內部特殊的微觀結構，進而建立符合巖體情況的顆粒流模型。利用所生成模型進行單軸壓縮模擬試驗，得出泊松比隨應力的變化曲線。基于所建模型，對巖體在低荷載水平作用下的負泊松比現(xiàn)象機理進行分析。研究結果表明：當軸壓為0~3 MPa時，巖體試樣泊松比為負值，最小泊松比為?0.059。模擬結果與室內試驗結果大致吻合，驗證了該顆粒流模型用于研究巖體負泊松比現(xiàn)象的可靠性。

巖體；負泊松比；PFC；顆粒流模型；微觀結構

西蒙?泊松定義材料負的橫向應變與縱向應變之比為泊松比[1]。在通常情況下，材料的泊松比為正值，即材料在縱向拉伸時橫向發(fā)生收縮。而負泊松比效應是指在材料縱向拉伸時，橫向發(fā)生膨脹；受壓縮時橫向發(fā)生收縮。根據(jù)經(jīng)典彈性力學，錢偉長等[2]從理論上嚴格證明在等溫條件下各向同性材料的泊松比取值范圍為?1~0.5，說明理論上存在負泊松比(negative Poisson’s ratio)物質。目前，在高分子材料領域已經(jīng)陸續(xù)研制出一些具有負泊松比效應的材料。Lakes[3]把1個長×寬×高為110 mm×38 mm×38 mm的普通聚氨酯泡沫放入1個長×寬×高為75 mm×25 mm× 25 mm的鋁制模具中，經(jīng)過一系列室內試驗處理后，所得泡孔單元呈內凹(re-entrant)結構，并測出其泊松比為?0.7。目前，有關內凹結構研究較多的是六邊形蜂窩材料，除此之外還有其他一些結構模型[4]具有負泊松比效應，如盧子興等[5]在前人研究基礎上提出了一種部分內凹、部分規(guī)則的六邊形二維多胞材料力學模型，并隨之預測了模型泊松比與剛度系數(shù)隨角度的變化關系，求解出其剪應變量，最終通過能量法計算出該力學模型的本構方程。而關于巖石類材料負泊松比的研究較少。王讓甲等[6?7]在對一組不同類別巖樣進行單軸壓縮試驗時，發(fā)現(xiàn)出現(xiàn)泊松比為負的情況；基于國內某海底隧道工程的巖體室內實驗，涂忠仁等[8]在低荷載水平作用下觀察到花崗巖泊松比為負值的現(xiàn)象，并利用Grima等[9?10]在對無機晶體材料研究時提出的“晶粒轉動”模型進行了初步解釋。但前人并未對巖體微觀結構進行深入分析。由于巖石是典型的非均質、非連續(xù)材料，采用連續(xù)介質力學理論和數(shù)值模擬方法對巖體進行分析明顯存在著不足，而離散單元法作為求解非連續(xù)介質力學問題的一個重要數(shù)值分析方法，卻可以從微觀到宏觀對巖體的破壞過程進行模擬分析[11]。為此，本文作者在前人的研究內容和室內試驗基礎上，從巖體內部的微觀結構角度出發(fā)，采用PFC程序數(shù)值模擬方法對巖體負泊松比現(xiàn)象進行研究 分析。

1 研究方法

根據(jù)材料泊松比定義，用公式表達即為

材料的彈性模量、體積模量、剪切模量、縱波速度p和橫波速度s與材料泊松比存在如下關系：

1.1 PFC顆粒單元

基于PFC程序，采用顆粒單元構建計算模型，設定顆粒間細觀接觸參數(shù)，模擬介質的基本特性隨應力環(huán)境的變化。當顆粒間黏結受外部作用影響發(fā)生破壞時，顆粒相互分離，實現(xiàn)對巖樣破壞過程的模擬。在模擬顆粒黏結破壞過程中，PFC程序提供了2種基本的顆粒黏結模型[12](如圖1所示)：接觸黏結模型和平行黏結模型，兩者合稱為黏結顆粒模型(BPM)[13]。本文所采用的是接觸黏結模型。

(a) 接觸黏結模型；(b) 平行黏結模型

由于BPM存在拉壓強度比過高及內摩擦角偏低等缺陷，為克服該問題本文采用簇單元模型(CPM，如圖2所示)來模擬巖體顆粒，該模型使用不規(guī)則的異形顆粒代替BPM介質中的圓形顆粒，通過增大顆粒之間的咬合及摩擦效應，克服BPM中顆粒形狀失真所帶來的問題[14]。

圖2 簇單元模型

1.2 負泊松比巖體PFC模型

不少科學工作者已經(jīng)研制出負泊松比材料，同時還提出了多種理論模型，并對其成因機理進行分析解釋。Grima等[9?10]在對一些無機晶體材料研究時發(fā)現(xiàn)，某些晶體內部通過彼此之間的鉸鏈連接在一起，并在微觀結構上呈方格狀排列，據(jù)此提出了負泊松比無機晶體材料的“晶粒轉動”模型，如圖3所示。

根據(jù)圖3所示的“晶粒轉動”模型，假設每個晶體單元的長度為，之間夾角為，則在每個單元內有

(a) 壓縮；(b) 拉伸

根據(jù)式(1)的定義，泊松比瞬時變化量為

由圖3和式(6)可得

聯(lián)立式(5)~(7)，求得該無機晶體材料的理論泊松比值瞬時變化量為?1。

基于圖3所示的無機晶體的“晶粒轉動”理論模型，對室內試驗中所觀察到的巖體負泊松比現(xiàn)象進行研究分析。假定部分特殊巖體具有類似圖3所示理論模型微觀結構，并將該模型轉換成可以在離散元軟件PFC中使用的顆粒流模型來研究這一現(xiàn)象，具體轉換步驟如下：1) 將“晶粒轉動”模型中的方格狀排列轉換成巖石顆粒的圓形球狀排列，即將晶體顆粒轉換成PFC中的簇單元模型(CPM)；2) 將晶體顆粒間起到連接作用的鉸鏈轉化成PFC程序中所提供的接觸黏結模型；3) 將試樣內巖石大顆粒間的空隙用較小顆粒進行填充，進而得到更符合巖體真實情況的顆粒流模型；4) 為避免模擬試驗計算結果失真，需消除PFC模型中的懸浮顆粒[15]。

對照“晶粒轉動”理論模型(見圖3)，最終得到的負泊松比巖體顆粒流理論模型如圖4所示。

2 顆粒流模擬與實例驗證

2.1 顆粒流模型

在已有室內花崗巖試樣單軸壓縮實驗的基礎上來進行顆粒流模擬試驗。按照國內有關試驗規(guī)程規(guī) 定[16]：抗壓試驗應采取高徑比為2.0~2.5的標準規(guī)則試件，為方便建模，本試驗中的巖石試樣高×直徑為120 mm×60 mm。

(a) 壓縮；(b) 拉伸；(c) 變換后壓縮；(d) 變換后拉伸

本顆粒流模型中，定義巖體試樣的墻體共4個，編號為1號、2號、3號和4號，所包圍的矩形范圍 (高×寬)為120 mm×60 mm。其中，1號和2號墻體為加載墻，3號和4號墻體在前期生成模型時起約束作用，在后期進行單軸壓縮模擬試驗時將其刪除。單軸壓縮模擬試驗PFC數(shù)值試樣通過一系列的微觀顆粒構建而成。其中：用于替換BPM的簇單元(CPM)大圓顆粒半徑=3.62 mm，小圓顆粒半徑=0.77 mm，2個小圓顆粒圓心距=4.20 mm(的含義見圖2)；填充物的顆粒半徑采用min~max的高斯分布，其中min=0.38 mm，max=0.40 mm；顆粒之間采用接觸黏結模型。對于給定尺寸的CPM模型在區(qū)域內部的指定位置均勻生成，總數(shù)為72個；而對于給定半徑范圍的充填物顆粒在空隙區(qū)域內隨機生成，顆粒總數(shù)為6 734個；PFC模型中總的顆粒數(shù)為6 806個。

為保證顆粒生成效果及效率，先生成小直徑顆粒，再將其半徑放大復原[17]。本次顆粒流模擬試驗中將與周圍顆粒接觸數(shù)量小于3的顆粒定義為懸浮顆粒，其消除步驟為：1) 搜索懸浮顆粒，將其固定并賦予懸浮顆粒速度為0 m/次；2) 將懸浮顆粒的半徑擴大20%使其與周圍顆粒充分接觸；3) 循環(huán)一定步數(shù)達到局部受力平衡狀態(tài)，若懸浮顆粒間最大接觸力小于指定容忍值時，則將此懸浮顆粒更改為普通顆粒；4) 當懸浮顆粒間最大接觸力大于指定容忍值時，將顆粒半徑縮小10%，重復步驟3)，以此消除PFC模型中所有懸浮顆粒。最終得到的巖體試樣單軸壓縮顆粒流模擬試驗PFC模型如圖5所示。

圖5 顆粒流模擬PFC 模型

2.2 模型細觀參數(shù)

本次研究中模型介質細觀參數(shù)校核自廈門東通道海底隧道中巖體泊松比出現(xiàn)負值的XASD Ⅲ級圍巖(花崗巖)巖樣[8]，其相關宏觀力學參數(shù)：單軸抗壓強度為50~80 MPa；彈性模量為4.600 GPa；泊松比為0.050。經(jīng)反復調試匹配，最終得到PFC數(shù)值模型的宏觀力學參數(shù)：單軸抗壓強度為58 MPa；彈性模量為 4.230 GPa，泊松比為0.056。此時模型所取細觀力學參數(shù)如表1所示。模擬結果與試驗結果基本一致，可將該組細觀力學參數(shù)用于實際模型計算[18?20]，利用該PFC數(shù)值模型能準確地對試樣泊松比變化情況進行研究分析。

表1 PFC模型細觀力學參數(shù)

注：接觸黏結法向和切向剛度分別為3.26×108Pa/m和1.27×108Pa/m；法向和切向強度分別為(2×107±5×106) Pa和(2×107±5×106) Pa。

2.3 模擬過程及結果

花崗巖是一種堅硬脆性巖石，本文中PFC模型試樣的單軸抗壓強度為58 MPa，而本次模擬主要目的是研究低荷載水平作用下巖體的負泊松比現(xiàn)象，所以，本次單軸壓縮模擬試驗只涉及巖體破壞過程中的裂隙閉合階段和彈性變形的初級階段。通過前期模擬計算，確定墻體上的最終法向荷載為6.5 MPa為宜。即在本次顆粒流模擬過程中，當1號和2號墻體上的法向荷載達到6.5 MPa時就停止加載。在顆粒流模擬過程中，由于PFC程序中對加載墻只能施加速率，因此，需要通過不斷調整1號和2號墻體的加載速率以控制墻體受力，從而可達到向PFC數(shù)值模型施加連續(xù)、穩(wěn)定法向荷載的目的。最終得到PFC模型試樣的應力?應變曲線、應變關系曲線和泊松比變化曲線如圖6~8所示。模擬試驗中縱向應變原始實測數(shù)據(jù)為負值，為方便作圖，將其取為相反數(shù)。

從圖6~8可知：1) 模擬試驗的應力?應變曲線與常規(guī)室內單軸壓縮試驗前期低荷載作用下的應力?應變曲線特征相同，變化規(guī)律相吻合；2)段為巖體裂隙閉合階段，段巖體已處于彈性變形階段；3) 當PFC模型試樣所受荷載較小(本試驗中小于3 MPa)時，巖體泊松比為負值；4) 當PFC模型試樣所受荷載較大(本試驗中大于3 MPa)時，巖體的泊松比逐漸變大，并恢復為正值。

圖6 模擬試驗應力?應變曲線

圖7 模擬試驗應變關系曲線

圖8 模擬試驗泊松比變化曲線

2.4 對比試驗

為驗證基于顆粒流分析方法的負泊松比巖體PFC模型的可靠性，將模擬實驗和室內單軸壓縮試驗[8]進行對比分析。對比試驗中巖石試樣為中粗粒黑云母花崗巖，試樣縱高×直徑為83 mm×50 mm，密度為 2 408 kg/m3。對比試驗中相關測試結果如表2所示，對比試驗的應變關系曲線如圖9所示。

2.5 對比結果

綜合顆粒流模擬試驗結果與室內對比試驗，結果發(fā)現(xiàn)：在低荷載水平作用下，部分特殊巖體會出現(xiàn)泊松比為負值的情況。為進一步驗證負泊松比巖體顆粒流模型的有效性，將模擬試驗與室內試驗的泊松比變化曲線進行對比分析，其結果如圖10所示。

表2 單軸壓縮測試結果

圖9 對比試驗應變關系曲線

1—模擬試驗；2—室內試驗

從圖10可以看出：1) 2條曲線走向、趨勢、變化特征基本相同；2) 當荷載大于3 MPa時2條曲線基本吻合；3) 室內試驗時，巖體出現(xiàn)負泊松比現(xiàn)象的荷載區(qū)間為0~2.3 MPa，最小泊松比為?0.023(相應軸壓為1 MPa)，模擬試驗此區(qū)間為0~3 MPa，最小泊松比為?0.059(相應軸壓為0.23 MPa)；4) 由于室內試驗應力不能連續(xù)加載，每次加載量間隔為0.5 MPa，所以，不排除荷載為0~1 MPa的區(qū)間內會存在更小的負泊松比，而這可能會使顆粒流模擬試驗結果和室內對比試驗結果擬合程度更好。

3 負泊松比現(xiàn)象形成機理分析

基于本文所提出的負泊松比巖體顆粒流模型，假定部分特殊巖體具有類似于圖4所示的理論模型微觀結構，則可對室內試驗和模擬試驗中所觀測到的負泊松比現(xiàn)象進行合理解釋。

基于上述假定，對于巖體單軸壓縮試驗，由圖4和圖8可知：階段由于應力比較小，巖體內部結構變化緩慢，巖石顆粒間發(fā)生相對旋轉、錯動，內部孔隙被壓密，巖樣向內出現(xiàn)“頸縮”的現(xiàn)象，所以，這一階段測量得到的橫向應變?yōu)樨撝担嬎愕玫降牟此杀纫矠樨撝?；階段應力較大，巖體內部的大部分孔隙壓密后又呈現(xiàn)出正常的變形規(guī)律，橫向應變開始逐漸增大，并恢復為正值，之后計算得到的泊松比也逐漸增大，并恢復為正值。

4 結論

1) 研究了具有負泊松比效應的無機晶體材料，參照其特殊的“晶粒轉動”模型，進而建立符合巖體實際情況，具有負泊松比效應的顆粒流(PFC)模型。

2) 基于顆粒流分析方法，利用負泊松比巖體PFC模型進行了數(shù)值模擬試驗，試驗中出現(xiàn)與已有室內試驗相同的結果，即部分特殊巖體會出現(xiàn)負泊松比效應。

3)將顆粒流模擬試驗與室內試驗進行對比，兩者的泊松比變化特征相吻合，低荷載時巖體泊松比為負，高荷載時恢復正值，兩者結果具有高度相似性，表明本文所提出的負泊松比巖體PFC模型可用于研究巖體的負泊松比效應。

4) 基于部分特殊巖體具有類似于本文所提出負泊松比巖體顆粒流理論模型微觀結構的假定，對巖體的負泊松比現(xiàn)象進行分析，即低壓壓縮時模型顆粒間發(fā)生相對旋轉、錯動，巖樣橫向出現(xiàn)“頸縮”現(xiàn)象，橫向應變?yōu)樨撝?，泊松比也為負值?/p>

[1] 史煒, 楊偉, 李忠明, 等. 負泊松比材料研究進展[J]. 高分子通報, 2003(6): 48–57. SHI Wei, YANG Wei, LI Zhongming, et al. Progress of materials with negative Poisson’s ratio[J]. Chinese Polymer Bulletin, 2003(6): 48–57.

[2] 錢偉長, 葉開沅. 彈性力學[M]. 北京: 科學出版社, 1956: 76–79. QIAN Weichang, YE Kaiyuan. Elastic mechanics[M]. Beijing: Science Press, 1956: 76–79.

[3] Lakes R. Foam structures with a negative Poisson’s ratio[J]. Science, 1987, 235(2): 1038–1040.

[4] Smith C W, Grima J N, Evans K E. A novel mechanism for auxetic behavior in reticulated foams: missing rib foam model[J]. Acta Materialia, 2000, 48(17): 4349–4356.

[5] 盧子興, 趙亞斌. 一種有負泊松比效應的二維多胞材料力學模型[J]. 北京航空航天大學學報, 2006, 32(5): 594–597. LU Zixing, ZHAO Yabin. Mechanical model of two-dimensional cellular materials with negative Poisson’s ratio[J]. Journal of Beijing University of Aeronautics and Astronautics, 2006, 32(5): 594–597.

[6] 王讓甲, 樊冀安, 高學之. 巖石負泊松比的探討[J]. 探礦工程, 1996(4): 17–19. WANG Rangjia, FAN Jian, GAO Xuezhi. Note on rock with negative Poisson’s ratio[J]. Exploration Engineering, 1996(4): 17–19.

[7] 王讓甲. 巖石聲波分級和巖石動彈性力學參數(shù)的分析研究[M]. 北京: 地質出版社, 1997: 96–101. WANG Rangjia. Study on classification of rock’s sound wave and parameters of rock’s dynamic-elastic mechanics[M]. Beijing: Geology Publishing House, 1997: 96–101.

[8] 涂忠仁, 楊強. 巖體負泊松比試驗研究[J]. 巖土力學, 2008, 29(10): 2833–2836. TU Zhongren, YANG Qiang. Test research on negative Poisson’s ratio of rock mass[J]. Rock and Soil Mechanics, 2008, 29(10): 2833–2836.

[9] Grima J N, Evans K E. A novel mechanism for generating auxetic behavior in reticulated foams: missing rib foam model[J]. Journal of Material Science Letter, 2000, 19(17): 1563–1565.

[10] Grima J N, Alderson A, Evans K E. An alternative explanation for the negative Poisson’s ratios in α-cristobalite[J]. Materials Science & Engineering A, 2006, 423(2): 219–224.

[11] 蔣明鏡, 陳賀, 劉芳. 巖石微觀膠結模型及離散元數(shù)值仿真方法初探[J]. 巖石力學與工程學報, 2013, 32(1): 15–23. JIANG Mingjing, CHEN He, LIU Fang. A microscopic bond model for rock and preliminary study of numerical simulation method by distinct element method[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(1): 15–23.

[12] Cho N, Martin C D, Sego D C. A clumped particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(7): 997–1010.

[13] Potyondy D O, Cundall P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329–1364.

[14] 余華中, 阮懷寧, 褚衛(wèi)江. 巖石節(jié)理剪切力學行為的顆粒流數(shù)值模擬[J]. 巖石力學與工程學報, 2013, 32(7): 1482–1490. YU Huazhong, RUAN Huaining, CHU Weijiang. Particle flow code modeling of shear behavior of rock joints[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(7): 1482–1490.

[15] Itasca Consulting Group. PFC2D(particle flow code in 2 dimensions) fish in PFC2D[M]. R Minneapolis, USA: Itasca Consulting Group, 2008: 20–22, 83–85.

[16] 趙文, 曹平, 章光. 巖石力學[M]. 長沙: 中南大學出版社, 2010: 23–25. ZHAO Wen, CAO Ping, ZHANG Guang. Rock mechanics[M]. Changsha: Central South University Press, 2010: 23–25.

[17] 周喻, 吳順川, 張曉平. 巖石節(jié)理直剪試驗顆粒流宏細觀分析[J]. 巖石力學與工程學報, 2012, 31(6): 1245–1256. ZHOU Yu, WU Shunchuan, ZHANG Xiaoping. Macro- and meso-analyses of rock joint direct shear test using particle flow theory[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(6): 1245–1256.

[18] 周喻, 吳順川, 焦建津, 等. 基于 BP 神經(jīng)網(wǎng)絡的巖土體細觀力學參數(shù)研究[J]. 巖土力學, 2011, 32(12): 3821–3826.ZHOU Yu, WU Shunchuan, JIAO Jianjin, et al. Research on mesomechanical parameters of rock and soil mass based on BP neural network[J]. Rock and Soil Mechanics, 2011, 32(12): 3821–3826.

[19] 吳順川, 周喻, 高斌. 卸載巖爆試驗及PFC3D數(shù)值模擬研究. 巖石力學與工程學報, 2010, 29(增2): 4082–4088.WU Shunchuan, ZHOU Yu, GAO Bin. Study of unloading tests of rockburst and PFC3D numerical simulation[J]. Journal of Rock Mechanics and Engineering, 2010, 29(Suppl 2): 4082－4088.

[20] 王莉, 譚卓英, 周喻, 等. 煤樣巴西劈裂試驗中聲發(fā)射特性的細觀模擬[J]. 北京科技大學學報, 2014, 36(3): 275–282. WANG Li, TAN Zhuoying, ZHOU Yu, et al. Mesoscopic simulation of the acoustic emission characteristic of coal samples in Brazilian test[J]. Journal of University of Science and Technology Beijing, 2014, 36(3): 275–282.

(編輯 羅金花)

flow code model and simulation of rock mass with negative Poisson’s ratio

JIN Aibing, WANG Kai, YANG Zhenwei, MENG Xinqiu

(School of Civil and Environment Engineering, University of Science and Technology Beijing, Beijing 100083, China)

In order to study the specific phenomenon of negative Poisson’s ratio of rock mass in the perspective of microstructure, a method of numerical simulation using particle flow theory and particle flow code (PFC) was proposed. An applicable PFC model of negative Poisson’s ratio of rock mass was constructed by studying the specific internal microstructure of inorganic crystals materials. A Poisson’s ratio curve of rock specimen was obtained in a numerical experiment of rock uniaxial compression using the PFC model. Based on the PFC model, the mechanism of negative Poisson’s ratio of rock specimen at low load level was interpreted. The results show that Poisson's ratio is negative when the axial compression is 0 to 3 MPa, and the minimum value is ?0.059. The results of numerical simulation agree with the test data. The reliability of the PFC model used to study negative Poisson’s ratio of rock mass is verified.

rock mass; negative Poisson’s ratio; particle flow code(PFC); particle flow code model; microstructure

10.11817/j.issn.1672-7207.2015.09.035

TU45

A

1672?7207(2015)09?3422?07

2014?09?24；

2014?11?01

國家自然科學基金資助項目(51074014)；中央高?；究蒲袠I(yè)務費專項資金資助項目(FRF–SD–12–002A) (Project(51074014) supported by the National Natural Science Foundation of China; Project(FRF–SD–12–002A) supported by the Fundamental Research Funds for the Central Universities)

金愛兵，博士，副教授，從事巖土工程災害治理研究；E-mail: jinaibing@ustb.edu.cn