梁紹敏 馮云田 趙婷婷 王志華

* (太原理工大學(xué)機械與運載工程學(xué)院,太原 030024)

? (斯旺西大學(xué)辛克維奇計算工程中心,英國斯旺西 SA1 8EP)

引言

顆粒材料在自然界和工程領(lǐng)域普遍存在,如谷物、砂石和土壤等.在水利、港口和交通等工程領(lǐng)域中,堆石、砂石和砂土等顆粒材料得到了廣泛應(yīng)用,在堤壩、隧道和地基等工程中,大載荷作用下顆粒材料容易發(fā)生破碎現(xiàn)象,且顆粒材料的破碎行為會給工程和建筑建設(shè)帶來極大的影響,如三峽水利樞紐工程中的填筑材料,其破碎率可能會達20%.因此,顆粒材料的破碎行為已經(jīng)引起人們的關(guān)注,并且逐漸成為顆粒材料力學(xué)性質(zhì)研究的新課題[1].此外,顆粒的破碎行為會引起顆粒系統(tǒng)級配的改變,從而影響其物理力學(xué)性質(zhì),大顆粒破碎后產(chǎn)生的小顆粒填充到顆粒系統(tǒng)孔隙中也會導(dǎo)致系統(tǒng)孔隙減小從而引起變形[2].同時,分析細觀尺度下顆粒的破碎現(xiàn)象與宏觀尺度下材料的非線性行為之間的關(guān)系,對顆粒破碎的宏、細觀機理揭示具有重要意義[3].由此看來分析顆粒破碎過程既有實際的工程意義,又有理論的研究價值.

顆粒破碎是指顆粒材料在外載荷作用下發(fā)生開裂、崩解或磨耗,從而導(dǎo)致級配改變的現(xiàn)象,且受加載方式、外載荷大小以及顆粒形狀、尺寸、自身強度等影響[4].顆粒材料的破碎行為會導(dǎo)致顆粒形狀、尺寸、強度、配位數(shù)、級配以及密實度等發(fā)生變化,甚至引起顆粒系統(tǒng)運動狀態(tài)的變化.顆粒的破碎過程存在以下3 個特征狀態(tài): 起始狀態(tài)、臨界狀態(tài)和終止狀態(tài).把顆粒開始涉及破碎的初始狀態(tài)定義為起始狀態(tài),一般破碎的起始狀態(tài)可人為地合理指定,所以具有相對性;顆粒材料的應(yīng)力和應(yīng)變不再變化,而偏應(yīng)變繼續(xù)發(fā)生變化的狀態(tài)稱為臨界狀態(tài),臨界狀態(tài)是客觀的;當顆粒材料內(nèi)部不再有破碎現(xiàn)象發(fā)生的極端狀態(tài)稱為終止狀態(tài),所以終止狀態(tài)具有唯一性.顆粒破碎具有不可逆性、累加性、應(yīng)力路徑相關(guān)性和卸載不產(chǎn)生顆粒破碎的特點[4].

目前,對顆粒破碎的研究主要采用試驗的方法,研究內(nèi)容主要有破碎現(xiàn)象的描述、評價指標的建立、破碎影響因素分析以及破碎對顆粒力學(xué)特性的影響等內(nèi)容.國內(nèi)外已經(jīng)開展了大量針對顆粒破碎的試驗研究[5-6],并取得了一定的成果,如建立了顆粒破碎的評價指標,討論了影響顆粒破碎的因素,并在試驗結(jié)果的基礎(chǔ)上建立了顆粒彈塑性本構(gòu)模型[7].典型的工作有通過改變級配、相對密度、圍壓等條件對顆粒材料進行大型三軸剪切試驗,分析其強度、變形及剪脹特性[8].

顆粒破碎過程是一個長期的、漸進的過程,而實驗室試驗僅僅可以模擬導(dǎo)致顆粒破碎的自然條件或漫長加載歷史中的某一階段[9],無法完全還原自然界中真實顆粒的破碎過程.此外,由于試驗過程中無法滿足試驗材料和作用載荷的同一性,因此難以保證多次試驗初始條件的一致性,導(dǎo)致結(jié)果的重復(fù)性較低,這樣獲得的研究結(jié)果對研究真實破碎過程意義不大.即使可以保證起始應(yīng)力狀態(tài)的一致性,不同加載歷史或路徑也將導(dǎo)致顆粒破碎后的大小以及破碎過程的能量變化不同.此外,試驗研究還有成本高、耗時費力、難以揭示破碎的細觀機理并存在一定的危險性等缺點.數(shù)值方法的快速發(fā)展和計算機技術(shù)的飛速提高,使得通過數(shù)值方法研究顆粒破碎成為可能.與實驗室試驗相比,數(shù)值模擬更容易獲得物理參數(shù)或現(xiàn)象,如裂紋的萌生、擴展等,因此,數(shù)值模擬方法得到發(fā)展和推廣.當然,實驗室試驗對驗證理論和數(shù)值方法的有效性和準確性方面發(fā)揮了很重要的作用.

本文首先整理了目前研究顆粒破碎使用較多的數(shù)值分析方法,包括黏結(jié)?破碎法、碎片替代法、比例邊界有限元法、組合有限元?離散元法和近場動力學(xué)法,包括方法的提出、發(fā)展、實現(xiàn)過程以及方法的優(yōu)勢和存在的問題.然后從目前的主要研究成果著手,詳細講述了破碎過程的細觀機理研究和顆粒材料破碎在工程中的應(yīng)用.最后對目前有關(guān)顆粒破碎的研究進行總結(jié),并對未來的發(fā)展趨勢進行展望.

1 基于離散元法的顆粒破碎數(shù)值算法

采用數(shù)值方法研究顆粒破碎過程時,考慮顆粒的離散特性,基于離散元方法理論則是最為合理和有效的方法.這一節(jié)將主要介紹基于離散元方法原理的黏結(jié)?破碎法和碎片替代法.

1.1 黏結(jié)?破碎法

1.1.1 黏結(jié)?破碎法的提出與發(fā)展

2000年Robertson[10]首次提出顆粒的黏結(jié)模型(bonded particle model,BPM),詳細介紹了如何將單顆粒黏結(jié)成組合顆粒.2004 年P(guān)otyondy 等[11]對黏結(jié)顆粒通過二維和三維離散元軟件PFC 對巖石的破碎過程進行了仿真,隨后該方法被應(yīng)用于脆性顆粒材料破碎過程的研究中.黏結(jié)?破碎模型將一定數(shù)量的基礎(chǔ)顆粒通過黏結(jié)鍵形成顆粒簇,當基礎(chǔ)顆粒間的黏結(jié)作用失效時,顆粒發(fā)生破碎行為.

黏結(jié)?破碎法在一定程度上推進了顆粒破碎的研究,并在巖土工程中解決了一些問題.如基于黏結(jié)?破碎法,構(gòu)建土壤的顆粒簇,并通過壓縮試驗研究了顆粒材料破碎對土壤臨界狀態(tài)的影響[12];構(gòu)造砂土的顆粒簇,對砂土顆粒在高應(yīng)力下的破碎效應(yīng)進行模擬,通過改變顆粒的細觀參數(shù)、外載荷的加載速率、顆粒簇的數(shù)量以及黏結(jié)顆粒的平均粒徑等分析了破碎現(xiàn)象的影響因素[13];構(gòu)造堆石料的顆粒簇,通過雙軸剪切數(shù)值試驗,再次肯定了該方法模擬真實顆粒材料破碎特性的有效性[14];構(gòu)造不規(guī)則形狀的堆石顆粒簇,用以模擬堆石顆粒在破碎過程中的尺寸效應(yīng)[15];采用球形單元黏結(jié)成不規(guī)則的碎石料顆粒簇,改變應(yīng)力加載路徑分析顆粒材料破碎對碎石料應(yīng)力–應(yīng)變關(guān)系及體應(yīng)變特性的影響,再次證明了離散元方法能夠較為準確地獲得碎石顆粒的力學(xué)特征[16].Lu 等[17-18]采用不可破碎的“超級顆?！睒?gòu)造真實形態(tài)顆粒,將一定數(shù)量的顆粒與“超級顆?！蓖ㄟ^平行黏結(jié)鍵連接,通過判斷平行黏結(jié)失效與否來確定小顆粒與“超級顆?！钡暮想x,從而實現(xiàn)顆粒破碎過程模擬,并將其用于對道砟顆粒棱角斷裂過程的分析,破碎過程如圖1 所示.以上研究再現(xiàn)了顆粒材料破碎現(xiàn)象,對顆粒破碎過程的認識有一定的推動作用.

