巖質(zhì)顆粒壓縮破碎過程的離散元數(shù)值模擬研究

，，，，

(1.許昌學(xué)院 土木工程學(xué)院，河南 許昌 461000；2.河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京，210098)

天然巖質(zhì)顆粒材料是土木、交通、水利等領(lǐng)域廣泛應(yīng)用的一種建筑材料，如碎石基礎(chǔ)、鐵路宕渣、土石壩等.隨著我國基礎(chǔ)建設(shè)的快速發(fā)展，工程規(guī)模越來越龐大，目前難以完全應(yīng)用膠凝材料[1]，當(dāng)顆粒在較大外力作用下會產(chǎn)生破碎現(xiàn)象，造成顆粒材料級配發(fā)生變化，影響到顆粒骨架結(jié)構(gòu)、咬合程度、粒間摩擦等一系列物理性能，改變材料的變形和剪切強(qiáng)度[2]，關(guān)系到工程的安全使用.如巴西坎潑斯諾沃斯面板堆石壩，設(shè)計高度達(dá)202 m，由于未充分考慮壩體底部顆粒破碎的影響，導(dǎo)致壩體因變形過大而發(fā)生破壞的現(xiàn)象[3].

天然巖質(zhì)材料是由多種礦物晶粒、膠結(jié)物和孔隙組成的復(fù)合材料，內(nèi)部含有不同階次隨機(jī)分布的微觀裂紋和孔隙，宏觀尺度上又被各種地質(zhì)構(gòu)造面所切割，如斷層、弱面等，造成巖質(zhì)材料存在著宏-細(xì)-微觀上的非連續(xù)性.因此，巖質(zhì)材料在宏觀上接近連續(xù)但又非完全連續(xù)，微觀上近乎破斷但又非完全破斷的復(fù)雜材質(zhì)，從而使巖質(zhì)材料具有非線性、隨機(jī)性、流變性等復(fù)雜的力學(xué)行為[4].

離散元數(shù)值模擬方法可以根據(jù)材料的屬性，能提供詳細(xì)的細(xì)觀動態(tài)參數(shù)、力鏈、位置的變化等信息，為顆粒碎裂方式、破壞機(jī)理、級配變化等問題提供有效的研究手段[5].其中將離散的基本顆粒單元，通過黏結(jié)的方式組合成可斷裂的組合單元，模擬顆粒破碎的過程，是顆粒破碎問題研究的重要方法，如膠凝顆粒模型(Bonded-particle Model，簡稱BPM).Thornton等[6-9]首先采用可破碎的BPM單元，對沖擊破碎展開了一系列研究，發(fā)現(xiàn)BPM單元密實(shí)度對顆粒的斷裂方式和破壞方向有直接影響.Mcdowell等[10,11]為了得到單顆粒破碎強(qiáng)度與數(shù)理統(tǒng)計Weibull分布規(guī)律相吻合，采用了刪除BPM單元中主要顆粒的方法.并進(jìn)行了砂樣的一維壓縮試驗(yàn)?zāi)M分析，發(fā)現(xiàn)得到的屈服應(yīng)力比試驗(yàn)中石英砂小[12]，隨后大量學(xué)者采用該類破碎組合模型展開了一系列研究.但是，沒有對顆粒組合的細(xì)觀結(jié)構(gòu)進(jìn)行詳細(xì)的探討，該破碎組合模型從細(xì)觀結(jié)構(gòu)上吻合含孔隙率的多孔巖質(zhì)材料，而大部分顆粒屬于密實(shí)結(jié)構(gòu)材料，因此模擬顆粒壓縮破碎時會釋放過多的孔隙率.Nguyen等[5]為了解決釋放多余孔隙率的問題，采用多邊形凝聚單元(BCM)模型對巴西圓盤進(jìn)行了研究，但又產(chǎn)生新的問題，如無法檢測破碎時的裂紋、單元劃分不當(dāng)造成內(nèi)部應(yīng)力集中等. 綜上所述，為了全面的研究巖質(zhì)顆粒的細(xì)觀破碎機(jī)理及對宏觀現(xiàn)象的響應(yīng)等重要問題，BPM組合單元仍需要不斷的細(xì)化和改進(jìn).

本文基于離散元方法對巖質(zhì)顆粒壓縮破碎過程進(jìn)行數(shù)值模擬分析，數(shù)值模型的建立和計算采用離散元軟件PFC2D，主要目的是展示應(yīng)用離散元方法，建立顆粒壓縮破碎試驗(yàn)仿真計算模型的方法，以及通過數(shù)值模擬得到仿真試驗(yàn)的數(shù)值結(jié)果，歸納顆粒破碎每個階段的特征.

1 離散元模擬BPM單元的基本原理

顆粒流(PFC)[13]，為了提高計算效率，將每個顆粒簡化成一個圓形顆粒.為了模擬宏觀顆粒的破碎過程，可采用凝聚顆粒單元法(BPM)模型，構(gòu)建密實(shí)狀態(tài)結(jié)構(gòu)的顆粒破碎單元，將巖質(zhì)材料的礦物晶粒等效為細(xì)觀單元，將連接方式等效為黏結(jié)模型.

