楊曉娟，馬 剛，林明春，周 恒，陸 希，周 偉

(1.武漢大學水資源與水電工程國家重點實驗室，湖北 武漢 430072；2.武漢大學水工巖石力學教育部重點實驗室，湖北 武漢 430072；3.中國電建集團西北勘測設計研究院有限公司，陜西 西安 710065)

0 引 言

宏觀離散的顆粒材料作為水利工程、巖土工程、道橋工程的本質(zhì)載體，其基礎研究具有深遠的理論價值和工程意義[1]。自然界散粒體的地質(zhì)災害和水利、巖土工程的破壞失效現(xiàn)象，嚴重威脅著國民生命財產(chǎn)安全，需要從深層次審視顆粒材料的失穩(wěn)破壞問題?；谑Х€(wěn)時出現(xiàn)的宏觀變形模式可以分為兩類[2]，一類是變形局部化失穩(wěn)，應變集中在厚度很小的剪切帶、壓縮帶或膨脹帶內(nèi)，出現(xiàn)高度異質(zhì)性的局部組織；另一類是分散性失穩(wěn)，應變場連續(xù)，沒有特定的失效位置和局部組織，呈現(xiàn)整個區(qū)域均勻破壞的狀態(tài)。然而正是這類在沒有任何可見宏觀模式的情況下發(fā)生和發(fā)展的失穩(wěn)模式，被證明與邊坡崩塌、泥石流、壩體潰決、靜力觸探和沉樁過程等密切相關[3- 4]。

自Castro等[5]發(fā)現(xiàn)飽和松砂在不排水剪切條件下的靜態(tài)液化，并且在特定體變控制下，密實顆粒材料也有可能發(fā)生分散性失穩(wěn)[6]，國內(nèi)外學者圍繞分散性失穩(wěn)開展了豐富且卓有成效的試驗、理論和數(shù)值等研究。大量研究表明，初始密度[7]、排水條件[2]、加載模式[8]、圍壓水平[9]等內(nèi)外部因素顯著影響著分散性失穩(wěn)的發(fā)生，然而考慮顆粒級配影響的研究并不多見。顆粒級配被廣泛證明與顆粒材料的宏觀力學強度和細觀變形特性聯(lián)系緊密[10-11]，因此分析顆粒級配對分散性失穩(wěn)起動的影響具有十分重要的意義。

圖1 數(shù)值試樣

近年來，復雜網(wǎng)絡理論為人們提供了一種研究復雜系統(tǒng)的角度和方法，它將成千上萬個構(gòu)成系統(tǒng)的元素視為節(jié)點，將元素之間的聯(lián)系或相互作用視為邊，將真實系統(tǒng)簡化為具有拓撲性質(zhì)的網(wǎng)絡，目前這種網(wǎng)絡科學的方法已被成功應用于顆粒材料復雜力學特性的研究[12-14]，如顆粒破碎[15]、堵塞相變[16]等，為顆粒系統(tǒng)的拓撲結(jié)構(gòu)和動力學提供了新穎的視角，深刻地補充和擴展了傳統(tǒng)方法。

本文利用離散元軟件LIGGGHTS進行不同顆粒級配的顆粒體系在等體積應變加載路徑下的三軸壓縮試驗，根據(jù)Hill二階功準則開展失穩(wěn)分析，構(gòu)建顆粒接觸網(wǎng)絡并采用復雜網(wǎng)絡的理論和工具研究其演化規(guī)律，通過網(wǎng)絡特征量演化的分析加深對顆粒材料分散性失穩(wěn)的認知，尤其是探究顆粒級配與顆粒材料失穩(wěn)現(xiàn)象的相關性，揭示級配對體系宏觀力學行為的影響機制。

1 等體積應變加載路徑的三軸試驗

1.1 數(shù)值試樣

在給定立方體空間中分別隨機生成4種不同級配的5 000個初始不重疊的圓球顆粒，顆粒粒徑分別滿足1.0d50、0.8～1.2d50、0.6～1.4d50、0.4～1.6d50間的均勻分布，平均粒徑d50為5 mm。為了控制4種試樣的初始密實程度的一致性，分別在制樣過程中改變顆粒間摩擦系數(shù)為0.1、0.2、0.27、0.5，采用各向同性壓縮制樣，直至達到預定圍壓值1 MPa，最終生成試樣的初始孔隙率分別為39.45%、40.12%、39.60%、39.61%，試樣尺寸分別為82 mm×82 mm×82 mm、83 mm×83 mm×83 mm、86 mm×86 mm×86 mm、90 mm×90 mm×90 mm，如圖1所示。