圖1 黏結(jié)顆粒構(gòu)造示意圖[17-18]Fig.1 Schematic diagram of bonded particles[17-18]

以上研究主要是采用球體作為基礎(chǔ)單元通過黏結(jié)鍵構(gòu)造顆粒簇,為了更接近真實顆粒形態(tài),以正四面體為核通過晶胞繁衍的方法生成堆石顆粒,對三軸剪切作用下堆石料的破碎過程進行研究[19].針對復(fù)雜顆粒形狀,湘潭大學(xué)龍志林團隊[20]采用凸多面體塊作為基本顆粒,通過可逆函數(shù)法對顆粒破碎強度變異性進行建模,如圖2 所示.此外,采用傳統(tǒng)的修正“巴西”準則作為破碎準則,再根據(jù)顆粒的接觸點和質(zhì)心可以確定顆粒內(nèi)部最終斷裂,一旦目標顆粒滿足斷裂準則,它就會被一系列與最終斷裂相一致的虛擬切割面切割成多個碎片.這樣可以避免局部應(yīng)力的產(chǎn)生及質(zhì)量和體積的不守恒.這種方法不必預(yù)定義碎片的模式[20-23].該研究將黏結(jié)基礎(chǔ)單元由球體發(fā)展到多面體,這使黏結(jié)顆粒的幾何形狀與真實顆粒更加接近.

圖2 破碎過程的實現(xiàn)[20]Fig.2 Realization of breakage process[20]

目前國內(nèi)外已經(jīng)發(fā)展了多種離散元模擬軟件,其中采用黏結(jié)?破碎法實現(xiàn)顆粒的破碎分析的主要有PFC、大連理工大學(xué)季順迎課題組研發(fā)的SDEM、南京大學(xué)劉春課題組的MatDEM 等.在處理多顆粒破碎問題時,EDEM 也是基于黏結(jié)?破碎模型實現(xiàn)的.

1.1.2 黏結(jié)?破碎法對顆粒破碎細觀力學(xué)機理的研究

Griffith 強度理論認為固體材料中包含眾多細小裂紋,在外力作用下這些裂紋的尖端會出現(xiàn)應(yīng)力集中現(xiàn)象,且這個應(yīng)力始終為拉應(yīng)力,稱之為誘導(dǎo)應(yīng)力,當這個應(yīng)力大于材料的抗拉強度時,尖端的裂紋將發(fā)生擴展并導(dǎo)致材料最終的破壞[3],這就是材料破碎機理研究的起源.借鑒材料的破壞機理研究方法,考慮顆粒材料的破碎機理,需要討論破碎產(chǎn)生的原因,包括引起破碎的外因、內(nèi)因,破碎發(fā)生的初始狀態(tài)、發(fā)展過程以及最終的狀態(tài).同時,要對破碎進行描述和衡量,還需要提出合理的破碎指標[24].顆粒的破碎機理直接決定破碎準則的建立,對破碎過程的模擬具有決定性作用[25].

通過關(guān)注顆粒內(nèi)部接觸力的變化情況以及破碎發(fā)生的位置,從而獲得顆粒破碎的演化規(guī)律,進而揭示顆粒破碎的細觀力學(xué)機理.Cil 等[26]采用數(shù)值方法記錄了黏結(jié)顆粒內(nèi)裂紋的產(chǎn)生和擴展過程,并與X-ray 掃描結(jié)果進行對比,發(fā)現(xiàn)模擬結(jié)果與掃描結(jié)果基本一致,如圖3 所示,這在一定程度上將顆粒材料破碎的現(xiàn)象研究發(fā)展到破碎過程的機理研究.Wang等[27-28]研究了顆粒材料破碎與顆粒材料剪切破壞行為之間的關(guān)系,如圖4 所示,通過微觀力學(xué)分析和機理演示,深入研究了平面應(yīng)變剪切條件下顆粒破碎對致密可壓試樣的應(yīng)力比、體應(yīng)變、塑性變形和剪切破壞行為的影響.首先,基于單個土壤顆粒真實斷裂行為對顆粒試樣進行建模,通過對比壓板壓縮模擬結(jié)果與物理實驗室試驗結(jié)果證明了模型的有效性,如圖5 所示;然后研究了影響顆粒材料破碎的主要因素,結(jié)果表明,體積膨脹率和峰值應(yīng)力比、塑性變形機制和剪切破壞模式隨土壤可壓性的變化是影響顆粒材料破碎的主要因素;最后,對破碎的機理進行了討論,從宏觀和微觀力學(xué)現(xiàn)象出發(fā),提出剪切帶和大規(guī)模體積收縮是致密試樣的兩種端部破壞模式,且這兩種模式分別由顆粒重排引起的膨脹和顆粒材料破碎引起的壓縮控制,在土壤可壓性和圍壓的中等范圍內(nèi),兩種破壞模式是組合和競爭關(guān)系,如圖6 所示[27].

圖3 單顆粒、三顆粒柱的SMT 和DEM 圖像以及不同載荷水平下實驗室規(guī)模的壓縮試驗[26]Fig.3 SMT and DEM images of single-particle and three-particle columns as well as laboratory-scale compression tests under different load levels[26]

圖4 具有膜型柔性側(cè)邊界的DEM 試件[27]Fig.4 DEM specimens with membrane-type flexible side boundaries[27]

圖5 單顆粒壓裂典型模擬結(jié)果[27]Fig.5 Typical simulation results of single particle fracturing[27]

圖6 黏結(jié)單元及其在20%軸向應(yīng)變下的顆粒位移和旋轉(zhuǎn)場[27]Fig.6 Bonding element and particle displacement and rotation field under 20% axial strain[27]

Wang 等[29]建立了水泥的橢球顆粒模型,采用黏結(jié)?破碎法分析了單顆粒的破碎機理,如圖7 所示.研究結(jié)果表明,高強度粒子碎裂成幾片,而低強度粒子碎裂成兩大塊.這說明顆粒的斷裂形態(tài)隨脆性的不同而變化,導(dǎo)致高強度顆粒出現(xiàn)脆性損傷,低強度顆粒出現(xiàn)延性損傷.對于高抗壓強度材料,由于顆粒脆性,損傷后無法承受載荷.對于低抗壓強度材料,由于顆粒的延展性,顆粒破碎不會立即完成,這些顆粒在損傷后仍能承受載荷.顆粒在剪應(yīng)力作用下膨脹,而孔隙率同時增加.這種剪應(yīng)力下的膨脹本質(zhì)是一種不穩(wěn)定的高勢能態(tài).顆粒的剪脹變形增大導(dǎo)致其結(jié)構(gòu)不穩(wěn)定,從而可能引起局部骨架結(jié)構(gòu)的坍塌.剪脹變形的積累是造成水泥橢球顆粒局部骨架結(jié)構(gòu)坍塌的原因之一.因此,在研究顆粒的破碎機理時,可通過單個顆粒的細觀和宏觀力學(xué)特性進行研究.

圖7 橢球體顆粒整體破碎過程演化(俯視圖,顆粒間線代表接觸黏結(jié))[29]Fig.7 Evolution of ellipsoidal particle disintegration process as a whole(top view,interparticle lines represent contact bonding)[29]

1.1.3 黏結(jié)?破碎法在工程中的應(yīng)用

由于黏結(jié)顆粒法可以構(gòu)造復(fù)雜的顆粒形態(tài),在實現(xiàn)顆粒破碎的模擬中可以對顆粒強度進行表征,因此,在分析此類工程問題時可以采用該方法.

鐵路道砟一般采用離散元方法建立數(shù)值模型,在列車的作用下會發(fā)生破碎.Lim 等[30]采用黏結(jié)顆粒法構(gòu)造道砟的三維非規(guī)則幾何形態(tài)顆粒,并進行了單顆粒壓碎試驗和側(cè)限壓縮試驗,研究結(jié)果與物理實驗進行對比從而證明模擬結(jié)果的準確性,如圖8 所示.

圖8 加載前黏結(jié)道砟顆粒的分布情況[30]Fig.8 Distribution of bonded ballast particles before loading[30]

