基于Voronoi圖的土工膜顆粒墊層隨機重構(gòu)模型

姜曉楨,田曉丹

(1. 南京水利科學(xué)研究院巖土工程研究所,江蘇 南京 210024;2.河海大學(xué)文天學(xué)院,安徽 馬鞍山 243031;3.河海大學(xué)水利水電學(xué)院,江蘇 南京 210098)

在土石壩壩面土工膜防滲結(jié)構(gòu)中,土工膜的墊層材料一般采用具有良好透水性的顆粒材料(碎石為主),以便及時排除膜下滲水[1]。但土工膜在水壓力作用下,容易被墊層中顆粒頂破或在顆粒間隙內(nèi)發(fā)生脹破。圍繞土工膜在工程中所存在脹破和頂破的破壞形式,目前使用的SL/T 235—1999《土工合成材料測試規(guī)程》提出了3種標準試驗來定性地衡量土工膜抵抗脹破和頂破的能力。由于是標準試驗,所以試驗中采用的試樣、夾具以及頂桿的尺寸是一定的,無法真實還原顆粒墊層上土工膜的真實受力條件。在頂破脹破試驗中,國外有學(xué)者采用呈等邊三角形排布的3個剛性錐形凸起來模擬凹凸不平的墊層表面,但是由于錐形凸起的間距和尺寸是固定的,所以試驗得出的結(jié)論無法推廣到實際應(yīng)用中[2-3]。事實上,顆粒墊層上土工膜頂破脹破現(xiàn)象除了與水壓力大小和土工膜自身物理力學(xué)特性有關(guān)外,墊層顆粒形狀、粒徑級配、孔隙率等是更為重要的影響因素,所以也有不少學(xué)者采用真實的顆粒材料進行試驗[4],但由于真實顆粒墊層中每一個顆粒的形狀、粒徑及位置是隨機的,同時其與顆粒墊層的孔隙率、粒徑級配等試驗條件密切相關(guān),所以采用真實顆粒進行試驗時,常常無法對試驗條件進行有效控制,同時受試驗組數(shù)和試樣大小的影響,想要獲得具有統(tǒng)計意義的試驗數(shù)據(jù),試驗的工作量較大。而數(shù)值模擬方法可以避免上述問題,特別是隨著近些年來離散元數(shù)值模擬方法的不斷進步,從單一的圓顆粒模型逐步發(fā)展出了可以考慮不同顆粒形狀的數(shù)值模型[5-9],給顆粒墊層的隨機重構(gòu)以及顆粒形狀參數(shù)的統(tǒng)計分析帶來了便利。本文提出了一種基于Voronoi圖的顆粒墊層隨機重構(gòu)模型,并與真實顆粒墊層進行比較分析,可為后續(xù)進一步研究土工膜在顆粒墊層上的頂破脹破概率與墊層顆粒形狀、粒徑級配和孔隙率之間的關(guān)系打下基礎(chǔ)。

1 隨機重構(gòu)模型的基本要求

a. 顆粒墊層隨機重構(gòu)模型中顆粒形狀、粒徑級配與孔隙率盡可能與真實顆粒相近。由于土工膜在顆粒墊層上一般是被不規(guī)則形狀的顆粒頂破,所以采用不規(guī)則的多邊形顆粒進行模擬比較合理。在級配良好、孔隙率較小的顆粒墊層中,小顆粒能夠有效地填充大顆粒之間的間隙,同時小顆粒之間的間隙尺寸也要比大顆粒之間的間隙尺寸小,增加了土工膜與墊層顆粒的接觸面積,這將降低土工膜在墊層顆粒上發(fā)生頂破脹破的可能性,所以粒徑級配與孔隙率也與土工膜的頂破脹破有關(guān)。

b. 由于真實的墊層顆粒之間不存在重疊,同時也不存在受力不平衡的懸浮顆粒,所以要求顆粒墊層的隨機重構(gòu)模型中,相鄰顆粒之間的重疊量要盡可能小,且每個顆粒均能通過相鄰顆粒之間的接觸關(guān)系實現(xiàn)受力平衡,需要盡可能提高算法精度從而減小顆粒間的重疊量,同時避免受力不平衡的懸浮顆粒的出現(xiàn)。

c. 模型中需要設(shè)置一定的隨機過程,隨機過程及其參數(shù)會對模型中的顆粒形狀、粒徑級配和孔隙率等產(chǎn)生影響,通過調(diào)整這些隨機過程的設(shè)置方法及其參數(shù),能夠達到控制顆粒形狀、粒徑級配和孔隙率等性狀的目的。

