賈明坤， 王 偉， 張 斌， 許文祥

(河海大學(xué) 力學(xué)與材料學(xué)院 固體力學(xué)研究所，南京 211100)

1 引 言

顆粒材料廣泛存在于自然界和工程應(yīng)用中，具有復(fù)雜的物理力學(xué)特性。顆粒體系的細(xì)觀組構(gòu)性能是影響其宏觀物理力學(xué)性能的決定性因素[1-5]，已經(jīng)有大量工作對顆粒體系的細(xì)觀組構(gòu)性能進行了研究。在試驗研究方面，學(xué)者們利用傳統(tǒng)物理試驗獲得顆粒材料的宏觀性能，同時結(jié)合CT成像等技術(shù)[6,7]探究顆粒體系內(nèi)部的細(xì)觀組構(gòu)信息，但試驗過程較復(fù)雜,且價格昂貴。近年來，隨著計算機科學(xué)技術(shù)的發(fā)展，數(shù)值模擬成為研究顆粒堆積行為與細(xì)觀組構(gòu)性能的重要手段，涌現(xiàn)出一大批優(yōu)秀的研究成果。Zhao等[8,9]采用松弛算法，對四面體和球柱體等形狀顆粒的隨機緊密堆積進行了系統(tǒng)研究，并給出了堆積組構(gòu)性能(堆積分?jǐn)?shù)和配位數(shù)等)與形狀參數(shù)的關(guān)聯(lián)性。Jiao等[10]提出了自適應(yīng)收縮元胞ASC(Adaptive Shrinking Cell)算法，獲得了柏拉圖粒子的最大隨機阻塞堆積結(jié)構(gòu)(maximally random jammed (MRJ) packings)，并討論了柏拉圖粒子的形狀對其MRJ堆積分?jǐn)?shù)的影響。趙仕威[2]開發(fā)非球形顆粒離散元算法，深入研究了顆粒形狀對四面體和超橢球顆粒在外力作用下的隨機堆積組構(gòu)性能的影響。袁野[11]采用多種填充算法對超橢球和球多面體顆粒(內(nèi)核為正四面體、立方體和正八面體)的無序阻塞堆積進行了研究，討論了堆積性質(zhì)與顆粒形狀的復(fù)雜關(guān)聯(lián)?？梢?，顆粒形狀是影響顆粒體系組構(gòu)性能的最重要因素之一，數(shù)值模擬為探究顆粒形狀對顆粒體系組構(gòu)性能的影響提供了強有力的技術(shù)手段。然而，目前已有成果主要集中在對規(guī)則形狀顆粒(球、橢球、超橢球和正則多面體等)體系組構(gòu)性能的研究，對于復(fù)雜凸多面體顆粒體系組構(gòu)性能及其形狀依賴性還未有報道。

數(shù)值研究顆粒體系組構(gòu)性能的前提是構(gòu)造顆粒的堆積結(jié)構(gòu)，其核心問題包括兩個,即建立單個顆粒幾何模型并開發(fā)相應(yīng)的接觸檢測算法以及開發(fā)高效的顆粒堆積算法。一方面，建立單個復(fù)雜非球形顆粒模型，通常有解析法和組合球法兩種。解析法能夠精確構(gòu)造顆粒模型，但需要開發(fā)復(fù)雜的接觸算法[12-14]實現(xiàn)顆粒之間的接觸檢測；組合球法[15]是將顆粒轉(zhuǎn)化為球體的集合，通過球體的重疊判斷實現(xiàn)顆粒的接觸檢測。另一方面，學(xué)者們也已經(jīng)開發(fā)了很多顆粒堆積算法，包括隨機有序堆積算法[16]、松弛算法[8,9,15]和離散元算法[17,18]。其中，隨機有序堆積算法是一種靜態(tài)堆積算法，算法流程雖然簡單，但獲得的顆粒堆積結(jié)構(gòu)的堆積分?jǐn)?shù)較低，無法達到實際顆粒材料的堆積分?jǐn)?shù)；離散元算法基于牛頓第二定律模擬顆粒的真實受力情況和運動過程，但獲得的堆積結(jié)構(gòu)與接觸力模型相關(guān)，模擬過程較為復(fù)雜；而松弛算法簡單高效，且能夠獲得很高的堆積分?jǐn)?shù)，廣泛應(yīng)用于顆粒堆積結(jié)構(gòu)的模擬。