對海冰與海洋結(jié)構(gòu)相互作用的研究,一般采用離散元方法建立海冰的數(shù)值模型,在海洋結(jié)構(gòu)的作用下海冰會發(fā)生破碎,因此,海洋工程中海冰的破碎是顆粒材料破碎的重要應(yīng)用領(lǐng)域.Lubbad 等[31]針對冰船相互作用過程,考慮了浮冰破碎后形成的小的浮冰大多會被推到一邊、旋轉(zhuǎn)或淹沒的復(fù)雜情況,建立了一個實時模擬船冰相互作用過程的數(shù)值模型,并將新的解析閉合解用于表示破冰過程.利用PhysX 求解計算域內(nèi)所有浮冰6 自由度剛體運動方程.圖9 為錐體結(jié)構(gòu)的破冰過程.其中浮冰采用多面體黏結(jié)單元模擬,破碎后的黏結(jié)單元不發(fā)生二次破碎,但也不會被刪除,將參與到離散元計算中并對結(jié)構(gòu)產(chǎn)生作用.針對海洋工程結(jié)構(gòu)物與海冰作用下產(chǎn)生的海冰破碎過程,季順迎團隊發(fā)展了球體[32]和多面體[33]離散元方法,建立了海冰的黏結(jié)破碎模型,并將其應(yīng)用于海洋工程的破冰過程中,如圖10所示.

圖9 錐體結(jié)構(gòu)的破冰過程[31]Fig.9 Ice breakage process of cone structure[31]

圖10 平整冰與海洋工程結(jié)構(gòu)物互相作用并發(fā)生破碎Fig.10 Interaction and fragmentation between flat ice and marine engineering structures

此外,徐永福[34]采用離散元軟件PFC-2D 通過黏結(jié)單元法構(gòu)造了粗粒土模型,并討論了直剪試驗中基礎(chǔ)顆粒間的黏結(jié)力、黏結(jié)顆粒的孔隙率和粗粒土試樣的孔隙率等因素對顆粒材料體應(yīng)變和剪切強度的影響.

1.1.4 黏結(jié)?破碎法的優(yōu)勢及存在的問題

黏結(jié)顆粒法能夠構(gòu)造復(fù)雜的顆粒形狀[15],并實現(xiàn)對顆粒破碎行為的力學(xué)模擬,再現(xiàn)了可破碎顆粒復(fù)雜的力學(xué)響應(yīng)[16],可以對顆粒強度的尺寸效應(yīng)和Weibull 分布特性進行一定程度的分析[35].然而,由于每個顆粒簇是由多個基礎(chǔ)顆粒黏結(jié)而成,模擬顆粒破碎過程將涉及大量的基礎(chǔ)顆粒從而影響計算效率.因此,基于黏結(jié)?破碎法的顆粒破碎分析方法在單顆?；蛐∫?guī)模顆粒集合體的破碎研究中具有一定的優(yōu)勢.黏結(jié)?破碎法中破碎后的顆粒無法發(fā)生二次破碎,這與真實的破碎過程不符.此外,在離散元模擬中,黏結(jié)顆粒的物理參數(shù)需要通過大量的試驗進行校正[14],在模擬顆粒破碎時需要預(yù)先設(shè)定破碎路徑,為了與真實破碎現(xiàn)象接近,需要大量的基礎(chǔ)顆粒數(shù)[36].因此,黏結(jié)?破碎法面臨的主要問題是計算效率問題,且該問題是無法通過方法優(yōu)化解決的,目前主要解決辦法是發(fā)展GPU 并行算法,通過提高計算機硬件來提高計算效率.

1.2 碎片替代法

1.2.1 碎片替代法的提出與發(fā)展

碎片替代法(fragment replacement method,FRM)是strm 等[37]在1998 年提出的.該方法主要實現(xiàn)過程如下: 當顆粒的受力情況滿足預(yù)設(shè)的破碎準則時,開始碎片替換,在將要發(fā)生破碎的顆粒所占空間中生成碎片并刪除原顆粒從而完成替換,如圖11 所示.替換后的子顆粒在非常短的時間步內(nèi)以較低的速度膨脹,這可能導(dǎo)致局部應(yīng)力增大.因此,需要控制子顆粒的膨脹量,這主要通過局部顆粒重疊量和顆粒速度來判定[38-39].

圖11 破碎替代法的破碎模式[39]Fig.11 Breakage mode of breakage substitution method[39]

碎片替代法研究顆粒破碎時,常用的破壞準則是八面體剪應(yīng)力破壞準則.該準則用八面體剪應(yīng)力來表示

式中,σ1,σ2,σ3分別為第1、第2、和第3 主應(yīng)力.離散元方法中一個顆粒的應(yīng)力張量 σij可定義為

式中,V為顆粒的體積,nc為該顆粒接觸的總數(shù),為接觸力,為接觸中心的支向量,σ1,σ2和 σ3則分別對應(yīng)于 σ11,σ22,σ33.在模擬過程中,當離散元計算得到的八面體剪應(yīng)力超過設(shè)定的容許八面體剪應(yīng)力q時,則認為顆粒破碎.De Bono 等[25]指出q=0.9σf=0.9Ff/d2,其中,σf為 顆粒的破碎強度,Ff為單顆粒徑向壓縮試驗中顆粒的峰值破碎力,d為顆粒的粒徑.根據(jù)Weibull 分布理論,每個顆粒的破碎強度 σf可表示為

式中,Ps(d),m,σ0,d0分別為顆粒存活概率、顆粒破碎強度的Weibull 模量、顆粒特征破碎強度和顆粒特征破碎強度對應(yīng)的顆粒粒徑.其中顆粒破碎強度的尺寸效應(yīng)通過 (d/d0)?3/m體現(xiàn).

?str?m 等[37]提出了3 個構(gòu)建碎片替代模式的標準: 盡量減少碎片數(shù)量,這是為了避免多個破碎發(fā)生后,碎片數(shù)目的迅速增加導(dǎo)致計算成本增大;合理布置碎片顆粒,以達到破碎時局部應(yīng)力突減的效果;替代破碎顆粒的碎片的粒徑分布應(yīng)盡可能與實際情況相符.?str?m 指出這3 個準則很難同時滿足.

Tsoungui 等[40]建立的碎片替換模式包含4 種尺寸12 個子顆粒,并采用二維離散元數(shù)值方法對側(cè)限壓縮試驗進行了模擬,結(jié)果顯示,該碎片替代模式能夠很好地模擬顆粒破碎的過程.Lobo-Guerrero等[41]采用碎片替代法,通過用許多不同大小顆粒的組合替換一個在拉力中失效的顆粒,模擬顆粒材料在不同雙軸應(yīng)力組合作用下的破碎演化,并將顆粒材料在這些條件下的破碎過程可視化.

此外,Ben-Nun 等[42]采用兩種新的替代模式對顆粒集合體在側(cè)限壓碎時的拓撲結(jié)構(gòu)對自組織特性的影響進行了分析.楊貴等[43]基于Ben-Nun 提出的替代模式發(fā)展了兩種新的碎片替代模式用于粗粒料破碎行為的研究.以上研究證明基于碎片替代法的顆粒破碎模擬可有效分析可破碎顆粒的力學(xué)響應(yīng).然而,以上均為二維替代模式,真實顆粒均處于三維環(huán)境中,且三維情況下顆粒材料的力學(xué)行為更加復(fù)雜,開展三維顆粒破碎研究具有更加重要意義.

近些年,三維情況下基于碎片替代法的顆粒破碎分析方法也得到了發(fā)展.例如,McDowell 等[44]通過在初始顆粒邊界內(nèi)放入有一定重疊量的軸對稱的等粒徑小顆粒的方式,構(gòu)建3 種碎片替代模式,用于側(cè)限壓縮試驗中顆粒細觀力學(xué)特性分析,如圖12 所示.Marketos 等[45]采用構(gòu)建的替代模式研究了砂巖顆粒的開裂過程,重點關(guān)注了壓縮帶的萌生和擴展過程.Tavares 等[46-47]根據(jù)應(yīng)力的施加位置人為設(shè)置母顆粒內(nèi)部的顆粒碎片,其中較大的碎片在垂直于引起破碎的應(yīng)力方向上,并且碎片顆粒在該方向上具有重疊量,較小的碎片排列在剩余的空間,通常較小的碎片與較大的碎片之間也有重疊量,如圖13 所示.初始狀態(tài)的重疊量是為了保證體積守恒,但這會引起母顆粒產(chǎn)生爆炸力,為了消除這個爆炸力,其中一個方法是使用松弛因子來限制由粒子重疊量計算得來的法向力.另一種簡單的方法是人為消除碎片生成時產(chǎn)生的重疊量,只允許母顆粒之間真正重疊量的存在.

圖12 碎片替代法研究壓縮作用下顆粒破碎的演化過程[44]Fig.12 Fragment substitution method for studying the evolution process of particle fragmentation under compression[44]

圖13 碎片的替代方案[47]Fig.13 Alternative scheme for fragments[47]

