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

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

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

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

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

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

1 離散單元模型

1.1 本構(gòu)模型

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

圖1 接觸模型[3]

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

(1)

(2)

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

Mr=krθ

(3)

(4)

(5)

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

(6)

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

(7)

1.2 數(shù)值模型

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

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

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

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

2 宏觀結(jié)果分析

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

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

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

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

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

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

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

(8)

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

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

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

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

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

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

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

其次，合同文本的語(yǔ)句具有如下特征：第一，多用復(fù)雜長(zhǎng)句以使語(yǔ)義表達(dá)更完整。第二，多用陳述句，因?yàn)殛愂鼍淇梢愿玫匾?guī)定和說(shuō)明商務(wù)合同雙方的責(zé)任和權(quán)利，避免歧義。第三，多使用被動(dòng)語(yǔ)態(tài)，因?yàn)楹贤瑑?nèi)容避免摻雜個(gè)人感情，被動(dòng)語(yǔ)態(tài)能更客觀公正的表達(dá)。

(9)

式(9)中：Nc為試樣內(nèi)部總接觸數(shù)：Ns為試樣內(nèi)部總顆粒數(shù)；Ns0為某顆粒周?chē)鸁o(wú)接觸的顆粒數(shù)；Ns1為某顆粒周?chē)佑|數(shù)目為1的顆粒數(shù)。

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

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

3.2 概率密度

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

(10)

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

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

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

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

3.3 法向接觸力

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

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

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

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

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

3.4 各向異性

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

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

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

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

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

4 結(jié)論

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

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

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