考慮堆石料顆粒形狀的大三軸數(shù)值模型構(gòu)造

郜家奇，遲世春

(1.大連理工大學(xué) 海岸與近海工程國家重點實驗室， 遼寧 大連 116024；2.大連理工大學(xué) 建設(shè)工程學(xué)部水利工程學(xué)院 工程抗震研究所， 遼寧 大連 116024)

堆石料作為高土石壩的主要建筑材料，其力學(xué)特性引起了研究者的廣泛關(guān)注，其中大型三軸剪切試驗是主要的研究手段。但是，室內(nèi)試驗結(jié)果受多種因素影響，可重復(fù)性差。隨著計算機技術(shù)的發(fā)展，離散元方法(DEM)被廣泛用于研究不同因素對堆石料力學(xué)特性的影響[1-3]。但是如何生成具有真實顆粒形狀和處于壓實狀態(tài)的堆石顆粒成為目前亟待解決的問題[4-6]。

采用DEM進行堆石料數(shù)值模擬，歷經(jīng)了由簡單到復(fù)雜、二維到三維的過程。對于堆石料顆粒形狀，早期一般使用單一圓盤或圓球，后發(fā)展為使用多個圓盤、圓球或者多邊形等來模擬。比如在橢圓形上隨機取點構(gòu)造不規(guī)則多邊形來描述顆粒形狀[7-9]；Azéma等[10]選用規(guī)則的五邊形研究顆粒形狀的影響；任樹林等[11]以類正方形、類三角形和類長方形三種形狀表征異形顆粒；邵磊等[12]以正四面體為核進行晶胞繁衍構(gòu)造不同形狀的堆石顆粒；Cheng等[13]提出在最密六方排列的近球體顆粒的基礎(chǔ)上通過移除一定數(shù)量的基本球來模擬真實的顆粒特征。這些模擬方法雖然在一定程度上可以描述顆粒的形狀特征，但是與真實的堆石顆粒形狀差距較大，也并未與扁平度、棱角度等形狀描述指標(biāo)相聯(lián)系。針對這一問題，Mollon等[14]采用傅里葉系數(shù)與平均粒徑之比定義了實數(shù)傅里葉描述符，采用傅里葉描述符的線性組合方程，用3個主傅里葉描述符表征顆粒形狀。方浩等[15]在Mollon等的基礎(chǔ)上深入探究了實數(shù)傅里葉描述符與平面顆粒形狀指標(biāo)參數(shù)之間的相關(guān)性，并證明了這種方法的重構(gòu)顆粒與實際顆粒形狀間具有較高的相似度。

構(gòu)造顆粒集合體，金磊等[16]基于擬振動法和分層復(fù)制法提出了適用于土石混合體的大三軸試樣建模方法，試樣相對密實度0.5；邵磊等[12]應(yīng)用落雨法和分層振搗法生成顆粒集合體，為提高計算效率，剔除了小粒徑組，放大了顆粒粒徑；史旦達等[17]采用重力沉積法建模三軸試樣，建模時步達300萬步；為了快速生成指定體積的均質(zhì)性良好的顆粒簇，Ciantia等[18]提出了模塊組裝法；Jiang等[19]在總結(jié)前人定點成樣法[20]、等向壓密法[21]和粒徑放大法[22]的基礎(chǔ)上提出了分層欠壓法制備均勻性較好的試樣，平均孔隙率0.30，適用于顆粒數(shù)目較少的模型；仝帆等[23]采用大三軸數(shù)值模型模擬了新疆某工程粗粒料顆粒破碎，顆粒最小粒徑為5 mm。上述建模方法雖然在構(gòu)建大三軸數(shù)值模型時具有一定的適用性，但是往往并不能滿足實驗室級配和孔隙率要求。

本文基于傅里葉描述符生成堆石單顆粒模型，結(jié)合堆石料形狀參數(shù)，對單顆粒模型進行優(yōu)化消減，最終建立了適用于堆石料大三軸模型數(shù)值計算的三維離散元顆粒模型。模型構(gòu)建采用Clump-ball混合的落雨式分層過壓法，重構(gòu)了粒徑跨度2 mm～60 mm、孔隙率0.22的大三軸數(shù)值模型。探討了如何使數(shù)值試樣既滿足實驗室級配粒徑跨度大的要求，又接近室內(nèi)試驗試樣的較小孔隙率。

1 不規(guī)則形狀顆粒模型

1.1 顆粒形狀描述指標(biāo)