試驗中采用周期性邊界，以避免邊界效應和剪切帶的形成[17]。采用Hertz-Mindlin非線性接觸模型和Mohr-Coulomb準則計算顆粒間接觸力[18]。本文DEM數(shù)值模擬所涉及的細觀參數(shù)如表1所示。

表1 DEM數(shù)值模擬細觀參數(shù)

1.2 加載方式

控制顆粒試樣3個正交方向的應變增量之間保持線性關系的加載路徑，稱為等比例應變加載路徑[8]。已有研究表明，等比例應變路徑可以有效地檢測顆粒材料的失穩(wěn)過程[19]。本文采用該路徑的1個特例——軸對稱的等體積應變加載路徑來進行DEM數(shù)值模擬，控制關系具體可表示為

(1)

式中，ε1、ε2、ε3分別為試樣在3個方向的應變。設置軸向恒定壓縮速度為1 mm/s。

2 不同級配顆粒材料的失穩(wěn)分析

(2)

在加載過程中dε1>0，因此當偏應力σ1-σ3出現(xiàn)峰值拐點時，二階功符號由正轉(zhuǎn)負，表示試樣失穩(wěn)的起始點。當偏應力σ1-σ3減小至0時，試樣完全喪失承載能力，表示最終整體分散性失穩(wěn)的發(fā)生。

不同級配的顆粒試樣在加載過程中的偏應力σ1-σ3隨軸向應變ε1的演化曲線以及p-q平面上的應力路徑如圖2所示。從圖2可知，顆粒級配較窄的兩個顆粒試樣(1.0d50、0.8d50～1.2d50)的偏應力會持續(xù)增大至試驗結(jié)束，可達到較高的偏應力水平，二階功保持為正值，p-q曲線呈現(xiàn)出典型的應變硬化路徑，試樣始終處于穩(wěn)定承載狀態(tài)；而顆粒級配較寬的兩個顆粒試樣(0.6d50～1.4d50、0.4d50～1.6d50)的偏應力存在峰值拐點，隨后減小為0，相應的二階功符號由正轉(zhuǎn)負，p-q曲線呈現(xiàn)出典型的應變軟化路徑，試樣最終發(fā)生整體性失穩(wěn)，抗剪強度消失。以上規(guī)律與相關研究的室內(nèi)試驗研究和數(shù)值模擬結(jié)果[21-22]一致，表明所開展的三軸剪切數(shù)值試驗是合理可靠的。

3 顆粒接觸網(wǎng)絡的演化

將顆粒系統(tǒng)的每一個顆粒質(zhì)心視為網(wǎng)絡的一個節(jié)點，將每一對顆粒間接觸視為連接節(jié)點的邊，由此構(gòu)建顆粒接觸網(wǎng)絡。在加載過程中，顆粒發(fā)生位置調(diào)整和重排列，顆粒間的接觸關系不斷調(diào)整，因此顆粒接觸網(wǎng)絡也處于不斷演化過程中。以表征復雜網(wǎng)絡的特征量為分析工具，量化顆粒接觸網(wǎng)絡在整個加載過程中的演化規(guī)律。

3.1 度與強度

作為復雜網(wǎng)絡中最基本的特征量，節(jié)點的度ki描述的是該節(jié)點周圍相連的邊數(shù)[23]，其含義與顆粒配位數(shù)這一傳統(tǒng)的細觀指標相同，將網(wǎng)絡中所有節(jié)點的度取平均值，作為網(wǎng)絡的平均度〈k〉，以反映顆粒系統(tǒng)的接觸密度。接觸網(wǎng)絡的度與度分布如圖3所示。

圖3 接觸網(wǎng)絡的度與度分布

