顆粒材料破碎行為數(shù)值分析方法研究綜述1)

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

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

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

引言

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

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

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

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

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

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

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

1.1 黏結(jié)?破碎法

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

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

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

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

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

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

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

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

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

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

圖3 單顆粒、三顆粒柱的SMT 和DEM 圖像以及不同載荷水平下實(shí)驗(yàn)室規(guī)模的壓縮試驗(yàn)[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)場(chǎng)[27]Fig.6 Bonding element and particle displacement and rotation field under 20% axial strain[27]

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

圖7 橢球體顆粒整體破碎過(guò)程演化(俯視圖,顆粒間線代表接觸黏結(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),在實(shí)現(xiàn)顆粒破碎的模擬中可以對(duì)顆粒強(qiáng)度進(jìn)行表征,因此,在分析此類(lèi)工程問(wèn)題時(shí)可以采用該方法.

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

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

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

圖9 錐體結(jié)構(gòu)的破冰過(guò)程[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 通過(guò)黏結(jié)單元法構(gòu)造了粗粒土模型,并討論了直剪試驗(yàn)中基礎(chǔ)顆粒間的黏結(jié)力、黏結(jié)顆粒的孔隙率和粗粒土試樣的孔隙率等因素對(duì)顆粒材料體應(yīng)變和剪切強(qiáng)度的影響.

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

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

1.2 碎片替代法

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

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

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

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

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

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

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

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

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

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

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

圖12 碎片替代法研究壓縮作用下顆粒破碎的演化過(guò)程[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]

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

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

1.2.2 碎片替代法的破碎準(zhǔn)則研究

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 比例邊界有限元法

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

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

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

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

2.1.2 比例邊界有限元法的破碎準(zhǔn)則

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

圖16 樣本點(diǎn)的分布Fig.16 Distribution of sample points

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

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

圖17 由Hock-Brown 準(zhǔn)則得到的破壞點(diǎn)Fig.17 Failure points obtained from the Hock-Brown criterion

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 內(nèi)聚力模型

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

FDEM 方法中,有限元離散化以后,將黏性界面單元插入到與每個(gè)粒子相關(guān)的有離散單元中,以考慮粒子的破碎過(guò)程.本節(jié)將重點(diǎn)介紹黏性界面的內(nèi)聚力模型(cohesive zone model,CZM)[9,80-82].內(nèi)聚力模型是1960 年Dugdale 等提出的[83-84],假設(shè)裂紋的尖端是由兩個(gè)裂紋界面組成的很小的內(nèi)聚區(qū),該內(nèi)聚區(qū)的本構(gòu)關(guān)系就是界面上作用牽引力T與兩裂紋面間相對(duì)位移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)聚力模型從未完全承載的點(diǎn)A開(kāi)始,T隨著U的增加而增加,隨之達(dá)到應(yīng)力最大值Tmax的點(diǎn)C,此時(shí)材料的應(yīng)力承載達(dá)到最大值,材料在該點(diǎn)處開(kāi)始出現(xiàn)初始損傷.隨著界面位移的繼續(xù)增大,應(yīng)力開(kāi)始下降,該階段為材料的損傷擴(kuò)展階段,點(diǎn)E為材料裂紋界面完全分離點(diǎn),此時(shí)材料的承載力降為0[85].內(nèi)聚區(qū)應(yīng)力的變化受內(nèi)聚力模型和裂紋界面位移的影響,因此,針對(duì)不同的材料,應(yīng)該選擇不同的內(nèi)聚力模型,同時(shí)通過(guò)選取適當(dāng)?shù)膮?shù),從而反映界面的強(qiáng)度、韌度等力學(xué)性能.

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

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

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

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

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

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

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

為了模擬巖石破裂過(guò)程,廣東工業(yè)大學(xué)高偉團(tuán)隊(duì)提出一種新的內(nèi)聚模型本構(gòu).在該本構(gòu)中,采用Mohr-Coulomb 準(zhǔn)則作為裂紋萌生準(zhǔn)則,建立壓縮和拉伸狀態(tài)下的牽引?分離模型,描述了巖石強(qiáng)壓力依賴性的力學(xué)特性.在粘結(jié)區(qū)完全破壞后,采用基于庫(kù)侖公式的摩擦模型來(lái)確定裂紋自由面之間的接觸相互作用.采用所提出的零厚度黏合單元本構(gòu),模擬了巴西盤(pán)、單軸壓縮兩個(gè)巖石破裂實(shí)例.比較仿真結(jié)果與實(shí)驗(yàn)結(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)聚力單元模擬巴西盤(pán)破裂和單軸壓縮巖石破裂實(shí)例[9]Fig.27 Simulation of brazilian plate rupture and example of rock fracture using cohesive elements[9]

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

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

3 近場(chǎng)動(dòng)力學(xué)方法及其他方法

3.1 近場(chǎng)動(dòng)力學(xué)方法

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

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

近場(chǎng)動(dòng)力學(xué)方法將固體材料離散為一系列具有體積和質(zhì)量的物質(zhì)點(diǎn),每個(gè)物質(zhì)點(diǎn)與其周?chē)欢▍^(qū)域內(nèi)所有其他物質(zhì)點(diǎn)之間存在相互作用力.如圖28所示.

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

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

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

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

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

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

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

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

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

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

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

圖32 撞擊速度為61 m/s 的破碎圖,以上分別是玻璃正面、底面和實(shí)驗(yàn)結(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]提出了一種用于顆粒狀和巖石狀材料斷裂的雙層近場(chǎng)動(dòng)力學(xué)(DH-PD)方法.與離散裂紋方法相比,DH-PD 方法不需要任何裂紋表面的表示形式,也不需要任何標(biāo)準(zhǔn)來(lái)處理復(fù)雜的裂紋形態(tài),裂紋路徑是模擬的自然結(jié)果.在該方法中作者提出了一種新的模擬壓縮裂縫接觸和約束侵徹條件的懲罰方法.該方法被應(yīng)用于地質(zhì)力學(xué)的一些基準(zhǔn)問(wèn)題,圖33 為采用雙層近場(chǎng)動(dòng)力學(xué)方法研究方形板裂紋擴(kuò)展的數(shù)值模擬結(jié)果.

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

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

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

3.1.3 近場(chǎng)動(dòng)力學(xué)方法的優(yōu)勢(shì)及存在的主要問(wèn)題

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

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

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

3.2 其他方法

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

4 總結(jié)與展望

4.1 總結(jié)

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

4.2 展望

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

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

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

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

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