考慮轉(zhuǎn)動阻抗對砂力學(xué)特性影響的離散元分析

王 壯， 段志波， 廖新超， 劉一鳴

(湖北工業(yè)大學(xué)土木建筑與環(huán)境學(xué)院， 武漢 430068)

離散單元法由Cundall等[1]于1979年提出，是一種從細(xì)觀出發(fā)研究物質(zhì)宏觀力學(xué)特性的方法，被廣泛應(yīng)用于砂土方面的研究。基于離散元法，程升等[2]通過雙軸壓縮模擬試驗，對南海軟黏土的細(xì)觀參數(shù)進行了標(biāo)定；史乃偉等[3]研究了臨界狀態(tài)下砂顆粒的運動狀態(tài)；蔣明鏡等[4]通過單粒組密砂的直剪模擬，分析了剪切帶形成的微觀機制。

傳統(tǒng)的離散單元模擬在生成大量顆粒來確保數(shù)值模型具有代表性的同時，為提高其模擬效率，將砂顆粒簡化為球形顆粒。顆粒之間的接觸作用常被認(rèn)為僅發(fā)生在一個點上，主要由顆粒重疊產(chǎn)生的正向力和由顆粒之間相對滑動產(chǎn)生的切向力組成，忽略了顆粒轉(zhuǎn)動的影響[5]。這導(dǎo)致了顆粒尺度細(xì)節(jié)的損失，使顆粒集合體材料的剪切強度低于實際值。1982年Oda等[6]通過對橢圓顆粒進行雙軸壓縮試驗發(fā)現(xiàn)決定顆粒材料力學(xué)機制及微觀變形機制的主要因素為顆粒之間相對轉(zhuǎn)動作用。眾多學(xué)者為考慮顆粒間轉(zhuǎn)動的影響對離散元法進行了改良，并引入轉(zhuǎn)動阻抗的概念[6-9]。Iwashita等[7]進一步拓展了離散元接觸模型，在其基礎(chǔ)上計入轉(zhuǎn)動分量。Jiang等[8]提出了能完整反映無黏性材料力學(xué)行為的抗轉(zhuǎn)動線性接觸模型，設(shè)定顆粒間的接觸為面接觸，并考慮了轉(zhuǎn)動阻抗的影響。Jiang等[9]研究表明，將接觸力矩納入離散單元接觸模型可以有效地獲得更真實的體積強度，并識別顆粒材料的變形行為。Zhao等[10]綜合比較了滾動阻力和顆粒形狀對顆粒材料內(nèi)部剪切誘發(fā)的織構(gòu)變化和各向異性的影響，發(fā)現(xiàn)滾動阻力模型可有效再現(xiàn)砂土中剪切誘導(dǎo)各向異性的主要特征。轉(zhuǎn)動阻抗的作用不可忽視且與砂土的強度以及變形有著密切聯(lián)系。

基于上述分析，針對轉(zhuǎn)動阻抗對砂土強度影響的研究已取得了一定的成果，然而在轉(zhuǎn)動阻抗微觀影響機理方面研究較少，同時二維研究導(dǎo)致顆粒細(xì)節(jié)缺失。鑒于此，基于離散元法，通過顆粒流程序PFC3D5.0，采用不同抗轉(zhuǎn)動系數(shù)和不同圍壓進行三軸模擬試驗，分析轉(zhuǎn)動阻抗對砂土強度和變形的影響，分別從平均配位數(shù)、概率密度函數(shù)、顆粒間法向接觸力和各向異性的角度分析轉(zhuǎn)動阻抗的細(xì)觀影響機制。以期為考慮轉(zhuǎn)動阻抗因素的巖土類材料離散單元建模提供理論參考。

1 離散單元模型

1.1 本構(gòu)模型

本構(gòu)模型的選取對于數(shù)值模擬的準(zhǔn)確性至關(guān)重要。本文研究采用了兩種本構(gòu)模型，砂顆粒與墻之間的接觸為線性接觸模型，砂顆粒與砂顆粒間的接觸為抗轉(zhuǎn)動線性接觸模型，其原理如圖1所示。

圖1 接觸模型[3]