顆粒形狀描述指標(biāo)按維度可分為三類：1D、2D、3D。1D形狀描述指標(biāo)包括長度、寬度、表面積、體積、直徑、周長。2D指標(biāo)一般包含縱橫比、凹凸度、球度、圓度、整體規(guī)則性、棱角性等。3D指標(biāo)包含：真球度、扁平度、三維球度、有效球度等。研究者一般選擇扁平度、伸長率、縱橫比、圓度、球度和規(guī)則度來描述顆粒形狀。圖1為顆粒的尺寸、最大內(nèi)切圓和最小外接圓示意圖。常用的形狀描述指標(biāo)定義為：

扁平度:e=S/I

(1)

伸長率：λ=I/L

(2)

縱橫比：γ=S/L

(3)

(4)

(5)

(6)

圖1 顆粒形狀示意圖

1.2 傅里葉函數(shù)法基本原理

構(gòu)造不規(guī)則形狀堆石顆粒的辦法主要包括：傅里葉函數(shù)法[14]、球諧函數(shù)法[24]以及POD正交分解法[25]。其中，以傅里葉函數(shù)法應(yīng)用范圍最廣，不僅可以模擬棱角度較低、光滑的卵石，而且可以模擬棱角度較高、粗糙的礫石。

如圖2所示，以二維為例，以形心O為圓心，將顆粒輪廓等角度θ劃分為N份，取為27=128，θ=2π/N。傅里葉函數(shù)法基本公式：

(7)

(8)

(9)

(10)

式中:r0表示顆粒平均半徑；An、Bn表示傅里葉系數(shù)；n稱為諧波序號；ri表示形心點O到輪廓點的距離;θi表示形心等分角。

圖2 傅里葉函數(shù)法基本原理

傅里葉描述符Dn：

(11)

Dn=2α·log2(n/3)+log2(D3)3

(12)

Dn=2β·log2(n/8)+log2(D8)8

(13)

式中:D1表示顆粒輪廓對中心點的偏移程度；D2表示顆粒伸長程度；D3～D7表示顆粒輪廓的不規(guī)則特性；D8～D64描述顆粒表面粗糙程度等細節(jié)特征。α和β分別取值-1.5和-2。僅需考慮D2、D3和D8的值，即可確定所有傅里葉描述符的值。

圖3為D2取值0.0、0.1、0.2，D3取值0.00、0.05、0.10以及D8取值0.000、0.015與0.030，在不同組合下的顆粒形狀。如圖3(a)所示，D2=D3=D8=0，模型為光滑球形。保持D3、D8值不變，隨著D2值的增大，顆粒形狀趨于扁平；保持D2、D8值不變，D3取值增大，此時顆粒不再光圓，形狀更有棱角；保持D2、D3值不變，D8取值增大，顆粒表面的粗糙度明顯增大。

圖3 傅里葉描述符對顆粒形狀的影響

分別統(tǒng)計D2、D3、D8參數(shù)值變化與扁平度、伸長率、縱橫比、圓度、球度之間的關(guān)系。如圖4(a)所示，隨著D2值的增大，顆粒扁平度、縱橫比和球度值都減小，表明顆粒形狀越來越扁平；圓度值增大。由圖4(b)可得D3與圓度正相關(guān)，與球度負(fù)相關(guān)。由圖4(c)可得D8對顆粒形狀影響較小，對圓度影響較大，成明顯的負(fù)相關(guān)關(guān)系。綜上，D2對顆粒形狀影響最大。

圖4 不同傅里葉描述符取值與扁平度、伸長率、縱橫比、圓度、球度關(guān)系

1.3 構(gòu)造顆粒形狀模板庫

通過改變D2、D3、D8的取值，可構(gòu)造出一系列不規(guī)則形狀顆粒。取D2=0.1，D3=0.05，D8=0.015時，顆粒扁平度為0.818、伸長率為0.767、縱橫比0.625、圓度0.316、球度0.714，處于堆石料實際顆粒形狀描述指標(biāo)的取值范圍[26]。傅里葉函數(shù)法構(gòu)造顆粒形狀，由批量離散點構(gòu)成。對于任意一個內(nèi)部節(jié)點，計算其到表面的距離，取最小值。遍歷所有內(nèi)部點，以內(nèi)部點為球心，以對應(yīng)的最小距離值為半徑生成Pebble單元，所有Pebble單元組合成一個Clump顆粒。如圖5為通過多個Pebble單元組合疊加而成的不規(guī)則形狀顆粒，提取顆粒基本信息，將數(shù)值模型導(dǎo)入PFC3D。