由于碎片替代法中替代顆粒的尺寸直接影響模擬效果,因此,眾多的專家學(xué)者對顆粒破碎過程的尺寸效應(yīng)進行了研究.Ciantia 等[48]采用碎片替代法進行了側(cè)限壓縮數(shù)值試驗,研究了顆粒破碎的尺寸效應(yīng).Zhou 等[49]借鑒等粒徑三碎片替代方法研究了主應(yīng)力對可破碎顆粒材料宏細觀力學(xué)響應(yīng)的影響.為了使碎片粒徑分布盡可能地符合實際情況,Zhou 等[50]又構(gòu)建了新的三維碎片替代模式,對可破碎顆粒破碎過程的尺寸效應(yīng)進行了研究.Cil 等[51]采用所提出的替代模式,進一步研究了顆粒破碎能量的尺寸效應(yīng).通過以上研究,對顆粒破碎過程中的尺寸效應(yīng)有了更深入的理解和探索.Tavares 等[46]認為顆粒的破碎程度可以通過一個參數(shù)t10來表征,t10表示小于母顆粒尺寸 1/10 的碎片的比例,并通過以下表達式來計算顆粒的破碎程度,其中,A和b′為由實驗數(shù)據(jù)擬合的模型參數(shù),E50b為發(fā)生破碎的顆粒的斷裂能的中位數(shù),為破碎過程中顆?？捎玫哪芰?t10值越高,子代顆粒尺寸分布越細.

EDEM 在考慮單顆粒破碎時是基于碎片替代法實現(xiàn)的.EDEM 在研究單顆粒破碎過程時可直接獲取粒徑分布、破碎率等數(shù)據(jù),直觀展示破碎效果,模擬速度非?？?但模擬結(jié)果目前無法體現(xiàn)破碎后的物料料型.

1.2.2 碎片替代法的破碎準則研究

目前,采用碎片替代法研究顆粒破碎過程時,破碎準則并沒有統(tǒng)一的標準.現(xiàn)存在的破碎準則均基于一定的假設(shè)并做了相應(yīng)的簡化.?str?m 等[37]認為基于碎片替代法的破碎準則主要可以歸結(jié)為兩類:一是基于應(yīng)力的判定準則,即通過對顆粒的應(yīng)力設(shè)定閾值,當應(yīng)力大于所設(shè)定的閾值時顆粒發(fā)生破碎;二是基于力的判定準則,即通過對顆粒的接觸力設(shè)定閾值,當最大接觸力大于所設(shè)定的閾值時顆粒發(fā)生破碎.其他的破碎準則均是在這兩種準則的框架下發(fā)展起來的.

假設(shè)最大接觸力引起顆粒局部應(yīng)力集中從而導(dǎo)致顆粒破碎[52-53],那么最大接觸力可以作為破碎準則.Lobo-Guerrero 等[41,54-55]在二維模擬中對顆粒受力情況進行了簡化,認為當接觸數(shù)目不大于3 時,顆粒的受力情況可以類比于巴西圓盤試驗.基于以上假定,采用顆粒的特征拉應(yīng)力 σt=2P1/(πLD) 作為顆粒破碎的判定條件,即 σt＞σtmax,也 即P1＞σtmaxπLD/2,其中P1為作用力,L為單位長度,D為盤直徑.

由于配位數(shù)對顆粒破碎有一定的影響[40,56],在力判定準則中引入配位數(shù)的影響,以上假定不再適用,即配位數(shù)大于3 時顆粒仍可能發(fā)生破碎[57].根據(jù)物理試驗結(jié)果建立球形顆粒在徑向加載時的最大流動剪切強度公式[58],當配位數(shù)在6～12 時該公式同樣適用[59].通過以上公式Ciantia 等[48]求得顆粒最大接觸力閾值,進而建立三維模擬中顆粒破碎的判定準則,即kmob≤k,其中kmob和k分別為顆粒的流動強度和固有強度,當流動強度大于固有強度時顆粒發(fā)生破碎.選取3 種不同破碎難易程度的材料,張科芬等[60]對其進行破碎模擬,證明以上破碎準則適用于離散元中球形顆粒的破碎模擬,并且能夠反映顆粒的真實破碎過程.Cil 等[51]基于顆粒材料特征強度的3 種尺寸效應(yīng)公式發(fā)展了3 種不同的力判定標準,并通過側(cè)限固結(jié)試驗研究了采用上述3 種力判定準則時顆粒宏細觀力學(xué)行為的不同.

根據(jù)均勻化思想發(fā)展了應(yīng)力判定準則假設(shè)外載荷作用下顆粒內(nèi)部的應(yīng)力均勻分布,是顆粒受力狀態(tài)的整體考量.半徑為R的顆粒受到任意力作用時,其平均應(yīng)力張量為是顆粒的接觸數(shù)目,分別為接觸點(c)處的單位法向量沿i方向的分量和沿j方向的接觸力,Vp為顆粒的體積.其中,以等效力張量為破碎準則的本質(zhì)是通過顆粒的特征強度判定顆粒是否發(fā)生破碎[57].顆粒的受力情況可以簡化為靜水壓力和偏應(yīng)力,采用等效二維Drucker-Prager 準則判定顆粒的破碎情況[40].然而,由于顆粒間的接觸十分復(fù)雜,當存在多個接觸力時顆?？赡芴幱诟哽o水壓力、低偏應(yīng)力的狀態(tài),在這種情況下顆粒的破碎現(xiàn)象不容易發(fā)生,此時平均應(yīng)力判定準則不在適用.針對以上情況,八面體剪應(yīng)力破碎準則被提出并通過單顆粒壓碎試驗獲得該準則中特征拉應(yīng)力的閾值[44,61];帶截斷的Mohr-Coulomb 準則也可有效解決平均應(yīng)力判定準則不適用的情況[49-50,62];楊貴等[43]通過對顆粒的拉裂破壞和剪切破壞情況提出了不同的破碎準則.

以上兩種判定準則在一定程度上可以反映顆粒的破碎機制,然而,由于其均在一定假設(shè)的前提下提出的,且進行了一定的簡化,所以無法完全描述顆粒復(fù)雜的破碎行為.此外,對破碎準則的選取并沒有特定的標準,選取不同的破碎準則獲得的顆粒系統(tǒng)孔隙率和級配的演化過程的正確性也并未得到驗證.此外,二維和三維問題存在較大差異.因此.針對三維問題提出一套全面的、準確的選取破碎準則的方法是十分必要的.

1.2.3 碎片替代法在工程中的應(yīng)用

由于碎片替代法可以實現(xiàn)二次破碎,且對于不考慮體積守恒和對破碎顆粒形態(tài)要求不嚴格的工程問題,可以采用該方法實現(xiàn)顆粒的破碎研究,同時,碎片替代法在大規(guī)模的破碎過程問題研究中也有較大的優(yōu)勢,因此在以下工程問題中得到廣泛應(yīng)用.

碎片替代法在巖土工程中具有廣泛的應(yīng)用,如通過模擬巖體破碎的過程來預(yù)測巖體的破碎區(qū)域、巖塊大小分布、巖體的穩(wěn)定性等[63];通過模擬地下開挖過程中巖土體的破碎行為,預(yù)測巖土體的破碎區(qū)域、位移和變形情況,這對于開挖工程的設(shè)計和施工具有重要意義[64];該方法還應(yīng)用于模擬巖土體在振動載荷下的破碎行為,如在振動臺試驗中模擬土體的破碎過程,研究振動對土體強度和穩(wěn)定性的影響;在模擬地下爆破引起的巖土體破碎過程中可獲得爆破對周圍巖土體的影響,如爆破震動對周圍地基的影響、爆破引起的巖層破碎區(qū)域等;此外,碎片替代法可以用于預(yù)測巖土材料的力學(xué)性質(zhì),如彈性模量、斷裂強度等,這對于巖土工程設(shè)計和材料選取具有重要意義[65].

Lobo-Guerrero 等[55]發(fā)展的碎片替代模式包含3 種不同大小共8 個子顆粒,并通過二維離散元數(shù)值方法研究了由于顆粒材料破碎引起的鐵路道砟級配的變化對鐵沉降變形的影響,如圖14 所示.

圖14 碎片替代法在鐵路道砟顆粒材料破碎中的應(yīng)用[55]Fig.14 Application of fragment substitution method in railway ballast particle breakage[55]

張科芬等[60]基于局部應(yīng)力集中的點載荷破碎準則,采用碎片替代法開展了3 種不同破碎難易程度材料的數(shù)值模擬.其中,為了保證破碎前、后顆粒之間的平衡,利用阿波羅填充和膨脹法,并引入尺寸因子來表征不同粒徑的顆粒強度.

1.2.4 碎片替代法的優(yōu)勢及存在的主要問題