目前常用的離散元顆粒模型建模方法大致可分為3種:①下落法[10-12]。下落法按顆粒級配以及形狀的要求建立有限數(shù)量的一組顆粒,并通過離散元模型施加重力,模擬顆粒自由下落堆積的過程。該方法最主要的問題在于無法控制自由下落堆積而成的顆粒堆積體的孔隙率,且建模的效率較低。②膨脹法[13-15]。膨脹法先在指定區(qū)域內(nèi)生成小顆粒,然后增大顆粒半徑使得顆粒發(fā)生接觸,并通過接觸力使得顆粒發(fā)生運動,充滿整個模型空間。該方法雖然能夠較好地實現(xiàn)孔隙率及粒徑級配的模擬,但其粒徑的增大系數(shù)難以控制,常常造成顆粒間有很大的重疊接觸,顆粒間以及顆粒與模型邊界之間存在較大的接觸力,顆粒間的接觸力瞬間釋放會造成大量顆粒飛溢出邊界,導(dǎo)致結(jié)果不準確。③純幾何法[16-17]。純幾何方法中,顆粒體的堆積純粹基于幾何計算,不模擬顆粒的動力特性,所以建模效率較高,但其同樣無法控制孔隙率。

2 基于Voronoi圖的隨機重構(gòu)模型

Voronoi圖,又稱泰森多邊形或Dirichlet圖,它是一組連接兩鄰點直線的垂直平分線所組成的連續(xù)多邊形,是一種重要的圖形幾何結(jié)構(gòu)[18]。Voronoi圖在隨機模型中有著廣泛的應(yīng)用,其基本概念可形象地描述為:在區(qū)域內(nèi)有n個火源(形核點),這n個火源同時點燃,并以相同速度向所有方向蔓延,當兩個火源擴散后相遇時熄滅,那么燃燒熄滅處所形成的圖便是Voronoi圖[19]。

2.1 形核點的位置

要生成Voronoi圖,首先要確定形核點的位置。顯而易見,當形核點在區(qū)域內(nèi)規(guī)則排布時,生成的Voronoi圖及其多邊形也是規(guī)則的,如圖1(a)所示,所以形核點的位置應(yīng)具有一定的隨機性,使得生成的Voronoi圖中的多邊形形狀及大小產(chǎn)生隨機性。文中形核點位置產(chǎn)生的隨機過程如下:

圖1 多邊形顆粒形核點位置示意圖

a. 在某二維區(qū)域([xmin,xmax], [ymin,ymax])內(nèi),以顆粒最大粒徑dmax為邊長生成正方形網(wǎng)格,正方形網(wǎng)格形心位置為形核點的初始位置。

b. 對于每個形核點位置進行隨機偏移,即形核點坐標增加一個隨機偏移量:

xi=xini+N(-ks,ks)

yi=yini+N(-ks,ks)

(1)

式中:xini與yini為初始位置坐標;N為隨機數(shù);ks為N的取值范圍。本文采用[-ks,ks]的均勻分布,ks的取值不能大于dmax,即形核點位置偏移后仍在正方形網(wǎng)格內(nèi),一般取dmax的1/4～2/3較為適宜。同時ks的大小會對多邊形的形狀和大小產(chǎn)生影響,如圖1所示,隨著ks取值的增大時,Voronoi圖的各個多邊形之間的差異也隨之增大,即隨機性也增大。

2.2 孔隙率的控制

確定形核點位置后產(chǎn)生的Voronoi圖中多邊形顆粒是閉合的,相互之間不存在孔隙,此時為顆粒墊層最緊密的狀態(tài),可通過一定的手段對每個多邊形進行變換使得多邊形顆粒之間產(chǎn)生孔隙。

通過觀察可知,由于Voronoi圖中相鄰的多邊形之間共有一條邊,所以多邊形之間沒有孔隙,可以通過在共有邊上隨機生成新的點作為兩側(cè)的多邊形頂點,從而“制造”出孔隙。如圖2所示,在一條共有邊上隨機生成4個點,在隨機生成條件下,點與點之間的距離有可能非常小,有的可能接近重合,此時,可產(chǎn)生不同的顆粒間接觸的情況:①顆粒P與顆粒M之間的情況(邊M1M2全部包含在邊P1P2內(nèi));②顆粒L與顆粒M之間的情況(邊L11L12與邊M11M12部分包含);③顆粒P與顆粒L之間(顆粒P以頂點P11與顆粒L的邊L1L2相接觸);④顆粒L與顆粒N之間接觸邊頂點L8與頂點N2相互重合;⑤顆粒M與顆粒N之間頂點N5與M9重合,頂點N4與M10也重合;⑥兩個顆粒以點對點的方式接觸,由于隨機條件下出現(xiàn)這種情況的概率非常小,所以未在圖2中繪出。