本文首先基于旋轉(zhuǎn)橢球面黃金螺旋網(wǎng)格構(gòu)造了一組復(fù)雜凸多面體顆粒模型(Polyκ - ng s)；然后采用松弛算法，并結(jié)合GJK算法和擴張多面體算法EPA(Expanding Polytope Algorithm)[20]，生成Polyκ - ng s多面體顆粒的隨機緊密堆積結(jié)構(gòu)；最后定量表征Polyκ - ng s多面體顆粒堆積體系的組構(gòu)性能，并探究形狀參數(shù)(長徑比κ和頂點數(shù)量ng s)的影響。

2 復(fù)雜凸多面體顆粒幾何模型

為了研究形狀更具一般性特征的凸多面體顆粒堆積體系的組構(gòu)性能，本文基于旋轉(zhuǎn)橢球面上的黃金螺旋網(wǎng)格構(gòu)造復(fù)雜凸多面體顆粒幾何模型，具體步驟包括，(1) 在單位球面上生成黃金螺旋網(wǎng)格；(2) 將單位球面上的黃金螺旋網(wǎng)格映射到旋轉(zhuǎn)橢球面上，生成新的復(fù)雜凸多面體。具體細(xì)節(jié)如下。

(1) 在單位球面上生成黃金螺旋網(wǎng)格。首先，在單位球面上進行黃金螺旋網(wǎng)格點采樣。黃金螺旋網(wǎng)格點采樣能夠在球體表面獲取均勻分布的采樣點，且采樣的數(shù)量可以是任意正整數(shù)。假設(shè)采樣點的數(shù)量是m，則在球坐標(biāo)系下第i個采樣點的天頂角和方位角坐標(biāo)公式[19]為

(1)

圖1 單位球面黃金螺旋網(wǎng)格

(2) 將單位球面上的黃金螺旋網(wǎng)格映射到旋轉(zhuǎn)橢球面上，生成新的復(fù)雜凸多面體。首先將獲得的黃金螺旋網(wǎng)格采樣點的天頂角和方位角坐標(biāo)代入旋轉(zhuǎn)橢球面方程，則可獲得旋轉(zhuǎn)橢球面上的黃金螺旋網(wǎng)格點，其直角坐標(biāo)公式為

(2)

式中κ為旋轉(zhuǎn)橢球的長徑比，ng s為黃金螺旋網(wǎng)格點數(shù)量；然后，保持單位球面黃金螺旋網(wǎng)格節(jié)點之間的幾何拓?fù)潢P(guān)系，則可以構(gòu)造出相應(yīng)的復(fù)雜凸多面體顆粒模型。本文選取的旋轉(zhuǎn)橢球面長徑比κ=0.4,0.6,0.8,1.0,1.5,2.0,3.0和4.0，選取的黃金螺旋網(wǎng)格點數(shù)量ng s=4,6,8,10，構(gòu)造出一組形狀按一定規(guī)律演變的復(fù)雜凸多面體顆粒，標(biāo)記為Polyκ - ng s。如圖2所示為κ=0.4,1.0,4.0的三組凸多面體顆粒?？梢?，基于黃金螺旋網(wǎng)格采樣點構(gòu)造的多面體，頂點數(shù)量可控，整體能夠維持對應(yīng)旋轉(zhuǎn)橢球體的形狀特征。

圖2 復(fù)雜凸多面體顆粒模型Polyκ - ng s(κ =0.4,1.0,4.0)