基于碎片替代法的顆粒破碎數(shù)值方法在一定程度上再現(xiàn)了顆粒的破碎過程,且該方法能夠?qū)崿F(xiàn)顆粒的二次破碎.碎片替代法為較大規(guī)模顆粒破碎過程分析和細觀尺度的機理研究提供了一定的思路.然而,由于碎片替代模式的主觀性導(dǎo)致模擬結(jié)果與真實的顆粒形態(tài)及破碎過程存在較大的差異;三維模擬中顆粒破碎的判定準則并不完善.對以上問題的有效解決將促進碎片替代法的發(fā)展和應(yīng)用.

2 基于離散元?有限元耦合方法的顆粒破碎數(shù)值算法

當參考連續(xù)介質(zhì)力學(xué)對材料破碎的研究方法并結(jié)合顆粒材料的離散特性,可采用離散元?有限元耦合的方法對顆粒破碎進行數(shù)值研究.以下將介紹比例邊界有限元方法、組合有限元?離散元方法以及典型的內(nèi)聚力模型.

2.1 比例邊界有限元法

2.1.1 比例邊界有限元法的提出與發(fā)展

比例邊界有限元法是由Song 等[66]提出的一種新的半解析方法,該方法是基于有限單元和邊界元方法提出的,同時擁有有限元法和邊界元法的優(yōu)點,已經(jīng)在一些工程領(lǐng)域中的動力學(xué)和靜力學(xué)問題研究中得到了應(yīng)用.比例邊界有限元法相對于有限元法和邊界元法計算過程更簡單有效,并可得到工程滿意的結(jié)果.Ooi 等[67]發(fā)展了適用于任意邊數(shù)的凸多邊形比例邊界有限元法.凸多邊形比例邊界有限元僅用一個多邊形單元就可以模擬具有任意邊數(shù)的凸多邊形顆粒,并得到其變形和應(yīng)力狀態(tài).

比例邊界有限元法實現(xiàn)顆粒的破碎過程時需要對每個顆粒的應(yīng)力?應(yīng)變進行完整分析,然后根絕準則來判定顆粒內(nèi)部的破壞點,當破壞點的數(shù)量達到一定比例后,認為該顆粒發(fā)生破碎.比例邊界有限元法認為顆粒破碎的路徑是直線,所以采用加權(quán)最小二乘法對破壞點進行擬合,將破碎點連接起來獲得破碎路徑,從而使顆粒破碎為兩部分.該方法不需要對破壞路徑的位置和方向進行預(yù)定義.破碎后產(chǎn)生的顆粒在下一時間步中參與計算.羅滔等[36]采用比例邊界有限元方法與離散元方法結(jié)合,對堆石料的雙軸壓縮實驗進行了模擬,結(jié)果如圖15 所示.

圖15 加載結(jié)束后不同圍壓下的顆粒破碎分布[36]Fig.15 Particle fragmentation distribution under different confining pressures after loading[36]

2.1.2 比例邊界有限元法的破碎準則

為了研究顆粒的破碎行為,需要建立相應(yīng)的破碎準則,由于比例邊界有限元方法已經(jīng)計算獲得了每個顆粒內(nèi)部的應(yīng)力場,因此,采用基于應(yīng)力的破壞準則判斷顆粒是否破壞.判斷破壞的流程為: 首先在顆粒內(nèi)部選取具有代表性的樣本點,如圖16 所示,基于Hoek-Brown 破碎準則判斷樣本點是否破壞;當破壞樣本點的數(shù)量達到一定比例時顆粒破壞;比例邊界有限元方法假設(shè)顆粒的破壞路徑為直線,并通過加權(quán)最小二乘法對破壞點進行擬合得到破壞路徑;沿著破壞路徑將顆粒分為兩個顆粒[68].

圖16 樣本點的分布Fig.16 Distribution of sample points

為了判斷樣本點是否破壞,應(yīng)用Hoek-Brown 破壞準則[68]

式中,σt為抗拉強度;σ1f是由式(4)得出的第1 主應(yīng)力,如果安全系數(shù)Fs＜1 則表示此點破壞,由Hoek-Brown 準則得到的破壞點如圖17 所示.

圖17 由Hock-Brown 準則得到的破壞點Fig.17 Failure points obtained from the Hock-Brown criterion

2.1.3 比例邊界有限元法的優(yōu)勢及存在的主要問題

在一定程度上,比例邊界有限元方法從細觀尺度揭示了顆粒破碎機制.比例邊界有限元方法計算顆粒內(nèi)部的應(yīng)力分布,從而實現(xiàn)顆粒的破碎過程,顆粒破碎產(chǎn)生的新的顆粒將直接參與到離散元和比例邊界有限元的計算中,比例邊界有限元方法對顆粒破碎的研究無需預(yù)定義任何子顆?；蛘呔W(wǎng)格重劃分[36].然而,該方法需要人為設(shè)計樣本點的分布,且其破碎路徑均為直線,因此,比例邊界有限元法的主觀性很強,很難客觀反映顆粒破碎真實過程,發(fā)展從細觀尺度反映顆粒破碎真實過程的數(shù)值方法迫在眉睫.

2.2 組合有限元?離散元方法

2.2.1 組合有限元?離散元方法的提出與發(fā)展

組合有限元?離散元方法(FDEM)是Munjiza等[69]在 20 世紀90 年代為解決混凝土加載中裂紋的萌生和擴展問題發(fā)展的數(shù)值方法.組合有限元?離散元方法首先對顆粒進行單元劃分,并通過接觸面單元對其進行連接,離散元法計算獲得顆粒之間的作用力作為每個顆粒應(yīng)力計算的邊界條件,從而實現(xiàn)單個顆粒的應(yīng)力?應(yīng)變分析,獲得顆粒的應(yīng)力狀態(tài)和變形情況.組合有限元?離散元法中顆粒沿有限單元的邊界破壞.

FDEM 采用FEM 理論作為顆粒的變形、壓裂破壞標準,通過DEM 概念來檢測新的接觸點,解決離散單元的平移、旋轉(zhuǎn)和相互作用問題[70].FDEM利用分布接觸力方法和罰函數(shù)方法解決接觸力問題,有效地處理了復(fù)雜形狀的變形問題[71].FDEM 在非線性彈性斷裂力學(xué)(NLEFM)框架下模擬裂縫的起裂和擴展過程,此外,在FDEM 中可以實現(xiàn)重網(wǎng)格技術(shù),但在建?？善扑轭w粒材料時,計算成本非常高.由于在模擬過程中沒有執(zhí)行重網(wǎng)格,網(wǎng)格拓撲也沒有更新,因此在加載過程中粒子不能繼續(xù)變小.

基于FDEM 發(fā)展的顆粒破碎仿真軟件主要有中科院力學(xué)所李世海團隊的GDEM,該軟件采用有限元和塊體離散元耦合的方法對巖土工程中的顆粒材料破碎問題實現(xiàn)了數(shù)值分析.基于多物理場斷裂分析軟件MultiFracS 的基本理論也是有限元?離散元耦合算法,該軟件在巖土材料斷裂及破碎數(shù)值建模方面有獨特優(yōu)勢,在隧道、邊坡、采礦、水利水電、爆破、非常規(guī)油氣及地熱開采等領(lǐng)域有重要應(yīng)用價值.

2.2.2 組合有限元?離散元方法在工程中的應(yīng)用

由于FDEM 方法在模擬顆粒破碎的過程時,對每個顆粒進行應(yīng)力分析,因此破碎過程與實際更加相符.在對單個顆粒的破碎分析中可有效揭示顆粒的破碎機理.因此針對考慮破碎機理的工程問題,以及需要模擬的破碎過程更加接近真實情況時,可采用該方法實現(xiàn).

組合有限元?離散元方法在巖土工程中得到廣泛的應(yīng)用,如嚴成增等[72]將其與數(shù)值圖像技術(shù)結(jié)合,系統(tǒng)分析了巖體的破裂機制;Ma 等[73]采用組合有限元?離散元方法建立了巖石顆粒的數(shù)值模型,如圖18 所示,并對巖石顆粒在軸壓作用下的斷裂過程進行了模擬,系統(tǒng)研究了破碎前顆粒的形貌對破碎后碎片的形態(tài)和尺寸分布的影響,如圖19 所示.此外,在組合有限元?離散元方法中引入基于應(yīng)力的破壞準則,并對堆石料的破碎過程進行模擬,模擬過程中無需預(yù)先定義或者假設(shè)顆粒的破壞路徑[74-75].Mahabadi等[76]使用三維FDEM 方法研究Opalinus clay 黏土頁巖地層中隧道周圍的開挖損傷區(qū)的發(fā)展.通過巴西圓盤試驗和單軸壓縮試驗,對三維FDEM 模型進行校準,結(jié)果如圖20 所示.