線性模型的接觸點可以想象為一對彈性彈簧，線性力由具有恒定法向剛度和切向剛度的線性彈簧產(chǎn)生，其力學(xué)原理為

(1)

(2)

在線性模型的基礎(chǔ)上，抗轉(zhuǎn)動線性接觸模型增加了轉(zhuǎn)動阻抗力矩[7]，其定義為

Mr=krθ

(3)

(4)

(5)

顆粒轉(zhuǎn)動控制方程為

(6)

內(nèi)部力矩隨著接觸點處接觸件的累積相對旋轉(zhuǎn)而線性增加。這種累積的最大極限值等于當(dāng)前法向力與抗轉(zhuǎn)動系數(shù)和有效接觸半徑的乘積。

(7)

1.2 數(shù)值模型

砂真三軸數(shù)值模擬試驗主要包括以下3個階段：①制樣階段，首先生成六面無摩擦剛性墻體，形成40 mm×40 mm×80 mm的長方體，然后在長方體內(nèi)按孔隙率0.35生成顆粒，不平衡力比達到10-5時停止，以平衡制樣產(chǎn)生的初始內(nèi)力；②固結(jié)階段，通過伺服程序控制六面無摩擦剛性墻體的相對移動，施加固結(jié)所需的目標(biāo)圍壓；③剪切階段，在伺服機制保持一定圍壓的同時，通過控制上下兩面墻體的相對移動實現(xiàn)對試樣的剪切。

模擬試驗方案：砂試樣建模完成后控制變量，分別取圍壓為100、200、300 kPa，抗轉(zhuǎn)動系數(shù)取0、0.1、0.2、0.3和0.4，共計15組試樣，進行固結(jié)排水三軸剪切模擬試驗。數(shù)值模型如圖2所示，細(xì)觀參數(shù)的選取如表1所示。

圖2 砂真三軸數(shù)值模型

表1 數(shù)值模擬試驗參數(shù)

2 宏觀結(jié)果分析

圖3為圍壓σc=100 kPa條件下，不同抗轉(zhuǎn)動系數(shù)時偏應(yīng)力-軸向應(yīng)變曲線，各曲線初始階段基本呈直線且斜率一致，表明轉(zhuǎn)動阻抗對試樣初始彈性模量影響不大；隨著剪切進一步進行，偏應(yīng)力相繼達到峰值，這是由于轉(zhuǎn)動阻抗抑制了砂顆粒間相對轉(zhuǎn)動，從而破壞面延后出現(xiàn)；達到峰值后，各曲線均表現(xiàn)出軟化現(xiàn)象，且抗轉(zhuǎn)動系數(shù)越大，軟化越明顯。

圖3 偏應(yīng)力隨軸向應(yīng)變變化曲線

圖4為偏應(yīng)力峰值隨抗轉(zhuǎn)動系數(shù)變化曲線，可以看出：抗轉(zhuǎn)動系數(shù)和圍壓的增大，都會引起偏應(yīng)力峰值增大；其中偏應(yīng)力峰值增幅隨著抗轉(zhuǎn)動系數(shù)的增大而逐漸減小，表明轉(zhuǎn)動阻抗的影響逐漸趨于飽和。

圖4 偏應(yīng)力峰值隨抗轉(zhuǎn)動系數(shù)變化曲線

圖5為100 kPa圍壓時體應(yīng)變變化曲線，表現(xiàn)為先剪縮后剪漲，且轉(zhuǎn)動阻抗對剪脹現(xiàn)象有明顯增強作用。

圖5 體應(yīng)變隨軸向應(yīng)變變化曲線

內(nèi)摩擦角是衡量顆粒材料抗剪強度的重要指標(biāo)[11]，采用峰值內(nèi)摩擦角φs來衡量砂的強度，其定義為

(8)

式(8)中：σ1為大主應(yīng)力；σ3為小主應(yīng)力。

由圖6可知，試樣峰值內(nèi)摩擦角均隨著抗轉(zhuǎn)動系數(shù)的增大而增大，且增長幅度逐漸減小，與蔣明鏡等[12]研究結(jié)果一致。

圖6 峰值內(nèi)摩擦角隨抗轉(zhuǎn)動系數(shù)變化曲線