從圖3a可知，對于級配較窄的兩組顆粒試樣，其網(wǎng)絡平均度經(jīng)過初期調(diào)整后保持基本穩(wěn)定的水平，而對于級配較寬的兩組顆粒試樣，其網(wǎng)絡平均度在臨界失穩(wěn)時驟降至零附近，說明試樣觸發(fā)整體性失穩(wěn)時，內(nèi)部有效接觸結(jié)構(gòu)完全崩塌，幾乎所有顆粒呈現(xiàn)“漂浮”狀態(tài)，無法形成傳遞外荷載的有效路徑，這就是試樣承載能力消失的根本原因。因此可以根據(jù)網(wǎng)絡平均度區(qū)分顆粒系統(tǒng)是否發(fā)生了整體性失穩(wěn)，并且可以很好地捕捉失穩(wěn)的時刻，表明借助復雜網(wǎng)絡的工具開展顆粒材料失穩(wěn)現(xiàn)象的研究是合理有效的。此外，雖然4組顆粒試樣的孔隙率相同，但由于顆粒級配的差異，它們的接觸密度有較大不同，隨著級配加寬，內(nèi)部接觸結(jié)構(gòu)越松散，越容易觸發(fā)失穩(wěn)。

為了衡量節(jié)點度的分布情況，圖3b給出了顆粒接觸網(wǎng)絡中度的標準差在加載過程中的演化。從圖3b可知，對于發(fā)生失穩(wěn)的試樣，其值最終驟降為0；對于穩(wěn)定承載的試樣，其值保持穩(wěn)定，這與網(wǎng)絡平均度的演化規(guī)律一致。此外，接觸分布情況與顆粒級配存在明顯相關性，級配越寬，度的標準差越大，則接觸分布的不均勻性越強；反之對于單顆粒的試樣，接觸分布最為均勻。這可能是由于粒徑相似的顆粒更傾向于擁有相似的接觸度，而粒徑差異越大，接觸度差異越大，接觸分布越不均勻。

為了更全面地捕獲顆粒材料的豐富行為，將接觸力的信息引入到接觸網(wǎng)絡中[24]?；仡櫳鲜龉?jié)點度的計算，實質(zhì)上是將每條邊的權重視為1，合計與該節(jié)點相連的所有邊的權重。現(xiàn)為了合并接觸力信息，將對應的接觸力絕對值作為每條邊的權重，類似的，合計與該節(jié)點相連的所有邊的權重，定義為節(jié)點的強度si，同理取所有節(jié)點的均值作為網(wǎng)絡的平均強度〈s〉，以衡量顆粒承載接觸力的平均水平。

4組顆粒試樣的網(wǎng)絡強度的平均值、標準差在加載過程中的演化曲線如圖4所示。從圖4可知，在初始時刻，由于圍壓相等，顆粒系統(tǒng)的承載狀態(tài)相同，因此各顆粒試樣的網(wǎng)絡平均強度相當，但強度分布情況仍有差異，級配越寬，強度的標準差越大，則強度分布的不均勻性越強。這可能是由于粒徑相似的顆粒更傾向于承載相似大小的接觸力，而粒徑差異越大，承載情況差異越大，強度分布越不均勻。隨著加載的進行，平均強度的演化規(guī)律與偏應力相似，對于兩個級配較窄的、未失穩(wěn)的顆粒試樣，其平均強度可持續(xù)增長至較高水平，導致其強度范圍加寬，強度分布也愈發(fā)不均勻；對于兩個級配較寬的、最終發(fā)生失穩(wěn)的顆粒試樣，其平均強度和標準差都在臨界失穩(wěn)瞬間驟降為0。

圖4 接觸網(wǎng)絡的強度及其分布

3.2 局部效率

網(wǎng)絡是由節(jié)點和邊的交替序列構(gòu)成的，因此對與路徑距離相關量的描述也尤為重要。一般將節(jié)點i和j之間可遍歷路徑的最少邊數(shù)定義為節(jié)點間距離dij，但往往認為，若兩個節(jié)點之間不存在可達路徑，則它們之間的距離為無窮大，因此在不完全連通(具有不止一個連通團簇)的顆粒接觸網(wǎng)絡中，距離特征量達無窮大，為了避免距離的這一發(fā)散問題，提出“效率”的概念[25]，定義節(jié)點間效率為距離的倒數(shù)，直觀地理解為：距離越遠，則傳輸信息的效率越低，反之亦然。將所有節(jié)點對間的效率取平均值，得到網(wǎng)絡全局效率。網(wǎng)絡全局效率與平均度的物理意義類似，可以反映顆粒接觸結(jié)構(gòu)的整體連通程度，不再展開分析?，F(xiàn)為了研究局部結(jié)構(gòu)的連通情況，引入局部效率[25]，與全局效率不同的是，它衡量的是以節(jié)點為中心的局部團簇的連通情況。統(tǒng)計某一節(jié)點周圍的鄰居節(jié)點集合，將該集合中所有節(jié)點對間的效率取平均值，得到該節(jié)點的局部效率Eloc，i，本質(zhì)上是由該節(jié)點的鄰居節(jié)點所構(gòu)成的子網(wǎng)絡的平均全局效率。將節(jié)點的局部效率取平均值，得到網(wǎng)絡局部效率Eloc，一般認為網(wǎng)絡局部效率越大，相鄰顆粒所組成的局部團簇排列越緊密，網(wǎng)絡的集團化現(xiàn)象越顯著。