圖2 多邊形顆粒頂點位置變換產(chǎn)生孔隙示意圖

由圖2可知,在共有邊上生成的新頂點位置決定了孔隙的大小,當新頂點均與原頂點重合時,模型的孔隙率最小為零,當顆粒與顆粒之間均是點對點接觸時,模型能達到最大的孔隙率,所以可通過控制新頂點產(chǎn)生的隨機過程,實現(xiàn)對顆粒墊層隨機重構(gòu)模型孔隙率的控制。本文采取以下方法進行孔隙率的控制:

a. 確定共有邊上產(chǎn)生新頂點中相距最遠的兩點間線段Lr，其長度為

lr=kr(μr,σr,0,1)lp

(2)

式中:lp為共有邊的長度;kr為取值范圍在0～1的隨機數(shù),本文采用正態(tài)分布產(chǎn)生kr;μr為kr的均值;σr為kr的標準差。

b. 通過確定線段Lr的起點來確定線段Lr的位置:

xp=x1+(lp-lr)kp(μp,σp,0,1)cosθ

yp=y1+(lp-lr)kp(μp,σp,0,1)sinθ

(3)

c. 在線段Lr內(nèi)再隨機產(chǎn)生2個點,將產(chǎn)生的新頂點(包含最大線段Lr的兩個端點)隨機劃分給共有邊兩側(cè)的多邊形,從而形成新的多邊形,若新多邊形之間相互不接觸,則重新進行隨機劃分,直至滿足相互接觸的條件。

根據(jù)上文可知,控制孔隙率的參數(shù)有2個,分別為kr和kp。隨機數(shù)的均值和標準差會對孔隙率產(chǎn)生影響,表1和表2分別給出了kr和kp在不同均值和標準差條件下的顆粒墊層的孔隙率。

表1 參數(shù)kr不同均值和標準差時的模型孔隙率

注:“/”表明在該組均值和標準差條件下,kr不在[0,1]內(nèi)的概率較高,故舍去。

表2 參數(shù)kp不同均值和標準差時的模型孔隙率

注:“/” 表明在該組均值和標準差條件下,kp不在[0,1]內(nèi)的概率較高,故舍去。

由于顆粒墊層的孔隙率主要由線段Lr的長度lr決定,所以從表1中可以看出,隨著kr的均值μr的增大,孔隙率不斷減小,而由于標準差σr只對kr的離散程度有影響,所以其對孔隙率的影響很小。從表2可以看出μp和σp取值大小與孔隙率的關(guān)系不大,說明Lr的起點位置與孔隙率的關(guān)系也不大,且從表2可以看出,其孔隙率基本與表1中μr=0.5一行的孔隙率相接近,說明只要當隨機過程的參數(shù)相同時,生成的隨機重構(gòu)模型孔隙率也相近,可見該模型的重復(fù)性較好。

2.3 粒徑級配的控制

對于多邊形顆粒來說,顆粒的幾何尺寸參數(shù)較多,不同文獻采用了不同的方法定義多邊形顆粒的粒徑[20-22]。本文采用多邊形顆粒邊界的最小Feret直徑作為顆粒粒徑,即顆粒邊界上不同方向上外切平行線間距的最小值,其值代表了顆粒所能通過的最小篩孔的孔徑。

從圖1可以看出,多邊形顆粒的形核點的位置不同,顆粒的粒徑就會變化。圖3為不同的ks時模型的粒徑級配曲線,隨著ks的增大,粒徑級配曲線略有變化,但變化的幅度很有限,粒徑分布范圍較窄?？梢娡ㄟ^改變ks,只能得到粒徑大小較為均勻的模型,并不能對模型的粒徑級配進行進一步的調(diào)整。本文通過設(shè)置顆粒的隨機分裂過程,產(chǎn)生粒徑較小的顆粒,以實現(xiàn)粒徑級配調(diào)整的目的。具體步驟如下:

圖3 不同ks值時粒徑級配曲線

步驟1通過顆粒粒徑來計算其分裂的概率Pf:

(4)

