顆粒材料剪脹性的微觀力學分析

秦建敏 ，遲璐璐

（1. 大連大學 土木工程技術研究與開發中心，遼寧 大連 116622；2. 大連理工大學 工業裝備結構分析國家重點試驗室，遼寧 大連 116024）

1 引 言

剪脹現象是由 Reynolds[1]于 1885年首先指出的，他把顆粒材料由剪切應變引起體積變形這一特性稱為剪脹性。1936年由Casagrande[2]應用于土力學領域，并分析了土體中摩擦角對體積應變的影響。Rowe[3]采用微觀分析方法，基于規則排列結構的研究，通過對內功率的無約束最小化提出了著名的膨脹函數，得到三維空間的表達其中，σ1、σ3分別為大、小主應力，φu為顆粒間摩擦角，V為顆粒體的體積，和分別為體積和軸向應變的變化率。在這一理論中，剪脹性被假設為應力比的單一函數，推導過程中顆粒系統被假設為等效連續體，不考慮顆粒接觸處靜態和動態約束條件的限制。文獻[4]中證實了Rowe的應力剪脹理論中的原理，并提出描述顆粒材料各向異性的一個概念性方法，發現有關結構各向異性演化的一些非常重要的事實，比如大變形后土體接近一個極限各向異性狀態。文獻[5]基于光彈材料及圓鋁棒的直剪試驗結果，從顆粒接觸點數分布變化規律出發，推導出滑動面上的應力剪脹性方程。文獻[6]用離散元方法研究了土體的應力剪脹方程。文獻[7]建立了一個“砂土樣品”的顆粒力學模型，探討了剪脹性的顆粒力學影響因素。但室內試驗資料證明，膨脹性不僅依賴于應力比，還與材料的內部狀態變量有關，如孔隙比。在同樣的應力比條件下，排列密實的顆粒體發生剪脹，疏松的顆粒體發生剪縮。至今，顆粒材料的剪脹性仍未得到充分解釋，現有的剪脹理論與試驗觀測結果仍有顯著背離，關于顆粒材料剪脹性的演化過程也鮮有報道。

由于顆粒材料所表現出來的變形特性取決于內部顆粒之間、顆粒與孔隙及其他邊界之間的相互影響，是其系統內部各種材料結構要素共同作用的結果。因此，本文將從顆粒材料的內部結構入手，用孔隙胞元（void cell）描述顆粒材料的最小單元，進行雙軸剪切試驗，給出剪切過程中孔隙胞元體積變形的空間分布和演化過程，探討顆粒材料的膨脹性對加載過程中應力比與內部結構的依賴關系，解釋排列密實顆粒材料在剪切過程中先壓縮后膨脹的微觀機制。

2 孔隙胞元系統

顆粒體系的宏觀變形和力學特性，通常不是單個顆粒對外荷載的反應，而是顆粒在空間排列起來對外荷載的反應。因此，在研究顆粒材料的剪脹性之前，必須先明確顆粒體系內部的微觀結構。為了描述顆粒體系的分布特性，需要根據顆?？臻g的幾何信息、粒間接觸關系將其劃分成幾何系統。二維顆粒體的微觀拓撲結構可以由顆粒（頂點）、接觸（邊）和孔隙胞元（面）構成。顆粒系統劃分成微元存在一定的任意性，以二維情況為例，點a和b是以O為中心的顆粒的兩個接觸點，如果點a和b處顆粒表面的切線相交，則接觸點a和b稱作是一組關聯接觸。a和b以及顆粒中心O形成一個半開放的子平面 a-O-b，相鄰顆粒的子平面連接形成一個閉合空間，如顆粒O、P、Q和R的子平面形成的閉合空間為OaPcQdRbO，稱為孔隙胞元，如圖1中陰影所示。以4 000個圓形顆粒為例，對顆粒體系中所有顆粒進行循環，生成孔隙胞元系統（void cell system），如圖2(b)所示。

從孔隙胞元的構建過程看出，孔隙胞元的邊界提供了傳遞接觸力的連續路徑，因此，孔隙胞元是承受外部應力時保持穩定的最小單元[8]。將顆粒材料視為具有不同特性孔隙胞元的集合，顆粒材料的力學響應可從孔隙胞元在承受外載時的響應來研究。

圖1 孔隙胞元的構建Fig.1 Construction of void cell

圖2 4 000個隨機排列的密實顆粒體模型Fig.2 Initial arrangement of dense granular assembly of 4 000 particles