可以看出，當(dāng)旋轉(zhuǎn)橢球面長徑比κ相等時，黃金螺旋網(wǎng)格點數(shù)量ng s的增大使得凸多面體顆粒形狀向?qū)?yīng)κ的旋轉(zhuǎn)橢球演變；當(dāng)黃金螺旋網(wǎng)格點數(shù)量ng s相等時，旋轉(zhuǎn)橢球面長徑比κ的變化使得凸多面體沿著旋轉(zhuǎn)橢球半軸a方向進行伸縮。同時，本文研究了Polyκ - ng s多面體顆粒的體積Vp和表面積Sp隨形狀參數(shù)的變化，如圖3所示，在相同的ng s下，Polyκ - ng s多面體的體積和表面積隨著κ的增大線性增大，這是因為κ的增大使得顆粒沿橢球顆粒長軸方向線性拉伸；在相同的κ下，Polyκ - ng s多面體的體積和表面積隨著ng s的增大而增大，逐漸接近相應(yīng)的橢球顆粒。另外需要說明的是，本文用于表征顆粒形狀的參數(shù)，即長徑比κ和頂點數(shù)量ng s無法直接用于表征工程實際中的一般顆粒，但可以找到與κ和ng s相對應(yīng)的形狀參數(shù)來表征一般顆粒的形態(tài)，如包圍盒尺寸和圓度等[21]。

圖3 形狀參數(shù)對Polyκ - ng s多面體顆粒體積Vp和表面積Sp 的影響

3 復(fù)雜凸多面體顆粒隨機緊密堆積

3.1 凸多面體顆粒隨機堆積的松弛算法

采用松弛算法模擬周期性邊界條件下復(fù)雜凸多面體顆粒的隨機緊密堆積，具體實現(xiàn)細(xì)節(jié)如下。

(3)

式中Vi為第i個顆粒的體積。

(2) 松弛機制。首先，對堆積結(jié)構(gòu)中的顆粒進行重疊檢測，并求解每一對重疊顆粒的接觸信息，包括接觸法向、接觸點和疊合量。本文采用GJK算法對凸多面體顆粒進行重疊檢測，并結(jié)合擴張多面體算法EPA(Expanding Polytope Algorithm)求解重疊顆粒對的接觸法向、接觸點和疊合量，算法細(xì)節(jié)可參見文獻[20]，本文不再贅述。如圖4所示，一個四面體和一個八面體重疊，采用GJK和EPA算法可求解出接觸法向n、接觸點C和疊合量δ。然后，根據(jù)接觸信息給顆粒施加平動和轉(zhuǎn)動位移矢量以減小顆粒之間的重疊。圖4的八面體獲得的平動位移矢量T和繞質(zhì)心B的轉(zhuǎn)動位移矢量R為

圖4 重疊凸多面體的接觸信息與松弛位移矢量

(4)

式中α為平動系數(shù)，β為轉(zhuǎn)動系數(shù)。當(dāng)一個顆粒和多個顆粒存在重疊時，需要將所有位移矢量疊加，獲得顆?？偟钠絼雍娃D(zhuǎn)動位移矢量。對所有顆粒重疊檢測、接觸求解并施加位移矢量為一個松弛循環(huán)，重復(fù)松弛循環(huán)可以不斷減小整個堆積結(jié)構(gòu)中所有接觸顆粒之間的重疊。

(3) 松弛過程的停止和堆積區(qū)域的擴張。每個松弛循環(huán)后，檢測所有接觸顆粒的疊合量是否小于預(yù)設(shè)值δmax，若滿足，則停止循環(huán)，獲得的結(jié)構(gòu)即為一個滿足周期邊界條件的凸多面體顆粒堆積結(jié)構(gòu)；若不滿足，則繼續(xù)松弛。重復(fù)T次松弛循環(huán)記為一個松弛迭代周期，每個迭代周期后，擴大堆積區(qū)域，即增大立方體容器邊長，增長量為

(5)

mor=δ/[0.5(Re q 1+Re q 2)]

(6)

式中δ為兩重疊顆粒的疊合量，Re q 1和Re q 2分別為兩顆粒的等效半徑，等效半徑Re q和等效直徑De q是用于度量非球形顆粒尺寸的參數(shù)，