PFC中有兩種接觸剛度模型：線性接觸和非線性Hertz-Mindlin接觸，同時提供了剛度本構(gòu)模型編譯接口，可以根據(jù)需要編譯新的本構(gòu)關(guān)系.PFC中有兩種黏結(jié)模型：接觸黏結(jié)和平行黏結(jié).選用點(diǎn)接觸的接觸黏結(jié)模型模擬，細(xì)觀單元間剛度模型選取線性接觸.

線性接觸模型是通過串聯(lián)的方式定義兩個相互接觸球形顆粒，接觸法向剛度由下式計算：

(1)

接觸切向剛度Ks由下式計算：

(2)

對于線性接觸模型，法向割線剛度模量kn等于法向切線剛度模量，即

(3)

式中：Fn和Un分別為兩球體顆粒法向作用力和位移重疊量；Kn由公式(1)給出.

如果法向拉力大于法向黏結(jié)強(qiáng)度，黏結(jié)發(fā)生破壞，法向和切向接觸力都為零.而當(dāng)切向作用力大于切向黏結(jié)強(qiáng)度時，黏結(jié)也會發(fā)生破壞，但接觸力不會發(fā)生改變.如果切向作用力沒超過最大摩擦力，法向力是壓力.在顆粒接觸點(diǎn)處相對位移和法向、切向接觸力分量之間的關(guān)系如圖1所示：

圖1 接觸點(diǎn)上本構(gòu)特性

2 離散元模擬的基本過程

2.1 細(xì)觀參數(shù)的標(biāo)定

對于離散單元法，細(xì)觀參數(shù)的標(biāo)定是復(fù)雜和繁瑣的過程.主要原因在于細(xì)觀參數(shù)的定義和尺度與宏觀參數(shù)存在很大的差異性，對于大部分離散材料，很難用明確的數(shù)學(xué)關(guān)系式表達(dá)宏-細(xì)觀參數(shù)之間的關(guān)系.如法向剛度kn、切向剛度ks等，影響材料屬性的參數(shù)非常多，因此需要合理的選擇標(biāo)定方法.

對于確定模型細(xì)觀物理參數(shù)值主要分為理論公式法和宏觀參數(shù)反演法.理論公式法參數(shù)標(biāo)定主要適應(yīng)于規(guī)則的連續(xù)介質(zhì)材料[14].而對于大部分材料，尤其離散顆粒材料，則需要采用反演法完成細(xì)觀參數(shù)的標(biāo)定.實(shí)際中應(yīng)先確定模擬材料的力學(xué)行為，然后與現(xiàn)場試驗(yàn)結(jié)果進(jìn)行對比，通過不斷對參數(shù)反復(fù)進(jìn)行校核，直至達(dá)到與真實(shí)材料力學(xué)效應(yīng)逼近的效果，確定合理的參數(shù)值[15].

2.2 宏觀顆粒模型的建立

對于離散元數(shù)值模擬過程中，計算時間是需要考慮的一個重要問題.針對密質(zhì)結(jié)構(gòu)材料，新的BPM單元的構(gòu)成細(xì)觀單元具有兩方面的特點(diǎn)：粒徑相同和最密實(shí)排列，該結(jié)構(gòu)避免了顆粒破碎過程中BPM單元內(nèi)部孔隙率的多余釋放，并極大提高BPM單元的生成效率.由于幾何關(guān)系的影響，采用在一個顆粒周圍依次增加顆粒數(shù)目來生成預(yù)定的BPM單元，會積累越來越大的誤差，并且生成形狀各異的宏觀顆粒非常困難，無法滿足研究需要.為了解決以上問題，提出了一種新的BPM單元生成方法.步驟如下：

圖2 宏觀顆粒生成示意圖

(1)首先，根據(jù)要生成BPM單元的尺寸自動設(shè)置足夠多細(xì)觀單元的列和行數(shù).生成緊實(shí)顆粒組合(孔隙比e=0.102 6)如圖2所示，幾何關(guān)系為：

(4)

式中：e為孔隙比；Aa為孔隙面積；A為總面積；As為固體面積；xi為橫坐標(biāo)；yi為縱坐標(biāo)；ri為細(xì)觀單元半徑.

可以通過坐標(biāo)轉(zhuǎn)換(公式(5))，生成具有一定傾角α的試樣，滿足內(nèi)部結(jié)構(gòu)變化要求，內(nèi)部結(jié)構(gòu)為傾角的宏觀顆粒.

(5)

(2)根據(jù)預(yù)定設(shè)置的列和行，根據(jù)需要自動檢索定位每一行單元的坐標(biāo)，根據(jù)需要確定每一行的單元數(shù)目，依次檢索每一行，完成模型試樣.

2.3 宏觀顆粒壓縮破碎過程的數(shù)值模擬

圖3 顆粒壓縮示意圖