綜上所述，轉(zhuǎn)動阻抗可以有效地提高試樣的剪切強度，彌補傳統(tǒng)模擬試驗中將砂顆粒模擬為球體而導(dǎo)致抗剪強度低于實際值的缺陷。

3 細(xì)觀力學(xué)特性

3.1 力學(xué)配位數(shù)

配位數(shù)是試樣中某個顆粒與周圍顆粒的激活接觸數(shù)目，即法向接觸力大于零的接觸，是判斷顆粒材料內(nèi)部穩(wěn)定性的重要指標(biāo)。一般采用平均配位數(shù)表征顆粒集合體材料內(nèi)部顆粒的接觸狀態(tài)。本研究選用平均力學(xué)配位數(shù)[13]，其定義為

其次，合同文本的語句具有如下特征：第一，多用復(fù)雜長句以使語義表達更完整。第二，多用陳述句，因為陳述句可以更好地規(guī)定和說明商務(wù)合同雙方的責(zé)任和權(quán)利，避免歧義。第三，多使用被動語態(tài)，因為合同內(nèi)容避免摻雜個人感情，被動語態(tài)能更客觀公正的表達。

(9)

式(9)中：Nc為試樣內(nèi)部總接觸數(shù)：Ns為試樣內(nèi)部總顆粒數(shù)；Ns0為某顆粒周圍無接觸的顆粒數(shù)；Ns1為某顆粒周圍接觸數(shù)目為1的顆粒數(shù)。

圖7為不同抗轉(zhuǎn)動系數(shù)下平均力學(xué)配位數(shù)隨軸向應(yīng)變的變化曲線，曲線均呈現(xiàn)先微弱增加后下降至某一穩(wěn)定值趨勢。圖7中配位數(shù)增加是由于試樣內(nèi)部孔隙的壓縮，有效接觸變多；隨著剪切進一步進行，試樣發(fā)生剪脹，因此配位數(shù)降低；最終趨于穩(wěn)定表明砂粒間形成了穩(wěn)定的力鏈結(jié)構(gòu)，有效接觸不再改變。抗轉(zhuǎn)動系數(shù)越大，曲線出現(xiàn)拐點時的軸向應(yīng)變越大，且配位數(shù)的穩(wěn)定值越小，這是由于轉(zhuǎn)動阻抗抑制了顆粒間的相對轉(zhuǎn)動，進而延后了穩(wěn)定力鏈結(jié)構(gòu)的形成，且減小了平均配位數(shù)。圍壓對于配位數(shù)的影響符合一般規(guī)律，圍壓越大，配位數(shù)越大。

圖7 配位數(shù)隨軸向應(yīng)變變化曲線

3.2 概率密度

法向接觸力可以用概率密度函數(shù)(probability density function，PDF)進行統(tǒng)計量化。Azéma等[14]根據(jù)平均法向接觸力〈fn〉，將法向接觸力fn分為強、弱接觸力。對應(yīng)的接觸力分別隨歸一法向接觸力fn/〈fn〉呈指數(shù)分布和近乎冪函數(shù)分布狀態(tài)，其表達形式為

(10)

式(10)中：α、β為擬合參數(shù)。

圖8反映了σc=100 kPa時峰值狀態(tài)下，轉(zhuǎn)動阻抗對法向接觸力概率密度分布的影響規(guī)律，從圖8(a)可知，不同抗轉(zhuǎn)動系數(shù)條件下概率密度均隨著法向接觸力的增大而減小，當(dāng)法向接觸力較小時，轉(zhuǎn)動阻抗對其概率密度影響不明顯；而較大的法向接觸力的概率密度隨著抗轉(zhuǎn)動系數(shù)的增大而增大。

圖8 法向接觸力概率密度分布

圖9 不同抗轉(zhuǎn)動系數(shù)對應(yīng)的擬合參數(shù)

3.3 法向接觸力

圖10為峰值狀態(tài)下平均法向接觸力隨抗轉(zhuǎn)動系數(shù)的演化曲線，由圖10可知，平均法向接觸力隨抗轉(zhuǎn)動系數(shù)增加呈增長趨勢，且圍壓越大，增長趨勢越明顯。

圖10 平均法向接觸力隨抗轉(zhuǎn)動系數(shù)變化曲線