(7)

式中VP為顆粒體積。

需要強調(diào)的是，在松弛算法中，只考慮了顆粒之間的法向接觸和排斥作用，故只適用于無摩擦顆粒系統(tǒng)。

3.2 松弛算法參數(shù)取值及可靠性驗證

圖5 松弛過程中堆積分?jǐn)?shù)p隨顆粒之間最大疊合量δm的變化(Poly4 - 10)

最后，本文保持其他參數(shù)不變，δmax=10-3De q，T=100，γ=0.01，檢測平動系數(shù)α和轉(zhuǎn)動系數(shù)β對Poly4 - 10多面體堆積分?jǐn)?shù)的影響，結(jié)果如圖6所示?？梢钥闯?，對于某一固定的平動系數(shù)α，松弛堆積分?jǐn)?shù)隨著轉(zhuǎn)動系數(shù)β的增大，先增大后減??；而對于不同的平動系數(shù)α，松弛獲得的最高堆積分?jǐn)?shù)相差不大。對于Poly4 - 10多面體，當(dāng)平動系數(shù)α取0.75，轉(zhuǎn)動系數(shù)β取1.00時，可以獲得最高堆積分?jǐn)?shù)為0.678，堆積結(jié)構(gòu)如圖7所示，圖7還展示了該結(jié)構(gòu)在z=0平面的截面示意圖，可見，顆粒之間不存在不可接受的重疊，且結(jié)構(gòu)滿足周期邊界條件，證明3.1節(jié)描述的松弛算法能夠獲得復(fù)雜凸多面體顆粒的緊密堆積結(jié)構(gòu)。

圖6 平動系數(shù)α和轉(zhuǎn)動系數(shù)β對松弛堆積分?jǐn)?shù)p的影響(Poly4 - 10)

圖7 Poly4.0 - 10多面體最高堆積分?jǐn)?shù)對應(yīng)的堆積結(jié)構(gòu)及截面

試算發(fā)現(xiàn),平動和轉(zhuǎn)動系數(shù)對最終堆積分?jǐn)?shù)的影響最大，故本文保持其他參數(shù)不變，δmax=10-3De q，T=100，γ=0.01，改變平動和轉(zhuǎn)動系數(shù)，來探索顆粒能獲得的最高堆積分?jǐn)?shù)。

4 復(fù)雜凸多面體隨機緊密堆積的組構(gòu)性能

本文研究Polyκ - ng s多面體顆粒隨機緊密堆積的組構(gòu)性能及其形狀依賴性。根據(jù)3.2節(jié)，在松弛算法中，對于同一種形狀的凸多面體顆粒，不同的算法參數(shù)取值能夠獲得不同堆積分?jǐn)?shù)的結(jié)構(gòu)。而本節(jié)研究的復(fù)雜凸多面體隨機緊密堆積，不作特殊聲明，默認(rèn)指的是當(dāng)前松弛算法能獲得的Polyκ - ng s多面體顆粒最高堆積分?jǐn)?shù)對應(yīng)的緊密堆積結(jié)構(gòu)。首先獲取松弛算法下Polyκ - ng s多面體顆粒的隨機緊密堆積結(jié)構(gòu)；然后表征結(jié)構(gòu)中顆粒體系的組構(gòu)性能，并討論組構(gòu)性能與形狀參數(shù)之間的關(guān)聯(lián)性。

4.1 顆粒位置分布與方向分布

首先對每種Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)的顆粒位置分布進行檢測，在全局坐標(biāo)系x，y和z方向上分別作堆積結(jié)構(gòu)的截面，每個方向的截面數(shù)量取100，截面等間距，計算與截面相交的顆粒數(shù)量，并求平均值，將每個截面上與截面相交的顆粒數(shù)量與平均值之比記為該截面上顆粒的數(shù)量密度比。圖8展示了圖7所示的Poly1.0 - 4多面體隨機緊密堆積結(jié)構(gòu)的截面數(shù)量密度比，可以看出三個坐標(biāo)軸方向的數(shù)量密度比均在0.85～1.2之間波動，其他形狀的Polyκ - ng s多面體的隨機緊密堆積結(jié)構(gòu)的截面數(shù)量密度比均表現(xiàn)出和圖8相同的規(guī)律?？梢?，對松弛算法獲得的凸多面體顆粒隨機緊密堆積結(jié)構(gòu)，3個坐標(biāo)軸方向上的顆粒位置分布是均勻的。