(3)

(4)

式中，Vi為節(jié)點i的鄰居節(jié)點集合；N為網(wǎng)絡中的節(jié)點總數(shù)；N=|V|；ki為節(jié)點i的度值。

接觸網(wǎng)絡的局部效率的演化如圖5所示。從圖5可以看出，隨著級配的加寬，網(wǎng)絡局部效率減小，局部結(jié)構(gòu)的連通密度降低，越容易發(fā)生全局化的失穩(wěn)，這與整體結(jié)構(gòu)連通性的規(guī)律一致。此外，對于兩個級配較窄、未失穩(wěn)的顆粒試樣，其網(wǎng)絡局部效率均有一定的上升，而對于兩個級配較寬、最終發(fā)生失穩(wěn)的顆粒試樣，其網(wǎng)絡局部效率在整體失穩(wěn)前出現(xiàn)持續(xù)明顯的下降，這表明在整體接觸結(jié)構(gòu)瓦解之前，局部結(jié)構(gòu)的聚集程度已經(jīng)出現(xiàn)顯著的降低，換言之，局部結(jié)構(gòu)連通性的持續(xù)降低，在一定程度上預示了最終分散性失穩(wěn)的發(fā)生。

圖5 接觸網(wǎng)絡的局部效率的演化

圖6 接近中心性的分布(初始狀態(tài))

3.3 節(jié)點中心性

由于顆粒級配對顆粒系統(tǒng)的拓撲結(jié)構(gòu)和力學行為影響顯著，為了研究不同粒徑的顆粒在系統(tǒng)中的角色定位，需要關注一個節(jié)點在網(wǎng)絡中的重要程度與其顆粒粒徑的關系。其中一個典型的特征量是節(jié)點的接近中心性CC(closeness centrality)[23]，它是所有可達節(jié)點n-1到目標節(jié)點的距離均值的倒數(shù)，即

(5)

式中，Ui為目標節(jié)點i所屬連通子圖的節(jié)點集合；n為該連通子圖的節(jié)點數(shù)，n=|Ui|。為了避免由于可達節(jié)點過少而導致接近中心性偏大的問題，提高在完全連通圖和不完全連通圖的適用性和可比性，通常運用wf修正，根據(jù)目標節(jié)點所屬連通子圖的規(guī)模進行歸一化，即

(6)

從上式可以看出，節(jié)點的接近中心性衡量了目標節(jié)點與其他可達節(jié)點之間的平均接近程度，值越高暗示著節(jié)點的中心性越高。

4組顆粒試樣在初始狀態(tài)時接近中心性的分布情況如圖6所示。從圖6可以看出，級配范圍越寬，接近中心性的值越小，進一步印證了接觸結(jié)構(gòu)的松散，同時看到，接近中心性分布的不均勻性增強，漂浮無接觸(數(shù)值為0)的顆粒增多，并且在同一個顆粒試樣中，粒徑較小的顆粒更傾向于較小的接近中心性值，這表明從與其他顆粒的接近程度的角度來看，小粒徑顆粒發(fā)揮的重要性不及大粒徑顆粒。

在許多現(xiàn)實世界的網(wǎng)絡系統(tǒng)中，節(jié)點對之間的路徑存在多條，信息往往沿著其中最短的，即最有效的，可用的路徑傳播。因此從最短路徑角度來量化節(jié)點的重要程度是極具現(xiàn)實意義的，節(jié)點介數(shù)BC(betweenness centrality)[23]就是這樣一個代表性的特征量，定義為網(wǎng)絡所有節(jié)點對的最短路徑中經(jīng)過目標節(jié)點的數(shù)量占比，為了便于結(jié)果在不同規(guī)模的網(wǎng)絡中進行比較，再根據(jù)除目標節(jié)點外的節(jié)點對總數(shù)(N-1)(N-2)/2進行歸一化，其計算表達式為