式中:dmax為顆粒的最大粒徑;dmin為顆粒的最小粒徑;di為該顆粒的粒徑。在0～1的范圍內(nèi)通過均勻分布取一個隨機數(shù)Pi,若Pi

步驟2在需分裂的多邊形顆粒內(nèi)以形心為原點建立一個坐標系,x軸與水平方向的夾角φ為0°～360°之間的一個隨機數(shù)(通過均勻分布產(chǎn)生)。

步驟3該坐標系的x軸與y軸會與多邊形顆粒的輪廓產(chǎn)生4個交點,通過連接多邊形顆粒形心與這4個交點形成4條線段,并在每條線段上通過均勻分布隨機取1個點,共產(chǎn)生4個隨機點,將每個需分裂的多邊形顆粒內(nèi)產(chǎn)生的隨機點和原有的形核點一起作為新的形核點繪制Voronoi圖,從而實現(xiàn)顆粒的一次隨機分裂。

由于顆粒分裂是隨機地在原有的顆粒內(nèi)部產(chǎn)生新的形核點,所以重新生成的模型中勢必產(chǎn)生了更小粒徑的顆粒,同時由于顆粒是按概率分裂,所以原有的大粒徑顆粒也能有部分得以保留,從而使得模型的粒徑級配曲線變得長而緩,不均勻系數(shù)明顯增大。若一次分裂后的粒徑級配曲線仍不滿足要求,則可重復(fù)顆粒分裂過程,直至滿足要求(見圖4,圖中Cu為不均勻系數(shù),Cc為曲率系數(shù))。當然也有可能存在所需的粒徑級配曲線恰好落在某一次分裂前后顆粒粒徑級配曲線之間的情況,此時,可按式(5)調(diào)整顆粒分裂的概率:

(5)

式中:kf為顆粒分裂控制系數(shù),取值范圍0～1,其值越小,代表大粒徑顆粒分裂的概率越小,分裂后級配曲線中大粒徑的顆粒占比越大,從而實現(xiàn)對模型級配進一步優(yōu)化調(diào)整的目的。圖5為不同kf時顆粒粒徑級配曲線。

圖4 顆粒隨機分裂后粒徑級配曲線

圖5 不同kf值時顆粒粒徑級配曲線

圖6為在長10 m、高2 m的二維范圍內(nèi)采用本文方法生成的顆粒墊層隨機重構(gòu)模型,由于控制孔隙率的操作只是在共有邊上產(chǎn)生新的頂點,對顆粒的最小Feret直徑幾乎沒有影響,所以顆粒墊層隨機重構(gòu)模型可先進行粒徑級配控制,得到滿足要求的級配曲線后,再按給定的孔隙率進行孔隙率控制,最終得到級配和孔隙率都滿足要求的顆粒墊層隨機重構(gòu)模型。

圖6 顆粒墊層隨機重構(gòu)模型的建立過程

3 模型顆粒與真實顆粒的形狀比較

3.1 顆粒形狀參數(shù)選擇

為了對比隨機重構(gòu)模型中顆粒與真實顆粒的形狀,首先要確定評價顆粒形狀的量化參數(shù)。關(guān)于顆粒形狀表達與評價方法的研究比較多,在尺度上可概括為3個層次:①整體輪廓(球狀、柱狀、片狀等);②棱角性(磨圓程度);③表面紋理[20-22]。對于顆粒墊層上土工膜的頂破脹破現(xiàn)象而言,其尺度一般只到第②層次,所以本文選擇了扁平度e(第①層次)和磨圓度q(第②層次)這2個形狀參數(shù)進行分析比較,計算公式如下:

(6)

(7)

式中:Fmax分別為最大、最小Feret直徑;PE為等效橢圓(扁平度和面積與顆粒都相同的橢圓)的周長;Pc為顆粒輪廓的實際周長。

3.2 真實顆粒形狀參數(shù)采集

目前數(shù)字圖像分析方法作為顆粒形狀測量的手段已經(jīng)廣泛應(yīng)用于巖土力學(xué)、材料力學(xué)等領(lǐng)域[20-22],本文采用高清數(shù)碼相機采集真實顆粒的數(shù)字圖像,并利用MATLAB軟件中的Image tool工具箱[23]對圖像文件進行分析從而獲得真實顆粒的形狀參數(shù)。圖像分析方法的基本原理是通過設(shè)置灰度的閾值,將顆粒與圖像的背景分離,形成黑白二值圖像,從而提取顆粒輪廓邊界,并對邊界進行分析,最終得到所需的形狀參數(shù),這就需要顆粒顏色與背景顏色產(chǎn)生較大的反差,才能保證顆粒輪廓邊界的準確采集和形狀參數(shù)的準確計算。圖7為部分真實顆粒的形狀參數(shù)采集過程,其中顆粒來源于浙江某抽水蓄能電站面板堆石壩料場,為灰?guī)r爆破料,顆粒數(shù)量為122個。