圖8 圖7所示堆積結(jié)構(gòu)沿坐標(biāo)軸方向的截面顆粒數(shù)量密度比

另外，本文計算了每種Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)的徑向分布函數(shù)g(r)，如圖9所示。

圖9 基于松弛算法的Polyκ - ng s顆粒隨機緊密堆積結(jié)構(gòu)的徑向分布函數(shù)

可以看出，當(dāng)ng s相等時，Polyκ - ng s多面體顆粒堆積結(jié)構(gòu)的徑向分布函數(shù)呈現(xiàn)相似的規(guī)律，先波動后衰減為1，當(dāng)κ越接近1.0時，波動越明顯，可見相應(yīng)的結(jié)構(gòu)在位置分布上存在長程順序；當(dāng)κ相等時，隨著ng s的增大，相應(yīng)堆積結(jié)構(gòu)的徑向分布函數(shù)的波動越明顯，可見當(dāng)顆粒形狀接近相應(yīng)的橢球時，位置長程有序性更加明顯。

然后對每種Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)的顆粒方向分布進行檢測。使用Polyκ - ng s多面體在局部坐標(biāo)系下的x軸正方向單位矢量表征顆粒方向。先根據(jù)結(jié)構(gòu)中每個顆粒的方向，計算其局部坐標(biāo)系下的矢量(1,0,0)在全局坐標(biāo)系下的矢量，然后將該矢量分別投影到xOy，xOz和yOz平面，并計算投影矢量方向的概率分布。圖10為不同形狀Polyκ - ng s多面體的方向矢量在xOy平面投影矢量方向的概率分布?？梢钥闯?，松弛算法獲得的隨機緊密堆積結(jié)構(gòu)中顆粒的方向分布不均勻，長徑比κ和頂點數(shù)量ng s都對顆粒方向分布有一定的影響。長徑比κ是主要影響因素，對于ng s相同的Polyκ - ng s多面體顆粒，κ越遠(yuǎn)離1，堆積結(jié)構(gòu)中的顆粒方向分布越不均勻。Polyκ - ng s多面體的方向矢量在xOz和yOz平面投影矢量方向的概率分布呈現(xiàn)的規(guī)律與圖10相似。

圖10 基于松弛算法的隨機緊密堆積結(jié)構(gòu)中Polyκ - ng s顆粒的方向矢量在x Oy平面投影的概率分布

4.2 堆積分?jǐn)?shù)

堆積分?jǐn)?shù)是顆粒材料最重要且使用最廣泛的組構(gòu)性能參數(shù)，指的是顆粒體積與堆積結(jié)構(gòu)體積之比。圖11展示了基于松弛算法獲得的Polyκ - ng s多面體顆粒的最高堆積分?jǐn)?shù)?？梢钥闯觯?dāng)ng s相等時，Polyκ - ng s多面體顆粒的最高堆積分?jǐn)?shù)隨著長徑比κ的增大先增大后減小，當(dāng)κ=1.0時達到最大；對于κ相等的Polyκ - ng s多面體顆粒，當(dāng)ng s=4時，顆粒最高堆積分?jǐn)?shù)最低，當(dāng)ng s≥6時，ng s對顆粒最高堆積分?jǐn)?shù)的影響不顯著；當(dāng)κ=1.0且ng s=6時，堆積分?jǐn)?shù)達到最大值0.720。

圖11 基于松弛算法獲得的Polyκ - ng s多面體顆粒的最高堆積分?jǐn)?shù)

4.3 配位數(shù)與顆粒接觸類型