法向接觸力的空間分布可以采用三維球面直方圖[15]直觀描繪。圖11中每一根棱條表示在其方向上的歸一化局部平均法向接觸力

圖11為100 kPa圍壓時各特征狀態(tài)下的三維直方圖，直觀反映了局部平均接觸力分別在初始狀態(tài)、峰值狀態(tài)和最終狀態(tài)的空間分布情況。初始狀態(tài)時，三維直方圖均近似球體，即法向接觸力基本表現(xiàn)為各向同性；峰值狀態(tài)時，球面直方圖沿豎向拉長，表示加載方向形成強接觸力，而隨著抗轉(zhuǎn)動系數(shù)的增大，球面直方圖更加細(xì)長，其尺度范圍隨之增大，表明抗轉(zhuǎn)動阻抗增大會增強其各向異性程度，使強接觸力抵抗軸向加載的能力更強；相比于峰值狀態(tài)，最終狀態(tài)較為粗短，即其各向異性在峰后呈現(xiàn)近似軟化趨勢。

圖11 局部平均法向接觸力三維演化

3.4 各向異性

剪切導(dǎo)致砂顆粒間接觸狀態(tài)不斷變化，砂顆粒間的接觸重分布進而導(dǎo)致組構(gòu)各向異性的演化。Oda[16]結(jié)合結(jié)構(gòu)各向異性與排列密度定義了組構(gòu)張量，通過分布密度函數(shù)E(Ω)和接觸法向各向異性系數(shù)ac來表征組構(gòu)張量。

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

圖12和圖13分別為不同圍壓條件下，接觸法向各向異性系數(shù)(ac)和法向接觸力各向異性系數(shù)(an)隨軸向應(yīng)變的變化曲線。剪切前期，由于試樣發(fā)生膨脹變形，ac和an的演化規(guī)律十分相似，均表現(xiàn)為急劇增加；達到峰值后，ac趨于穩(wěn)定，這是由于試樣演化出的接觸法向各向異性程度達到上限，即試樣內(nèi)部已經(jīng)形成較為穩(wěn)定的接觸結(jié)構(gòu)，這與配位數(shù)最終趨于穩(wěn)定相互匹配；而an隨剪切進一步進行顯示出軟化趨勢，且抗轉(zhuǎn)動系數(shù)越大，軟化趨勢越明顯，與應(yīng)力應(yīng)變規(guī)律相對應(yīng)，表明法向接觸力對抗剪強度具有顯著的貢獻。對比圖12與圖13可知，整個剪切過程中，an均大于ac，這是由于法向接觸力大小的非均勻分布導(dǎo)致的。隨著抗轉(zhuǎn)動系數(shù)的增大，ac與an均隨之增大，且兩曲線拐點對應(yīng)的軸向應(yīng)變亦隨之增大，表明轉(zhuǎn)動阻抗增大了剪切誘導(dǎo)各向異性的演化程度，進而提高了砂的抗剪強度及穩(wěn)定性。

圖12 接觸法向各向異性系數(shù)隨軸向應(yīng)變變化曲線

圖13 法向接觸力各向異性系數(shù)隨軸向應(yīng)變變化曲線

4 結(jié)論

(1)轉(zhuǎn)動阻抗對砂體強度及變形規(guī)律存在顯著的影響，隨著抗轉(zhuǎn)動系數(shù)增大，峰值強度及峰值內(nèi)摩擦角均增大，試樣先剪縮后剪脹的變形特征更加顯著。宏觀結(jié)果與前人研究成果趨勢一致，驗證了本模型的可行性。

(2)轉(zhuǎn)動阻抗的作用降低了砂體的平均力學(xué)配位數(shù)，但同時增大了砂試樣內(nèi)部強接觸力的概率密度及其大小，而強接觸力主要起抵抗外部荷載的作用，因此提高了砂的抗剪強度。

(3)隨著抗轉(zhuǎn)動系數(shù)的增大，接觸法向和法向接觸力各向異性均隨之增大，法向接觸力各向異性增強主要是強接觸力各向異性的增加。從而使砂土內(nèi)部力鏈結(jié)構(gòu)更加穩(wěn)定。