為了方便說明，以圖3所示的二維等徑顆粒體的規則排列為例進行分析，假設顆粒為剛性體，單個顆粒的變形通過接觸點變形來刻畫，不考慮顆粒的破碎。研究規則排列的等徑顆粒體，有助于揭示顆粒介質力學行為的內在機制。對于規則排列，其力學特性可通過分析代表性孔隙胞元（如圖中陰影區域所示）來展開討論。記其豎直方向和水平方向的長度分別為l1和l2，它們的比值l1/l2稱為孔隙胞元的形狀因子，定義 α= arctan(l1/l2)以描述周期性結構的幾何特征，β表示接觸點發生滑動時的角度。作用于試樣邊界的力與內部接觸力平衡，接觸力包括法向和切向分量，F11、F22以及F12、F21分別為作用在顆粒上的法向、切向力。顆粒間的摩擦力提高了顆粒體系在承受外荷載時保持穩定的能力。規定以周期性結構受壓為正，豎直方向為大主應力方向，記為σ11；水平方向為小主應力方向，記為σ22。

圖3 二維規則排列的顆粒體Fig.3 Two dimensional granular assembly with regular arrangement

當孔隙胞元受豎直方向的剪切荷載時，即σ11/ σ22＞ 1，l1/l2減小，孔隙胞元體積發生改變。僅僅從幾何上考慮，當孔隙胞元為各向同性時，其面積最大；各向異性越強，面積越小。以四邊形為例，孔隙胞元的面積為，當α=45°時，即l1=l2為正方形，面積最大。因此，受剪切荷載作用（ σ11/σ22＞1）時，在基本構型不變的前提下，對于l1/l2>1的孔隙胞元，面積增大，孔隙胞元呈膨脹性；而對l1/l2<1的孔隙胞元，接觸點的滑動導致孔隙胞元面積的減小，為壓縮的。圖3中陰影部分孔隙胞元的剪切過程如圖4所示，即從圖4(a)～(c)經過先膨脹后壓縮的過程。圖4(a)和圖 4(c)為相同孔隙比條件下孔隙胞元的兩種最基本狀態，對于α＞45°時的狀態圖4(a)，大主應力方向與孔隙胞元的方向平行，結構相對穩定，需要承受較大的主應力比才能使原有結構發生破壞，演變為圖4(b)的狀態，孔隙胞元面積增大；而在α＜45°的圖 4(c)狀態，孔隙胞元的結構不穩定，只需要較小的主應力比就能使α繼續減小，孔隙胞元的孔隙比也減小。因此，孔隙胞元的體積變化與所受應力比和形狀因子都有關系?？紫栋男螤钜蜃硬荒茉黾拥綗o窮大，當形狀因子太大時，孔隙胞元可能不能承受外部荷載。在持續的剪切荷載作用下，孔隙胞元的構型發生改變，這時在豎向有新的接觸形成，接觸力的傳遞路徑因而發生改變，微觀結構重新排列。

圖4 四邊形孔隙胞元的剪切變形Fig.4 Shear deformation of quadrangular void cell

3 離散元雙軸剪切試驗

本文采用基于孔隙胞元的離散元方法，對顆粒體試樣進行雙軸剪切試驗，詳細理論介紹見文獻[9]，限于篇幅，本文不再重復介紹。高質量的顆粒體試樣必須滿足以下3個條件：①顆粒體須滿足統計上的各向同性和均勻性；②顆粒體初始狀態時的重疊量不宜過大；③在初始狀態，必須保證顆粒體處于平衡。本文采用壓縮邊界法生成顆粒體集合：首先，將顆粒體邊界范圍按照一定的倍數放大，根據顆粒的級配曲線生成一系列的顆粒；然后，將周圍邊界壓縮，得到設定的孔隙比；最后，使顆粒在非平衡力的作用下自由運動，并在給定的邊界范圍內達到平衡狀態。取如圖2中排列致密的4 000個粒徑不等的圓形顆粒體為研究對象，對該試樣進行雙軸剪切試驗。試樣為高、寬度均為60 mm的正方形，顆粒的粒徑在（0.45～1.40）D50之間隨機分布，其中D50=1.1 mm，為顆粒體的中值粒徑，表示當小于某粒徑的土粒質量累計百分比為50%時的粒徑。不均勻系數Cu=1.69。試樣的初始狀態為緊密、隨機、各向同性排列，初始孔隙比為0.154，顆粒的初始重疊量為 1.01×10-2D50。顆粒體四周均采用周期性邊界，顆粒重度為20 kN/m3，顆粒之間的摩擦系數為0.5（顆粒間摩擦角φ=26.57°），法向接觸剛度為2.5×108N/m，切向剛度與法向剛度的比值為0.25，時間步長取 1×10-6s。在雙軸加載過程中，顆粒體試樣的高度以恒定的速度減小，采用伺服應力控制原理使水平應力σ22保持常數值。