單顆粒的壓縮破碎強(qiáng)度分析是研究顆粒破碎問題的重要內(nèi)容，由于顆粒屬性的千差萬別，同時又受加載方式、形狀、尺寸、應(yīng)力狀態(tài)的影響，造成顆粒破碎強(qiáng)度值離散性非常大.巖質(zhì)材料破碎主要為脆性破壞，即使同種材料同樣粒徑的情況下，測得的強(qiáng)度值有時也會超過平均強(qiáng)度值的倍，偏差遠(yuǎn)遠(yuǎn)大于彈塑性連續(xù)介質(zhì)材料[16].

為了分析天然顆粒破碎的過程，采用上節(jié)方法生成27行由517個單元組成的內(nèi)部密實(shí)，外部形狀為多棱的宏觀顆粒，加載前內(nèi)部排列具有一定夾角，如圖3所示.采用上節(jié)反演法確定單元的細(xì)觀參數(shù)，最終得到細(xì)觀參數(shù)，詳見表1.得到力-位移曲線與文獻(xiàn)[17]一致(圖4所示)，破碎最大作用力為7.36 kN，和試驗(yàn)結(jié)果7.34 kN, 誤差控制在0.3%以內(nèi).

表1 單顆粒壓縮模型參數(shù)

圖4 力-位移曲線對比圖

離散元方法可以根據(jù)需要，通過Fish語言編譯，并調(diào)用Fishcall程序，監(jiān)測黏結(jié)斷裂的位置、坐標(biāo)、數(shù)目、類型等細(xì)觀數(shù)據(jù)，分析顆粒破碎過程的變化規(guī)律.對于該類顆粒，壓縮破碎的基本過程可歸納為四個階段，如圖5、6所示：

(1)預(yù)壓階段.如圖5(a)所示，主要針對不規(guī)則形狀的顆粒，由于平整度不夠或者擺放位置的影響.在該階段與加載扳的接觸一般為點(diǎn)接觸，接觸點(diǎn)部位作用力相對較大且不均勻，一般不發(fā)生顆粒破碎或少量破碎.使擺放不平或接觸面偏小的顆粒發(fā)生轉(zhuǎn)動，主要特征是；作用力小，位移很大，而規(guī)則的顆粒或內(nèi)部強(qiáng)度低的顆粒沒有此階段.

(2)穩(wěn)定加載階段.如圖5(b)所示，顆粒與加載邊界的加載方式轉(zhuǎn)變?yōu)槊婕虞d，顆粒內(nèi)部作用力分布均勻，并相對于顆粒中心對稱.對于不規(guī)則顆粒，隨著加載作用力的增大，會發(fā)生局部破碎現(xiàn)象.并且伴隨少量的法向力鏈斷裂，切向鏈斷裂數(shù)目呈線性增加.對于規(guī)則形狀的顆粒力-位移一般成線性關(guān)系.

(3) 顆粒破碎階段.如圖5(c)所示，該階段從接觸強(qiáng)度觀察到法向力鏈斷裂數(shù)目快速增加，裂紋上下擴(kuò)展貫通，對于不同屬性顆粒破碎方式有所不同.加載作用力達(dá)到最大值，顆粒斷裂.破碎的細(xì)觀原因是由于法向黏結(jié)斷裂斷裂數(shù)目的突變.

(4)破碎完成階段.如圖5(d)所示，顆粒分裂成多個碎塊或散體顆粒，顆粒承載作用力急劇降低或變形加速變大.黏結(jié)斷裂增加到一定數(shù)目達(dá)到穩(wěn)定值，該單顆粒已完全破碎.

圖5 加載過程破碎示意圖

圖6 接觸點(diǎn)破碎數(shù)目-位移曲線

通過離散數(shù)值方法模擬顆粒破碎過程，通過反演法對細(xì)觀參數(shù)進(jìn)行標(biāo)定，能夠得到和試驗(yàn)相吻合的宏觀響應(yīng)；從細(xì)觀角度能很好的反應(yīng)顆粒破碎的過程，并發(fā)現(xiàn)顆粒的法向黏結(jié)斷裂數(shù)目的突變是顆粒承載力失效的重要特征，而切向黏結(jié)斷裂數(shù)目，在宏觀顆粒完全失去承載力之前基本保持線性增長.

3 結(jié)論

應(yīng)用離散元方法建立了模擬巖質(zhì)顆粒的數(shù)值模型，實(shí)現(xiàn)了較好的仿真模擬.數(shù)值模擬可以監(jiān)測壓縮破碎過程任意階段的作用力、黏結(jié)斷裂位置、類型、數(shù)量等物理量的變化過程，能彌補(bǔ)室內(nèi)試驗(yàn)限于技術(shù)原因無法得到某些需要的物理量的不足，提供強(qiáng)大的數(shù)據(jù)支持.

仿真模擬能得到與室內(nèi)試驗(yàn)相一致的力-位移曲線及相應(yīng)的壓縮破碎階段.

細(xì)觀顆粒的法向黏結(jié)斷裂數(shù)目的突變是宏觀顆粒承載力失效的重要特征.