圖7 真實墊層顆粒圖像轉(zhuǎn)化為黑白二值圖像

3.3 形狀參數(shù)對比

圖8 顆粒扁平度與磨圓度散點圖

由于隨機重構(gòu)模型的顆粒較多,所以采用隨機抽樣的方法抽取122個模型顆粒與真實顆粒進行對比。圖8為模型顆粒與真實墊層顆粒的扁平度和磨圓度散點圖。由圖8可知，扁平度主要分布在1～2之間,其中落在1.2～1.6之間居多,模型顆粒樣本的扁平度均值為1.67,標準差為0.45,真實顆粒樣本的扁平度均值為1.44,標準差為0.28。磨圓度主要分布在0.85～0.95之間,且主要集中在0.90附近,磨圓度的分布范圍要比扁平度窄,模型顆粒的磨圓度均值為0.908 4,標準差為0.037,真實顆粒的磨圓度均值為0.885 1,標準差為0.049。

根據(jù)扁平度和磨圓度的定義可知,扁平度為最大Feret直徑與最小Feret直徑之比,所以其值必然大于1,磨圓度為等效橢圓周長與顆粒實際周長之比,其值必然小于1。由于顯著性檢驗中參數(shù)檢驗方法的前提是樣本來源于正態(tài)分布(分布范圍為[+∞,-∞])的總體[24],雖然目前尚未能完全確定扁平度與磨圓度符合何種分布,但從分布范圍來說顆粒扁平度和磨圓度均不是正態(tài)分布,所以對模型顆粒與真實顆粒的形狀參數(shù)統(tǒng)計分布的顯著性檢驗只能采用非參數(shù)檢驗方法。本文采用Ansari-Bradley方法[25]對模型顆粒與真實顆粒的扁平度和磨圓度進行顯著性檢驗,顯著性水平均為0.05。檢驗結(jié)果顯示:扁平度的顯著性檢驗指標p=0.081 6>0.05,接受假設(shè),即模型顆粒樣本與真實顆粒樣本的扁平度來源于同一總體分布;磨圓度的顯著性檢驗指標p=0.353 0>0.05,接受假設(shè),即模型顆粒樣本與真實顆粒樣本的磨圓度也來源于同一總體分布。

4 結(jié) 論

a. 參數(shù)ks用于控制Voronoi圖中形核點位置,可生成無孔隙的初始模型,且初始模型顆粒粒徑級配與參數(shù)ks之間關(guān)系不大,初始模型的粒徑較為均勻。

b. 通過在初始模型中多邊形顆粒的共有邊上設(shè)置新的頂點,可生成帶有孔隙的模型,模型的孔隙由參數(shù)kr和kp控制的隨機過程來產(chǎn)生,其中參數(shù)kr控制的隨機過程決定了模型孔隙率的大小。

c. 由于初始模型是根據(jù)顆粒的最大粒徑建立的,所以調(diào)整顆粒粒徑級配的過程是按顆粒粒徑大小來確定顆粒分裂的概率,顆粒分裂過程可不斷重復(fù),直至獲得滿足要求的粒徑級配,其中參數(shù)kf用于對顆粒分裂概率進行控制,以實現(xiàn)級配曲線在一次分裂中的微調(diào)。

d. 由于孔隙率的控制過程對顆粒的粒徑影響很小,所以一般先進行粒徑級配控制,再進行孔隙率控制,最終實現(xiàn)模型顆粒的粒徑級配以及孔隙率與真實顆粒的相近或相同。

e. 扁平度和磨圓度兩個形狀參數(shù)的比較表明,模型顆粒與真實顆粒的形狀較為相似,同時由于孔隙率和粒徑級配的可控制性,所以本文建立的顆粒墊層隨機重構(gòu)模型基本能夠還原真實的顆粒墊層,為進一步研究顆粒墊層上土工膜的頂破脹破奠定基礎(chǔ),也可為不規(guī)則顆粒離散元方法的建模提供參考,但本文模型只是二維隨機重構(gòu)模型,三維隨機重構(gòu)模型的構(gòu)建還有待進一步的研究。