圖18 巖石顆粒的FDEM 數(shù)值模型[73]Fig.18 Numerical model of rock particles[73]

圖19 FDEM 方法研究巖石顆粒在軸壓作用下的斷裂結(jié)果[73]Fig.19 FDEM method for studying the fracture process of rock particles under axial compression[73]

圖20 采用FDEM 方法進行巴西圓盤和單軸壓縮試驗[76]Fig.20 Brazilian disc and uniaxial compression tests using FDEM method[76]

Rougier 等[77]將離散裂紋與塑性變形相結(jié)合,并采用FDEM 方法模擬二維泰勒桿的斷裂過程,如圖21 所示.其中變形采用FDEM 方法描述,塑性變形在材料嵌入坐標系中表述,乘法分解和塑性流動在拉伸空間中求解.該工作將FDEM 中的裂紋和破碎標準有效結(jié)合.

圖21 FDEM 方法模擬二維泰勒桿的斷裂過程[77]Fig.21 FDEM method simulation of the fracture process of twodimensional Taylor bars[77]

針對脆性材料的斷裂和破碎響應(yīng),Chen 等[78]采用FDEM 方法對夾膠玻璃在硬體沖擊下的斷裂和破碎過程進行了模擬,如圖22 所示,并討論了玻璃、層間和層間界面的破壞模型,證明了該模型在模擬夾層玻璃破裂時是可靠的.Lilja 等[79]采用FDEM 方法建立海冰的三維數(shù)值模型,如圖23 所示,該模型由同轉(zhuǎn)動、黏性阻尼Timoshenko 梁有限元組成的平面梁晶格與形成實際海冰的剛性離散單元的質(zhì)量質(zhì)心相連接.通過質(zhì)心Voronoi 鑲嵌程序生成網(wǎng)格,由于內(nèi)部存在阻尼,基于晶格的結(jié)構(gòu)的機械響應(yīng)依賴于應(yīng)變速率和尺寸.

圖22 夾膠玻璃在硬體沖擊下的斷裂和破碎響應(yīng)[78]Fig.22 Fracture and fragmentation response of laminated glass under hard impact[78]

圖23 采用FDEM 方法建立海冰的三維數(shù)值模型[79]Fig.23 Establishing a three-dimensional numerical model of sea ice using FDEM method[79]

2.2.3 組合有限元?離散元方法的優(yōu)勢及存在的主要問題

FDEM 方法可有效模擬顆粒的破碎過程,由于該方法對每個顆粒進行應(yīng)力分析,因此破碎過程與實際更加相符.在對單個顆粒的破碎分析中可有效揭示顆粒的破碎機理.然而,由于模擬過程中需要對每個顆粒進行有限單元劃分,且顆粒的破碎只能沿著網(wǎng)格進行,為了獲得更為真實的破碎過程,需要網(wǎng)格精細化,這均導(dǎo)致計算成本非常大[36].這將限制了其在工程中的應(yīng)用,對實際工程問題的模擬無法起到推動作用.

2.3 內(nèi)聚力模型

2.3.1 內(nèi)聚力模型的提出與發(fā)展

FDEM 方法中,有限元離散化以后,將黏性界面單元插入到與每個粒子相關(guān)的有離散單元中,以考慮粒子的破碎過程.本節(jié)將重點介紹黏性界面的內(nèi)聚力模型(cohesive zone model,CZM)[9,80-82].內(nèi)聚力模型是1960 年Dugdale 等提出的[83-84],假設(shè)裂紋的尖端是由兩個裂紋界面組成的很小的內(nèi)聚區(qū),該內(nèi)聚區(qū)的本構(gòu)關(guān)系就是界面上作用牽引力T與兩裂紋面間相對位移U之間的關(guān)系.內(nèi)聚力模型和裂紋尖端內(nèi)聚區(qū)的分布[85]如圖24 所示.

圖24 內(nèi)聚力模型和裂紋尖端內(nèi)聚區(qū)[85]Fig.24 Cohesion model and crack tip cohesive zone[85]

內(nèi)聚力模型從未完全承載的點A開始,T隨著U的增加而增加,隨之達到應(yīng)力最大值Tmax的點C,此時材料的應(yīng)力承載達到最大值,材料在該點處開始出現(xiàn)初始損傷.隨著界面位移的繼續(xù)增大,應(yīng)力開始下降,該階段為材料的損傷擴展階段,點E為材料裂紋界面完全分離點,此時材料的承載力降為0[85].內(nèi)聚區(qū)應(yīng)力的變化受內(nèi)聚力模型和裂紋界面位移的影響,因此,針對不同的材料,應(yīng)該選擇不同的內(nèi)聚力模型,同時通過選取適當?shù)膮?shù),從而反映界面的強度、韌度等力學(xué)性能.

內(nèi)聚力模型可分為基于位移的內(nèi)聚力模型和基于勢能的內(nèi)聚力模型.在基于位移的內(nèi)聚力模型中,裂紋上下表面之間的有效牽引力可由有效分離位移的函數(shù)表示,通常稱之為牽引力分離法則.常用到的牽引力分離法則包括線性軟化、雙線性軟化、指數(shù)和梯形法則等.不同法則的區(qū)別在于與之間函數(shù)關(guān)系的不同.將有效牽引力與內(nèi)聚強度σmax歸一化處理的結(jié)果如圖24(b)所示.基于位移的內(nèi)聚力模型出現(xiàn)的裂紋擴展問題,可以在基于勢能的通用內(nèi)聚力模型中得到解決.基于勢能的內(nèi)聚力模型采用三次多項式表示法向牽引力,采用線性關(guān)系式表示切向牽引力[85].

早期非線性斷裂研究一般認為,當內(nèi)聚區(qū)尺寸小于裂紋和試樣尺寸時,內(nèi)聚力模型理論與Griffith的能量平衡理論等效.Dugdale[83]認為內(nèi)聚應(yīng)力的分布在數(shù)值上與材料的屈服強度相等,但這有悖于實際物理事實,Barenblatt[84]考慮內(nèi)聚力分子尺度的特征,把內(nèi)聚應(yīng)力看作內(nèi)聚區(qū)裂紋面各點處裂紋張開位移的函數(shù),但具體的解析式無法準確給出,在大多數(shù)情況下內(nèi)聚力仍被假設(shè)為常數(shù).基于以上工作,Hillerborg 等[86]考慮顆粒的拉伸強度,采用有限元方法建立脆性材料的內(nèi)聚力模型并對其斷裂過程進行模擬,實現(xiàn)了原有裂紋的增長、新裂紋的萌生和擴展研究,對斷裂細節(jié)進行了完整的描述.Needleman[87]通過高次多項式函數(shù)對延性材料的斷裂情況進行了模擬.Kolher 等[88]采用了分段函數(shù)描述內(nèi)聚力模型的方法對鎳鋁合金的剪切斷裂性能進行了模擬.

為模擬脆性巖體的斷裂行為,清華大學(xué)徐文杰團隊[89]提出了一種針對多面體離散元顆粒設(shè)計的內(nèi)聚斷裂模型(CFM),如圖25 所示.巖體被離散成一系列剛性多面體塊體,這些塊體沿邊界面的法向和剪切方向黏結(jié)在一起.黏性準則規(guī)定了CFM 的法向和剪切斷裂強度.通過在Blaze-DEM 框架中對多個圖形處理單元的并行處理,提高了CFM 在多面體模擬中的計算效率.

圖25 斷裂前、后的接觸模型[89]Fig.25 Model after fracture[89]

2.3.2 內(nèi)聚力模型在工程中的應(yīng)用

內(nèi)聚力模型在處理非線性和大變形問題時具有非常強的適應(yīng)性.隨著內(nèi)聚力模型的改進和發(fā)展,在多種材料的斷裂問題上具有獨特的優(yōu)勢,因此,在處理此類工程問題時可采用該方法.

Gui 等[90]在通用離散元編碼(UDEC)中,引入了結(jié)合拉伸、壓縮和剪切材料特性的內(nèi)聚斷裂模型,用于模擬巖石動力試驗中的壓裂過程,通過缺口半圓彎曲試驗和巴西盤試驗兩個數(shù)值算例,驗證了內(nèi)聚力模型對巖石動態(tài)破壞模擬的有效性.

為了模擬巖石破裂過程,廣東工業(yè)大學(xué)高偉團隊提出一種新的內(nèi)聚模型本構(gòu).在該本構(gòu)中,采用Mohr-Coulomb 準則作為裂紋萌生準則,建立壓縮和拉伸狀態(tài)下的牽引?分離模型,描述了巖石強壓力依賴性的力學(xué)特性.在粘結(jié)區(qū)完全破壞后,采用基于庫侖公式的摩擦模型來確定裂紋自由面之間的接觸相互作用.采用所提出的零厚度黏合單元本構(gòu),模擬了巴西盤、單軸壓縮兩個巖石破裂實例.比較仿真結(jié)果與實驗結(jié)果發(fā)現(xiàn)以上模型是可靠的[9],如圖26～圖27 所示.