(7)

式中，δjk為節(jié)點j、k之間的最短路徑數(shù)量；δjk(i)為節(jié)點j、k之間的最短路徑且經(jīng)過節(jié)點i的數(shù)量。節(jié)點介數(shù)衡量了目標節(jié)點對最短路徑的貢獻程度，對節(jié)點對之間信息交流的重要程度，值越高暗示著節(jié)點的中心性越高。在顆粒材料的物理背景下，節(jié)點介數(shù)越高的顆粒意味著其對有效傳力路徑的參與度越高，對邊界之間傳遞外荷載的影響程度越大。

圖7 節(jié)點介數(shù)的分布(初始狀態(tài))

4組顆粒試樣在初始狀態(tài)時節(jié)點介數(shù)的分布情況如圖7所示。從圖7可知，隨著級配的加寬，節(jié)點介數(shù)最大值總體呈現(xiàn)增大的趨勢，不均勻性增強，并且在同一個顆粒試樣中，小粒徑顆粒的介數(shù)值普遍比大粒徑顆粒的值小，即節(jié)點介數(shù)與粒徑成反比，這一相關性規(guī)律比接近中心性與粒徑的關系更顯著。在寬級配顆粒試樣中，正是小顆粒對有效傳力路徑的參與度較低，因此凸顯了大顆粒在傳力路徑中的重要性，導致節(jié)點介數(shù)的最值增大，這表明從對最短路徑的貢獻程度的角度來看，大顆粒構(gòu)成了傳遞外荷載的有效路徑，而小顆粒起到了填充作用。

節(jié)點的接近中心性和介數(shù)從兩個不同角度衡量了粒徑不同的顆粒的重要性，均說明了隨著級配的加寬，主要起填充輔助作用的小顆粒數(shù)量增多，導致大顆粒需要承擔更大的負載，不利于形成穩(wěn)定承載的整體結(jié)構(gòu)，因此更容易發(fā)生失穩(wěn)。

4 結(jié) 論

本文對不同級配顆粒材料進行了等體積應變加載路徑下的三軸剪切數(shù)值模擬，根據(jù)Hill二階功準則判定材料是否發(fā)生分散性失穩(wěn)，采用復雜網(wǎng)絡理論和工具研究顆粒接觸網(wǎng)絡的演化規(guī)律，探究顆粒級配對顆粒材料力學行為的影響機制。結(jié)論如下：

(1)對于不同級配的顆粒材料，級配較寬的顆粒體系更可能且更急劇出現(xiàn)分散性失穩(wěn)。當體系觸發(fā)整體性失穩(wěn)時，內(nèi)部有效接觸結(jié)構(gòu)完全崩塌，幾乎所有顆粒呈現(xiàn)“漂浮”狀態(tài)，無法形成傳遞外荷載的有效路徑，這就是顆粒體系承載能力消失的根本原因。

(2)級配對顆粒系統(tǒng)的拓撲結(jié)構(gòu)和力學行為影響顯著。隨著級配加寬，顆粒接觸結(jié)構(gòu)的連通密度降低，激勵的接觸力強度減小，但接觸分布和強度分布的不均勻性增強，這可能是由于粒徑相似的顆粒更傾向于擁有相似的接觸度和接觸力。

(3)局部結(jié)構(gòu)連通性的持續(xù)降低，在一定程度上可以預示最終分散性失穩(wěn)的發(fā)生。對于發(fā)生失穩(wěn)的顆粒體系來說，在整體接觸結(jié)構(gòu)瓦解之前，局部結(jié)構(gòu)的聚集程度已經(jīng)出現(xiàn)顯著降低；而對于穩(wěn)定承載的顆粒體系來說，局部結(jié)構(gòu)的聚集程度有所增強。

(4)不同粒徑的顆粒在顆粒材料中的角色定位不同，這影響著顆粒材料的力學行為。從與其他顆粒的接近程度的角度來看，小粒徑顆粒發(fā)揮的重要性不及大粒徑顆粒；從對最短路徑的貢獻程度的角度來看，大顆粒構(gòu)成了傳遞外荷載的有效路徑，而小顆粒起到了填充作用。隨著級配范圍的加寬，主要起填充輔助作用的小顆粒數(shù)量增多，導致大顆粒需要承擔更大的負載，不利于形成穩(wěn)定承載的整體結(jié)構(gòu)，因此觸發(fā)失穩(wěn)的風險更大。