圖5給出了排列密實的顆粒體在雙軸加載過程中應力比（ q=σ11/σ22）、孔隙比與軸向應變之間的關系曲線。從圖 5(a)可以看出，加載過程分為如下幾個階段，即：硬化階段（A- B- C）、峰值狀態（點C）、軟化階段（C- D）、臨界狀態。在加載初期，顆粒體的力學性能接近于彈性，隨著加載的進行，應力比逐漸增大，達到峰值后發生軟化，軟化現象的發生是由于密實顆粒之間有咬合力，剪切過程中咬合力被破壞，強度降低。隨著加載的進一步進行，應力比達到穩定狀態保持不變。從圖5(b)發現：在初始加載階段，顆粒體的孔隙比首先發生輕微的下降，然后隨著變形的增大，排列致密的顆粒體發生膨脹，并且剪脹現象一直伴隨在隨后的加載過程中（硬化－達到峰值－軟化），直至達到穩定態，孔隙比變化率很小，即達到臨界狀態。為了說明加載過程中顆粒材料的體積變化不僅與應力比有關，選出圖5(a)應力比-應變關系曲線中應力比相同的B和D兩點作為比較，詳細討論見第4部分。

4 剪脹性演化過程

下面詳細討論顆粒材料在細觀尺度上，即單個孔隙胞元尺度上的變形情況。對顆粒材料應變張量的微觀定義具有多樣性，大致分為兩類：①基于等效幾何圖形的應變定義，即顆粒體系可以視為等效連續體，假設相應顆粒中心的位移集中發生在接觸點處，在均勻場的假設下，建立位移梯度張量與顆粒位移以及相關幾何量的關系；②基于最優擬合的應變（best-fit strains），即基于平均位移梯度（應變）擬合出顆粒體系的位移場，尋找位移場與顆粒體系的特征位移偏差最小的平均位移梯度。在這里，特征位移可以選取顆粒中心的平動位移或接觸點處的相對平動位移。理論和計算結果表明，基于等效幾何圖形和顆粒平動位移的細觀應變定義與基于材料邊界變形的應變更為接近。本文用Kruyt等[10]提出的基于Satake的圖形系統定義的適用于二維顆粒物質的表達式。

圖5 應力比和孔隙比與軸向應變的關系曲線Fig.5 Stress ratio and void ratio vs. axial strain

二維顆粒體區域S的拓撲結構可以分成n個多邊形子結構 Si的集合，每個子結構表示由接觸顆粒的接觸支鏈圍成的孔隙胞元，如圖6所示。整個顆粒體區域S的變形率可由速度梯度L表示，L為位置x的函數，L在宏觀尺度上是均勻的，而微觀尺度（即每個孔隙胞元）上的有較大的變化。整個區域的平均速度梯度與每個孔隙胞元的平均速度梯度有如下關系：

Kuhn假設速度沿著孔隙胞元的邊是線性變化的，推導出具有m條邊的多邊形孔隙胞元Si的的表達式如下[11]：

圖6 粒子圖和對應第i個孔隙胞元的矢量圖Fig.6 Particle graph and the vectors in the corresponding i th void cell

圖7給出了初始排列為密實狀態的顆粒體在雙軸剪切過程中，孔隙胞元及其內部體積變形的演化過程。圖中ε22為軸向應變，q為應力比。從圖7中的離散元數值結果可以看出，局部孔隙胞元的變形具有明顯的非均勻性，這從圖中的顏色差異可以體現出來。以藍色填充的孔隙胞元代表其中的 Φi為正，表示該胞元內的變形為膨脹；以紅色填充的孔隙胞元中的變形率為負，該孔隙胞元內的變形為壓縮。比較圖7中各個不同加載階段發現，在雙軸加載過程中，顆粒體內部孔隙胞元的形狀發生改變，并存在由體積變化較大的孔隙胞元形成的窄帶。在初始加載狀態，孔隙胞元的形狀基本為各向同性，顆粒材料內部幾乎所有的孔隙胞元內體積應變都是壓縮的（如圖7(a)、(b)所示），這與前面圖5(b)中相同軸向應變時孔隙比的變化一致。結合圖4中四邊形孔隙胞元的3種基本狀態，在加載初期，圖4(a)狀態的孔隙胞元保持不變，而處于圖 4(c)狀態的孔隙胞元發生變化，持續加載使α變小，從而孔隙胞元面積變小，即孔隙比變小。初始狀態越密實的顆粒材料，處于圖4(c)狀態的孔隙胞元對于加載初期的剪縮現象的貢獻就越大，這就是在加載初期排列密實的顆粒材料會發生剪縮現象的原因。隨著加載的進行（如圖7(c)所示），孔隙胞元在豎直壓縮方向上變大，孔隙胞元的幾何形狀出現各向異性；接觸力較大的力鏈沿著近似平行于大主應力的方向（即豎直方向）排列，形成柱狀結構，荷載沿著相對較強的力鏈傳遞，此時孔隙胞元結構相對穩定，總體的應力-應變關系曲線為硬化狀態（對應圖5(a)）；孔隙胞元內部體積變形大部分為壓縮（對應圖5(b)中孔隙比與軸向應變關系曲線中也為壓縮），由壓縮變形大的孔隙胞元組成窄帶，窄帶中的局部孔隙胞元內壓縮變形略微大于周圍孔隙胞元的壓縮性。隨著加載的進行，顆粒體內部孔隙胞元的結構進一步改變，由于顆粒體在水平方向上持續膨脹，顆粒材料發生剪脹，孔隙胞元內部體積變形大部分為膨脹，由膨脹性大的孔隙胞元組成了窄帶（如圖7(d)所示），持續剪脹使內部很多孔隙胞元解體，顆粒體的力學響應進入軟化狀態。比較圖7(b)和圖7(d)，雖然顆粒體所處的應力水平相同，但是由于在加載過程中，顆粒體內部孔隙胞元的幾何結構發生了很大的變化，微觀尺度上體積變形也有很大的不同，其對應的宏觀尺度上力學響應也有明顯的區別。