顆粒數(shù)量直接影響計算效率。如圖6(a)所示，一個Clump顆粒由550個Pebble單元構(gòu)成，預(yù)估大三軸模型單元數(shù)量高達200萬個，計算效率無法保證。按照保留半徑較大值的辦法，把組成一個顆粒的單元數(shù)分別消減到200個、50個、10個時，如圖6所示。消減后的模型能基本刻畫堆石顆粒的不規(guī)則形狀，將由10個單元構(gòu)造的顆粒模型定義為三維單顆粒模型。在構(gòu)造大三軸數(shù)值模型時，堆石顆粒粒徑不同，因此在將一系列不同形狀的顆粒存儲為形狀模板之前，考慮到傅里葉函數(shù)法的尺度不變性[27]，可以對已構(gòu)造好的顆粒進行各向等比例粒徑歸一化處理，以方便后續(xù)建模調(diào)用時控制粒徑。

圖5 不規(guī)則形狀的單顆粒模型

圖6 消減后的單顆粒模型

2 大三軸數(shù)值模型構(gòu)造

關(guān)于顆粒集合體的建模，學(xué)者們提出了純幾何法[14]、膨脹法[28-30]、分層欠壓法[19]、自重沉積法[31-32]、模塊組裝法[18]等。純幾何法、膨脹法適用于圓盤或球體顆粒；自重沉積法則適用范圍較廣，建模過程能保證顆粒位置的隨機性。自重沉積法構(gòu)造模型達不到指定密實度，且顆粒粒徑分布不均勻(上疏下密)。對自重沉積法進行優(yōu)化，本文采用了落雨式分層過壓法構(gòu)造大三軸模型。具體建模流程如圖7所示。

圖7 大三軸建模流程圖

2.1 生成顆粒

Cho等[5]認(rèn)為顆粒形狀越不規(guī)則，臨界摩擦角越大，越難達到較小孔隙率；并且由實驗室統(tǒng)計得，顆粒粒徑越小，顆粒形狀越近似球形。因此，本文采用Clump-ball混合法來表征顆粒，粒徑大于10 mm的顆粒用不規(guī)則形狀的Clump表示，粒徑小于10 mm的顆粒采用Ball表示。

大三軸模型半徑0.15 m，高0.60 m。若直接同時生成全部顆粒，顆粒體積偏差較大，制樣速度慢。本文采用分層法，將模型分為6層，每層高0.1 m。如圖8所示，在4倍高 (0.4 m)壓力室內(nèi)隨機生成指定體積的堆石顆粒，顆粒位置隨機分布。連接鍵賦值之后，顆?；ハ鄰楅_。但是有時顆粒夾在兩顆粒之間，無法彈開；或者如圖9所示兩顆粒相互嵌套，形成“特殊顆粒”。遍歷接觸信息，定位特殊顆粒，刪除后按照同粒組、等體積原則重構(gòu)顆粒，投放至壓力室內(nèi)，見圖10。

圖8 生成指定體積顆粒

圖9 “特殊顆?！笔疽?/p>

圖10 壓力室內(nèi)顆粒無重疊

2.2 分層建模

指定體積的顆粒生成之后，添加重力場，顆粒在重力作用下向下移動。為了增加顆粒移動速度，減少顆粒觸底反彈位移，設(shè)置重力場加速度為100 m/s2。大粒徑顆粒接觸之后形成骨架，小粒徑顆粒填充骨架之間的空隙，一定程度上小顆粒向下運動的空間較大，此時顆粒分布“下密上疏”。如圖11所示，受顆粒形狀和不同顆粒排布方式的影響，顆粒自由落體之后，上表面不平整，模型達不到指定孔隙率。實驗室制樣使用錘子擊實調(diào)整顆粒密實度，本文在模型內(nèi)添加剛性平板墻體(墻體設(shè)置大剛度)，自0.13 m處由上到下壓實。墻體向下移動時，恢復(fù)重力場加速度為10 m/s2。

圖11 顆粒自由堆積

由實驗室伺服圍壓300 kPa壓縮顆粒時，模型往往無法達到指定分層高度0.1 m。本文采用過(超)固結(jié)顆?；貜椃?，圍壓擴大1 000倍，顆粒之間產(chǎn)生較大重疊，將顆粒擠壓在0.095 m的高度范圍內(nèi)，釋放圍壓至300 kPa，顆?；貜椫疗胶鉅顟B(tài)。圖12(a)表示在正常圍壓作用下的顆粒位置，此時，Ball向下運動的通道被封閉；增大圍壓, Clump1和Clump2繼續(xù)向下運動，接觸重疊量增大，圖12(b)Ball向下運動通道被打開，Ball在擠壓和重力作用下繼續(xù)向下運動；恢復(fù)圍壓300 kPa，Clump1和Clump2在反作用力作用下“回彈”，如圖12(c)所示。結(jié)合自然界土的超固結(jié)現(xiàn)象，應(yīng)用到大三軸建模，一方面顆粒為純剛體，恢復(fù)圍壓至正常應(yīng)力水平進行數(shù)值實驗，顆粒不會因為超固結(jié)應(yīng)力歷史影響其力學(xué)性質(zhì)和變形特性，另一方面滿足大三軸模型孔隙率的要求。此時，壓力室體積7.069×10-3m3，Clump顆粒體積4.632×10-3m3，占比81.9%；Ball顆粒體積1.023×10-3m3。顆?？傮w積5.655×10-3m3，實際孔隙率0.22。