圖12 基于松弛算法的隨機緊密堆積結(jié)構(gòu)的配位數(shù)概率分布

圖13 基于松弛算法的隨機緊密堆積結(jié)構(gòu)的平均配位數(shù)

除了配位數(shù)大小，對于多面體顆粒，顆粒接觸類型更能反應(yīng)顆粒堆積的緊密程度，本文對每種顆粒隨機緊密堆積結(jié)構(gòu)中的顆粒接觸類型數(shù)量進行統(tǒng)計。多面體的接觸類型可以劃分為6種，包括點-點(v - v)、點-邊(v - e)、邊-邊(e - e)、點-面(v - f)、邊-面(e - f)和面-面(f - f)，由于 v - v 和 v - e 接觸只出現(xiàn)在兩個顆粒只有一點接觸的情況，所以在堆積結(jié)構(gòu)中極少出現(xiàn)，圖14為Poly1.0 - 4多面體的 v - f,e - e,e - f和f - f 四種接觸類型。圖15展示了本文獲得的Polyκ - ng s多面體顆粒隨機緊密堆積結(jié)構(gòu)中各種接觸類型的接觸數(shù)?？梢钥闯?，雖然所有形狀顆粒堆積結(jié)構(gòu)的平均配位數(shù)都在7～8.5之間，但各接觸類型數(shù)量所占比例不同，f - f 接觸的占比通常能夠反映多面體顆粒堆積的緊密程度，f - f 接觸占比隨著κ的增大呈先增大后減小的趨勢，與圖11相應(yīng)的堆積分?jǐn)?shù)的規(guī)律相對應(yīng)；當(dāng)ng s和κ同時較大時，f - f 接觸占比非常高，但此時堆積分?jǐn)?shù)并不高，原因是當(dāng)ng s較大時，顆粒三角面數(shù)量較多但面積較小，即使較多的 f - f 接觸數(shù)也不一定有較大的有效接觸面積。

圖14 接觸類型(Poly1 - 4多面體)

圖15 基于松弛算法的隨機緊密堆積結(jié)構(gòu)的接觸類型數(shù)量

5 結(jié) 論

本文基于旋轉(zhuǎn)橢球面黃金螺旋網(wǎng)格構(gòu)造Polyκ - ng s多面體顆粒模型，然后通過松弛算法獲得隨機緊密堆積結(jié)構(gòu)，最后表征Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)的組構(gòu)性能，研究了形狀參數(shù)對組構(gòu)性能的影響，得出了以下主要結(jié)論。

(1) 結(jié)合GJK和EPA算法進行凸多面體接觸求解，松弛算法能夠用于獲得凸多面體緊密堆積結(jié)構(gòu)，平移系數(shù)和旋轉(zhuǎn)系數(shù)對堆積分?jǐn)?shù)影響最大。

(2) Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)中，顆粒位置分布均勻，當(dāng)長徑比κ接近1，頂點數(shù)量ng s增大時，堆積結(jié)構(gòu)表現(xiàn)出更強的位置長程有序性；顆粒方向分布不均勻，方向分布的不均勻程度主要受長徑比κ影響。

(3) Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)的堆積分?jǐn)?shù)與長徑比κ和頂點數(shù)量ng s有關(guān)。當(dāng)ng s相等時，最高堆積分?jǐn)?shù)隨著長徑比κ的增大先增大后減小；對于κ相等的Polyκ - ng s多面體顆粒，當(dāng)ng s=4時，顆粒最高堆積分?jǐn)?shù)最低；當(dāng)κ=1.0且ng s=6時，堆積分?jǐn)?shù)達到最大值0.720。

(4) 基于松弛算法獲得的Polyκ - ng s多面體隨機緊密堆積結(jié)構(gòu)的配位數(shù)分布符合高斯分布，頂點數(shù)量ng s影響配位數(shù)分布的峰值。不同形狀的Polyκ - ng s多面體的隨機緊密堆積結(jié)構(gòu)中的接觸類型分布不同，f - f 接觸占比隨著κ的增大呈先增大后減小的趨勢。