圖26 0 厚度黏合單元本構(gòu)[9]Fig.26 Constitutive law of zero thickness adhesive element[9]

圖27 采用內(nèi)聚力單元模擬巴西盤破裂和單軸壓縮巖石破裂實例[9]Fig.27 Simulation of brazilian plate rupture and example of rock fracture using cohesive elements[9]

2.3.3 內(nèi)聚力模型的優(yōu)勢及存在的主要問題

內(nèi)聚力模型的理論基礎(chǔ)是彈塑性斷裂力學(xué)理論,該模型認為斷裂問題是非線性邊值問題,不必考慮斷裂的起裂和擴展準則,因此,適應(yīng)性強,可以解決非線性和大變形問題.隨著計算機技術(shù)和有限元方法的發(fā)展,內(nèi)聚力模型得到了改進和發(fā)展,已經(jīng)可以解決多種材料的斷裂問題.然而,針對顆粒型材料,仍然存在計算效率的問題,因此,該方法在研究顆粒型材料的破碎問題中仍沒有得到推廣.

3 近場動力學(xué)方法及其他方法

3.1 近場動力學(xué)方法

3.1.1 近場動力學(xué)方法的提出與發(fā)展

近場動力學(xué)(peridynamics,PD)是一種新興的非局部性理論,針對模型中的非連續(xù)問題,該方法不需要任何處理,對模型中從連續(xù)到非連續(xù)的過程可以很好地處理[91-92],因此,該方法在研究離散介質(zhì)的破碎問題時有一定的優(yōu)勢.近場動力學(xué)方法是Silling等[93]于2004 年提出的,該方法基于長程非局部作用思想建立力學(xué)模型,通過求解空間積分型運動方程來描述物質(zhì)力學(xué)行為,實現(xiàn)了對連續(xù)介質(zhì)和非連續(xù)介質(zhì)力學(xué)行為的統(tǒng)一描述.此外,由于損傷概念的引入,近場動力學(xué)理論對材料的損傷和裂紋的形成過程可以有效描述.因此,在處理裂紋自發(fā)萌生和擴展等非連續(xù)問題時,近場動力學(xué)理論具有一定的優(yōu)勢[94].

近場動力學(xué)方法將固體材料離散為一系列具有體積和質(zhì)量的物質(zhì)點,每個物質(zhì)點與其周圍一定區(qū)域內(nèi)所有其他物質(zhì)點之間存在相互作用力.如圖28所示.

圖28 近場動力學(xué)理論的局部接觸式作用和長程非局部作用[91]Fig.28 Local contact type action and long-range non local action in peridynamics theory[91]

近場動力學(xué)理論是一種非局部連續(xù)介質(zhì)力學(xué)理論,采用積分的形式表示鄰域半徑 δ 內(nèi)物質(zhì)點x和x′之間的相互作用力,即鍵力.物質(zhì)點x的動力學(xué)方程為

近場動力學(xué)理論通常被分為鍵型(bond-based)和態(tài)型(state-based)兩種.其中鍵型近場動力學(xué)理論認為兩物質(zhì)點之間的相互作用僅包含兩物質(zhì)點間的鍵力.而態(tài)型近場動力學(xué)理論在分析兩物質(zhì)點間的相互作用時考慮了鄰域內(nèi)所有物質(zhì)點對鍵力的影響.顯然,鍵型近場動力學(xué)理論是態(tài)型近場動力學(xué)理論的一種特例[91].

計算程序編寫促進了近場動力學(xué)理論的發(fā)展和應(yīng)用.Silling 團隊在1998 年采用Fortran 90 語言開發(fā)了近場動力學(xué)方法的代碼EMU,并對其進行了測試,這是第一次將近場動力學(xué)理論程序化;Parks 等基于分子動力學(xué)軟件LAMMPS 開發(fā)了軟件PDLAMMPS,擴展了近場動力學(xué)理論的應(yīng)用;Macek 等在經(jīng)典動力學(xué)分析軟件ABAQUS 中引入近場動力學(xué)算法,促進了近場動力學(xué)理論的發(fā)展;此外,LS-DYNA 中也引入了近場動力學(xué)方法擴展程序包,并命名為DYNAPeridynamics,有效地結(jié)合了現(xiàn)有商業(yè)軟件和近場動力學(xué)方法在處理動力學(xué)問題的優(yōu)勢;近期,桑迪亞國家實驗室基于離散單元法開發(fā)了開源軟件Peridigm,該軟件兼容了Cubit 網(wǎng)格生成工具和Paraview 可視化程序,將近場動力學(xué)方法平臺化,具有人機交互性好、計算穩(wěn)定、后處理方便、實用性強等優(yōu)點[91];趙吉東團隊[95]開發(fā)了PD-NSCD 混合法的代碼,并應(yīng)用于顆粒材料連續(xù)破碎過程的研究,自主開發(fā)代碼實現(xiàn)近場動力學(xué)理論部分,使用開源代碼Project Chrono 進行基于NSCD 的離散建模.同時,通過MATLAB,Python,C 語言,Fortran 等語言進行近場動力學(xué)數(shù)值計算與分析也是部分學(xué)者在處理實際問題時采用的手段[96-97].

趙吉東團隊[95,98-99]開發(fā)的PD-NSCD 混合法,采用非光滑接觸動力學(xué)模擬顆粒系統(tǒng)的剛體運動以及顆粒間相互作用,利用近場動力學(xué)理論分析單個顆粒的破碎.基于近動力學(xué)和物理引擎提出一種新的混合計算框架來模擬機械載荷下的可破碎顆粒材料.在此框架下,利用近場動力學(xué)方法模擬單個粒子的破碎,利用基于非光滑接觸動力學(xué)方法的物理引擎模擬粒子的剛體運動和粒子間的相互作用.該混合框架可以對顆粒破碎過程進行模擬,并實現(xiàn)連續(xù)破碎過程中不規(guī)則顆粒形狀的表征和模擬,克服了許多現(xiàn)有方法面臨的缺陷和挑戰(zhàn).其工作為研究顆粒破碎的復(fù)雜過程提供了一種有效的方法,為今后研究可破碎顆粒材料的微觀結(jié)構(gòu)行為提供了參考.圖29 顯示了近場動力學(xué)分析初始化和構(gòu)造粒子的完整過程[98],其中,圖29(a)為受接觸力作用的顆粒,圖29(b)為對離散粒子進行的近動力學(xué)分析,圖29(c)未顆粒破碎的近動力學(xué)分析結(jié)果,圖29(d)為物理引擎中的模擬粒子,圖29(e)為破碎粒子的拆分視圖.圖30～圖31 為砂石顆粒連續(xù)破碎過程的數(shù)值模擬結(jié)果[95,98].

圖29 顆粒破碎建模過程示意圖[98]Fig.29 Schematic diagram of particle breakage modeling process[98]

圖30 球形顆粒模擬砂石顆粒連續(xù)破碎過程[98]Fig.30 Simulation of continuous breakage process of sand and stone particles with spherical particles[98]

圖31 多面體顆粒模擬砂石顆粒連續(xù)破碎過程[95]Fig.31 Simulation of continuous breakage process of sand and gravel particles using polyhedral particles[95]

3.1.2 近場動力學(xué)方法在工程中的應(yīng)用

目前,近場動力學(xué)理論已經(jīng)在多個領(lǐng)域得到了應(yīng)用[100].如混凝土壓碎[93,101]、金屬沖擊[102]、子彈穿透[103-104]、土體[105-107]和陶瓷破碎[108]、薄膜撕裂[109]以及玻璃沖擊[110]等.目前,特別是在巖石力學(xué)研究中,近場動力學(xué)方法已經(jīng)成功模擬了水壓致裂破壞過程[111],對于預(yù)制裂紋的萌生、擴展和破壞過程模擬,近場動力學(xué)方法也是有效的.通過近場動力學(xué)理論,Ha 等[112]、朱其志等[113]和王振宇等[114]研究了單條預(yù)制裂紋在單軸壓縮條件下裂紋的演化過程;Lee等[115-116]研究了含有預(yù)制裂紋巖石的翼型裂紋和次生裂紋的形成過程;Rabczuk 等[117]采用雙鄰域近場動力學(xué)模型對含多條隨機裂紋的巖石的破壞過程進行了研究.此外,Zhou 等[118-121]引入摩爾庫倫準則和極限張拉破壞準則,對預(yù)制裂紋的巴西圓盤[120-121]、三點彎曲[119,121]、單軸壓縮[119]以及雙軸張拉[118]進行了數(shù)值模擬.