2.3 模型組裝

重復(fù)前述步驟，構(gòu)造6層分層模型。計算機并行建模，提高效率。模型組裝有兩種方法，一是分層復(fù)制法，另一種是Brick組裝法。分層復(fù)制法，遍歷顆粒獲取顆粒位置坐標(biāo)x,y,z及半徑信息r，將第二層模型顆粒高度(z)坐標(biāo)增加1倍層高(0.1 m)，將第三層模型顆粒高度(z)坐標(biāo)增加2倍層高(0.2 m)，直至第六層。這種方法需要重新賦值顆粒的密度、法向剛度、切向剛度、摩擦系數(shù)等參數(shù)。Brick組裝法，將建成的每一層模型定義為“塊”，編號為Brick1、Brick2、…、Brick6?！皦K”不僅保留了模型顆粒的物理性質(zhì)和接觸參數(shù)，而且避免了重新生成顆粒的繁瑣。在定義“塊”時，指定模型的基點。通過基點重新賦值新的坐標(biāo)完成組裝。模型組裝后，各層模型之間分層明顯，需讓模型內(nèi)所有顆粒在重力場內(nèi)循環(huán)至平衡狀態(tài)。如圖13(a)、圖13(b)、圖13(c)表示模型分層組裝流程，圖13(d)、圖13(e)、圖13(f)表示建模完成后的小粒徑組、大粒徑組及全級配的大三軸數(shù)值模型。

圖12 超固結(jié)回彈法示意圖

2.4 模型驗證

大三軸試驗是研究堆石料力學(xué)特性的主要方法之一，實驗室制樣通常經(jīng)過篩料、砸樣、伺服等過程。通過篩料將堆石料按照粒徑大小分組，并按照室內(nèi)大三軸試驗指定級配[33]進行稱重，備料。在數(shù)值試樣建模時，在生成堆石顆粒時即進行粒徑分組，建模完成后利用Clump.vol函數(shù)編程求出各粒徑組體積，繪制數(shù)值模型的級配曲線并與指定室內(nèi)試驗級配做對比。如圖14所示，數(shù)值模型級配曲線與室內(nèi)試驗級配曲線基本重合。數(shù)值模型的孔隙率為0.22，接近室內(nèi)試驗的孔隙率0.20。說明本文提出的方法可以生成滿足顆粒形狀和孔隙率要求的堆石料數(shù)值試樣。

圖13 模型組裝流程

圖14 顆粒級配曲線

3 結(jié) 論

本文采用不規(guī)則形狀堆石料生成方法和Clump-ball混合的落雨式分層過壓建模法，構(gòu)建了顆粒形狀、顆粒級配、試樣孔隙率和建模過程都更接近實驗室大三軸試樣的數(shù)值模型。本文主要結(jié)論如下：

(1) 傅里葉函數(shù)法構(gòu)造三維離散元顆粒模型。采用傅里葉描述符快速構(gòu)造不規(guī)則形狀顆粒，分析傅里葉描述符D2、D3、D8對扁平度、伸長率、縱橫比、圓度、球度等形狀描述指標(biāo)的影響。D2對顆粒形狀影響最大，且與圓度成正相關(guān)關(guān)系，與顆粒扁平度、縱橫比和球度負(fù)相關(guān)；D3主要影響顆粒棱角度，取值越大，顆粒形狀棱角越多；D8對顆粒形狀影響較小，主要與表面粗糙度成正相關(guān)關(guān)系。

(2) 采用Clump-ball混合的落雨式分層過壓法來構(gòu)造大三軸數(shù)值模型。Clump-ball混合法解決了全級配粒徑跨度大、顆粒數(shù)目多的問題；落雨式重力沉積法保證了顆粒位置分布的隨機性及模型組構(gòu)的各向異性；在分層法提高模型構(gòu)造效率的基礎(chǔ)上，過壓回彈法使模型顆粒更密實均勻。

(3) 數(shù)值模型完全按照堆石料實驗室級配，目標(biāo)孔隙率0.20，實際孔隙率0.22。

(4) 數(shù)值模型制樣過程更接近實驗室真實制樣過程，接下來可以進行大三軸數(shù)值剪切試驗，為進一步探討顆粒形狀對堆石料力學(xué)特性等的影響作鋪墊。