圖7 顆粒體孔隙胞元內的體積變形Fig.7 Volume deformation in individual void cell of particulate assembly

5 結 論

以單個孔隙胞元描述顆粒材料的內部結構，分析了承受剪切荷載時其體積變化與孔隙胞元形狀、所受應力狀態的關系。對排列致密的顆粒材料進行離散元雙軸剪切試驗，得到了不同加載階段孔隙胞元系統的演化過程。

（1）微觀尺度上孔隙胞元內體積存在先壓縮后膨脹、形狀上先各向同性后沿主應力方向變長的現象，體積變形同時受孔隙胞元形狀及所受應力狀態二者的影響。

（2）在剪切過程中，孔隙胞元的體積變形存在局部化現象，由體積變形較大的孔隙胞元組成傾斜的窄帶。

（3）綜合孔隙胞元的受力分析和離散元數值結果表明，致密排列顆粒材料的剪脹性與微觀尺度上孔隙胞元的幾何結構及其內部的力鏈傳遞方式密不可分。

[1]REYNOLDS O. On the dilatancy of media composed of rigid particles in contact, with experimental illustrations[J]. Philosophical Magazine (Series 5), 1885,20: 469－481.

[2]CASAGRANDE A. Characteristics of cohesionless soils affecting the stability of slopes and earth fills[J]. Journal of the Boston Society of Civil Engineers, 1936, 23(1):257－276.

[3]ROWE P W. The stress dilatancy relation for static equilibrium of an assembly of particles in contact[J].Proceedings of the Royal Society of London (Series A):Mathematical and Physical Sciences, 1962, 269: 500－527.

[4]HOME M R. The behavior of an assembly of rotund, rigid,cohesionless particles[J]. Proceedings of the Royal Society of London (Series A): Mathematical and Physical Sciences, 1965, 286: 62－97.

[5]MATSUOKA H. A microscopic study on shear mechanism of granular materials[J]. Soils and Foundations, 1974, 14(1): 29－43.

[6]LIU S H, MATSUOKA H. Microscopic interpretation on a stress-dilatancy relationship of granular materials[J].Soils and Foundations, 2003, 43(3): 73－84.

[7]常在, 楊軍, 程曉輝. 砂土強度和剪脹性的顆粒力學分析[J]. 工程力學, 2010, (4): 95－104.CHANG Zai, YANG Jun, CHENG Xiao-hui. Granular mechanical analysis of the strength and dilatancy of sands[J]. Engineering Mechanics, 2003, 43(3): 73－84.

[8]秦建敏, 張洪武. 顆粒材料的微觀臨界狀態理論模型[J]. 巖土力學, 2010, 31(12): 3697－3703.QIN Jian-min, ZHANG Hong-wu. Microscopic theoretical model of critical state for granular materials[J].Rock and Soil Mechanics, 2010, 31(12): 3697－3703.

[9]張洪武, 秦建敏. 基于非線性接觸本構的顆粒材料離散元數值模擬[J]. 巖土工程學報, 2006, 28(11): 1964－1969.ZHANG Hong-wu, QIN Jian-min. Simulation of mechanical behaviors of granular materials by discrete element method based on mesoscale nonlinear contact law[J]. Chinese Journal of Geotechnical Engineering,2006, 28(11): 1964－1969.

[10]KRUYT N P, ROTHENBURG L. Micromechanical definition of the strain tensor for granular materials.ASME[J]. Journal of Applided Mechanics, 1996, 63(3):706－711.

[11]KUHN M R. Deformation measures for granular materials[C]//Mechanics of Deformation and Flow of Particulate Materials. New York: American Society of Civil Engineers, 1997: 91－104.