Hu 等[110]采用近場動力學(xué)方法分析了玻璃板的損傷模式,其中沖擊速度從61～200 m/s 不斷增加.將實驗結(jié)果與簡化模型的近場動力學(xué)模擬結(jié)果進行了比較,如圖32 所示.實驗觀察到的主要斷裂模式是由3 種沖擊速度的每一種測試的近場動力學(xué)模型捕獲的.該方法通過進一步改進邊界條件以實現(xiàn)對多層體系的沖擊造成的脆性損傷更準確的模擬.近場動力學(xué)計算模型揭示了玻璃層復(fù)雜脆性損傷演化的早期階段,以及邊界條件對動態(tài)破裂過程的影響.

圖32 撞擊速度為61 m/s 的破碎圖,以上分別是玻璃正面、底面和實驗結(jié)果[110]Fig.32 Fragmentation diagram with an impact velocity of 61 m/s,showing the front,bottom,and experimental results of the glass[110]

Rabczuk 等[117]提出了一種用于顆粒狀和巖石狀材料斷裂的雙層近場動力學(xué)(DH-PD)方法.與離散裂紋方法相比,DH-PD 方法不需要任何裂紋表面的表示形式,也不需要任何標準來處理復(fù)雜的裂紋形態(tài),裂紋路徑是模擬的自然結(jié)果.在該方法中作者提出了一種新的模擬壓縮裂縫接觸和約束侵徹條件的懲罰方法.該方法被應(yīng)用于地質(zhì)力學(xué)的一些基準問題,圖33 為采用雙層近場動力學(xué)方法研究方形板裂紋擴展的數(shù)值模擬結(jié)果.

圖33 雙層近場動力學(xué)研究方形板中裂紋的擴展[111]Fig.33 Double layer Peridynamics study on crack propagation process in square plates[111]

Shen 等[101]采用近場動力學(xué)理論分析混凝土結(jié)構(gòu)的損傷和漸進破壞.建立了混凝土柱單軸受壓非局部近動力模型,如圖34 所示.模擬結(jié)果表明,混凝土近動力體的裂縫是自發(fā)形成的.此外,模擬結(jié)果也展現(xiàn)了基于近場動力學(xué)理論建?；炷翐p傷累積和漸進破壞的優(yōu)勢.Shen 等的研究為分析復(fù)雜的不連續(xù)問題提供了一種新的思路.

圖34 單軸壓縮混凝土柱的近場動力學(xué)模擬[101]Fig.34 Peridynamics simulation of uniaxial compression concrete columns[101]

3.1.3 近場動力學(xué)方法的優(yōu)勢及存在的主要問題

近場動力學(xué)方法克服了基于局部作用思想建模和求解空間微分方程在處理非連續(xù)問題時遇到的奇異性困難,突破了經(jīng)典分子動力學(xué)方法在計算尺度上的局限,在宏觀和微觀非連續(xù)力學(xué)問題分析中具有明顯的優(yōu)勢[91].僅經(jīng)過十幾年的發(fā)展,近場動力學(xué)在模擬材料非連續(xù)破壞方面已展現(xiàn)出明顯的優(yōu)越性.因此,隨著近場動力學(xué)理論與應(yīng)用研究的不斷深入將促進巖石力學(xué)的進一步發(fā)展.

目前,近場動力學(xué)模型中的“邊界效應(yīng)”問題是需要解決的一個重要問題[122].對近場動力學(xué)模型中自由邊界附近的物質(zhì)點進行積分時,會出現(xiàn)鄰域不完整和受力不平衡的現(xiàn)象,這導(dǎo)致了“邊界效應(yīng)”問題的產(chǎn)生.無外力作用下發(fā)生自由膨脹時,由于積分域不完整導(dǎo)致模型自由邊界附近物質(zhì)點的鍵力分布不對稱,導(dǎo)致物質(zhì)點的受到非平衡力,降低了計算精度.鍵型近場動力學(xué)模型中的微模量c被視為常數(shù)也會引起“邊界效應(yīng)”問題.目前,主要采用以下辦法處理近場動力學(xué)模型中“邊界效應(yīng)”問題: (1)對自由邊界附近物質(zhì)點進行體積修正或虛構(gòu)邊界層;(2)進行力密度修正;(3)能量密度修正[91].然而,上述方法僅僅削弱了“邊界效應(yīng)”的影響,而沒有解決邊界問題的根本.

近場動力學(xué)方法的計算效率相對傳統(tǒng)數(shù)值計算方法較低,主要原因在于其在鄰域半徑范圍內(nèi)進行積分,且求解過程以顯式計算為主.因此,計算時間步需要足夠小才能獲得穩(wěn)定解,這將導(dǎo)致計算量增大.雖然近場動力學(xué)方法在解決非連續(xù)問題上具有優(yōu)勢,但其不完善性也限制了其在大規(guī)模問題上的應(yīng)用.為了提高計算效率,目前研究者們將近場動力學(xué)方法與傳統(tǒng)數(shù)值方法相結(jié)合.同時,隨著計算機技術(shù)的快速發(fā)展和近場動力學(xué)理論的不斷完善,計算效率問題正在逐步得到解決.

3.2 其他方法

除了以上論述的數(shù)值方法外,還有一些方法本文未詳細介紹,如Harmon 等[123]提出的水平集切割法,武漢大學(xué)王橋團隊[124-125]和周偉團隊[126]發(fā)展的相場理論方法,劉彪等[68]采用的邊界元方法,以及物質(zhì)點法等.這些方法都在顆粒破碎的研究過程發(fā)揮了重要的推動作用,為今后發(fā)展更為可靠、高效的數(shù)值方法奠定了基礎(chǔ).

4 總結(jié)與展望

4.1 總結(jié)

顆粒破碎是顆粒材料研究中的關(guān)鍵內(nèi)容.目前,顆粒破碎的研究已經(jīng)取得了一定的成果,并在工程領(lǐng)域得到廣泛的應(yīng)用,解決了一定的實際問題.然而,現(xiàn)有的方法大多是從現(xiàn)象出發(fā),缺乏對顆粒破碎機理的明確研究,因此,截止目前,對顆粒材料的破碎機理認識仍不清晰.此外,現(xiàn)有的方法一般預(yù)設(shè)了破壞路徑,這與真實的破碎過程存在一定的差異;有的方法僅可以破碎一次,破碎后的顆粒無法實現(xiàn)二次破碎,這與實際也是不符的;若不提前預(yù)設(shè)破壞路徑,且顆?？梢园l(fā)生多次破碎,那么計算效率將會非常低.面對以上問題,亟待提出一種從機理出發(fā)不需要預(yù)設(shè)破壞路徑且可以發(fā)生多次破碎的研究顆粒破碎的高效數(shù)值方法.

4.2 展望

針對顆粒破碎的數(shù)值研究,目前仍有一些問題存在,這將是今后研究的主要內(nèi)容和方向.

首先,為了探究顆粒破碎的內(nèi)在機理,則需要從破碎的本質(zhì)出發(fā),考慮顆粒真實破碎的內(nèi)在原因.從以上內(nèi)容不難看出,目前有關(guān)顆粒破碎的機理研究內(nèi)容較少,為了揭示顆粒破碎的內(nèi)在機理,必須將機理研究作為首要研究內(nèi)容.

其次,當從顆粒破碎的現(xiàn)象描述出發(fā),則需要考慮顆粒破碎程度的度量指標問題.對顆粒的破碎程度以及破碎后顆粒進行描述,這就需要有合理的描述指標來表示顆粒的破碎程度以及破碎后顆粒的相關(guān)性質(zhì).

然后,工程中顆粒材料破碎的計算效率問題: 從以上內(nèi)容不難看出,目前的研究方法均存在計算效率問題,尤其是面向工程實際時,需要考慮的顆粒數(shù)量往往會達到百萬甚至千萬級,面對如此巨大數(shù)量的顆粒,如何保證計算效率是今后工作的重點,也是將顆粒材料破碎應(yīng)用于工程的必要條件.

最后,顆粒破碎有效模擬的計算軟件問題: 面對機理研究以及工程應(yīng)用,如何更有效地研究顆粒破碎過程,研發(fā)有效的計算軟件將是解決這一問題最為直接和有效的方法.因此,研發(fā)具有自主知識產(chǎn)權(quán)的模擬顆粒破碎的計算軟件,也將是今后工作的重點,當然,這也是難點.