王秀菊，石 崇，李德杰，梁邦炎，嚴晨宇

（1.南京交通職業(yè)技術學院 建筑工程系，江蘇 南京 211188；2.河海大學 巖土工程科學研究所，江蘇 南京 210098；3.中交第四航務工程局有限公司，廣東 廣州 510000）

1 引 言

顆粒離散元方法是近年來廣泛應用于巖土試驗模擬、變形破壞機制分析的有力工具，但該方法的準確性與合理性卻受制于顆粒模型的構(gòu)建。因此，在采用顆粒離散元方法進行數(shù)值建模工程分析時，如何得到合理、可靠的計算模型是獲得良好計算結(jié)果的前提條件。

目前顆粒離散元最常用的建模方法是膨脹法[1-2]，該方法先在指定區(qū)域內(nèi)生成小顆粒，然后增大顆粒半徑產(chǎn)生接觸力，迫使顆粒運動，直到充滿模型區(qū)域。但采用該方法構(gòu)建模型時膨脹系數(shù)難以控制，常常造成顆粒間有較大的重疊量，顆粒間、顆粒與邊界約束間的內(nèi)力較高。雖然當顆粒間的內(nèi)力小于指定的黏結(jié)力或者有邊界約束時，仍可用以模擬連續(xù)介質(zhì)的性質(zhì)。但分析邊坡滑坡、介質(zhì)破壞等動力學行為時，若顆粒重疊量較大，顆粒間黏結(jié)力不足以約束顆粒，顆粒間應變能的瞬間釋放就會造成大量顆粒飛溢出邊界，產(chǎn)生不準確的結(jié)果[3-4]。

針對這一問題，國內(nèi)外諸多學者嘗試采用顆粒相切條件進行數(shù)值模型的生成，并探討了各類算法[5-10]，但對于按指定密度、尺寸分布相符合、適合于任意形狀區(qū)域，并且體系中的顆粒都能處于平衡狀態(tài)，同時所生成的顆粒體系具有較高的接觸精度并與邊界耦合完好的充填算法還沒有出現(xiàn)。本文基于顆粒離散元原理，利用圓形顆粒相切生成與優(yōu)化算法，提出了一種基于Delaunay 算法的任意形狀范圍顆粒填充方法，實現(xiàn)了更完美模型的生成，可作為顆粒離散元計算的前處理模塊使用。

2 圓形顆粒模型構(gòu)造要求

離散元法是分析不連續(xù)介質(zhì)力學行為的數(shù)值方法，由Cundall[1]首先提出。Cundall和Strack[2]研制的圓盤、圓盤生成程序則奠定了顆粒離散元方法的基礎。

顆粒離散元屬于顯示求解的數(shù)值方法，顆粒間接觸力的計算廣泛采用Kelvin 模型（見圖1），顆粒間接觸力與顆粒間的重疊量及顆粒剛度成比例關系，其中，法向剛度描述的是總法向力和總法向位移之間的關系，切向剛度描述的是剪切應力增量和剪切位移增量之間的關系。

式中：Kn為法向剛度；ks為切向剛度；為法向接觸力；Un為法向變形重疊量；為切向相對位移；為相對切向力；ni為接觸方向矢量在法向的投影；i為接觸編號。

由式（1）可知，計算時需要輸入顆粒的剛度參數(shù)，如果所生成模型的初始重疊量較大，那么在離散元計算中，就會產(chǎn)生初始擾動力（應變能），與真實的物理模型產(chǎn)生較大偏差。在顆粒剛度確定的情況下，顆粒間的接觸力與重疊量成正比關系。當顆粒剛度較大且若邊界約束或者顆粒間黏結(jié)力消失時，顆粒間應變能的釋放就會造成顆粒飛溢出邊界的現(xiàn)象，掩蓋了材料變形破壞機制的實際情況，因此是數(shù)值模擬必需避免的。

圖1 Kelvin 模型顆粒間接觸力示意圖Fig.1 Contact forces between particles in Kelvin model

由于顆粒元離散元法模擬材料力學行為的第1步就是初始模型生成，即用顆粒元離散研究對象，在感興趣的區(qū)域內(nèi)生成顆粒體系。為了使模擬結(jié)果接近真實的物理過程，合理的模型需滿足如下要求：

（1）所生成的顆粒體系堆積密度能夠反映模擬對象的真實情況。如在模擬土體的力學行為時，需要生成密度較小的顆粒體系，而用離散元法模擬巖石材料的破壞與損傷時，則要求生成密度相對較大的顆粒體系。

（2）所生成模型的顆粒尺寸分布滿足指定要求。離散元模擬結(jié)果因尺寸分布的不同會有較大差異，因而顆粒的尺寸分布應與模擬對象的尺寸分布對應。

（3）顆粒間的接觸精度足夠高。要求模型生成算法所生成的顆粒體系中相鄰顆粒間的重疊量足夠小，需要盡可能提高算法精度從而減小顆粒間的重疊量。

（4）體系中顆粒與邊界緊密接觸，完整耦合。要使位于邊界處的顆粒與邊界相切，并且邊界與顆粒體系間的孔隙盡可能小，由此便于外載荷的施加，否則模擬結(jié)果與實際物理模型會產(chǎn)生較大差異。

（5）體系中顆粒處于受力平衡狀態(tài)。在生成的初始構(gòu)形中，每個顆粒所受合力為0，否則，在離散元計算中，即使在不施加外荷載的情況下，個別顆粒也會在重力作用下發(fā)生運動，與實際情況不符。

（6）適用于任意形狀的邊界。在實際工程的模擬中，所研究對象的形狀是千差萬別的，因而，模型的生成算法必須能夠適應各種邊界條件。

3 復雜顆粒離散元模型構(gòu)建

在某二維區(qū)域內(nèi)，為了在復雜模型域內(nèi)生成半徑范圍為[Rmin,Rmax]的顆粒體系，本文設想首先從一個種子出發(fā)逐步擴展填充到整個區(qū)域，并對模型邊界及孔隙內(nèi)進行再填充，從而提高模型密度。

3.1 填充種子位置

種子，在此定義為區(qū)域內(nèi)填充顆粒的初始位置。對二維問題，種子是由3 個圓盤(2D)組成，圓盤兩兩相切，在局部上達到最大密度，種子中每個圓盤的尺寸都是按照用戶給定的尺寸分布規(guī)律給出的，并且完全位于給定的區(qū)域內(nèi)。種子是生成復雜顆粒離散元模型的基礎，任意復雜模型都是通過在種子周圍不斷生成新顆粒，即膨脹填充來實現(xiàn)的。

在2D 情況下，如果已知兩個圓形顆粒的中心坐標與半徑，兩個圓的中心坐標為P2(x2,y2)，對應半徑分別為 r1、r2，則給定半徑為r 的新顆粒的位置 P(x,y) 可通過求解下式得到：

在求解式（3）時，同時對兩個根進行檢測，需要同時滿足在邊界內(nèi)且與已生成的其他顆粒不相交，如果都滿足要求，則同時被接受。

當種子中心位置與構(gòu)成顆粒的尺寸確定后，則以初始位置為中心，涵蓋3 個顆粒圓的最小圓稱為種子域（見圖2），其半徑R 定義為種子半徑，可由式（4）確定。

圖2 種子位置構(gòu)成Fig.2 Generation of a seed position

種子位置可以隨機，為了保證種子的所有顆粒均位于給定區(qū)域內(nèi)部，采用3.5 節(jié)所述距離控制時要求種子中心對復雜邊界的距離d＞R 。

3.2 局部Delaunay 化

種子生成后，把組成種子的顆粒與種子位置作為參考點，生成Delaunay 三角網(wǎng)格[11-12]，在2D 情況下，Delaunay 網(wǎng)格為一組三角形，去除網(wǎng)格中包含種子位置的三角形，以及有一條邊滿足Lij-ri-rj≥2Rmax的三角形，保留下來的三角形3個頂點處的圓盤即是用來生成新顆粒的單元，新顆粒生成后，判別是否與現(xiàn)有顆粒以及邊界重疊，如果與現(xiàn)有顆粒重疊則放棄，否則接受，如果與邊界重疊，則按所示算法直接計算顆粒的半徑與位置，如果半徑在[Rmin,Rmax]中，接受，否則放棄。重復上述過程，直至沒有新顆粒產(chǎn)生為止。把上述新顆粒生成的過程稱為膨脹，如圖3 所示。

圖3 邊界耦合顆粒生成Fig.3 Generation of boundary coupling particles

為了提高計算效率，不可能在每次膨脹中以所有顆粒位置為參考點進行Delaunay 剖分，為此，如圖4 所示，把在同一次膨脹過程中生成的顆粒標定為一個代，這個代號隨著膨脹的順序遞增，把組成種子顆粒的代號定義為0，把與當前代號相差為一定值內(nèi)的序列一起作為波前，把這個差值稱之為前緣厚度。如圖5 所示，選用具有一定厚度的波前，是因為相鄰的兩個波前間存在一個重疊區(qū)域，對這個區(qū)域可以進行多次充填，從而對復雜形狀區(qū)域進行充填時可以克服鎖死現(xiàn)象，并且實現(xiàn)重復充填，前緣厚度越大，重復充填的概率越高，可以得到較高的充填密度，有助于滿足用戶指定的堆積密度與尺寸分布。當然這樣做將會犧牲計算效率。

圖4 生成代填充Fig.4 Filling process by generation control

圖5 前緣內(nèi)部顆粒的局部Delaunay 網(wǎng)格Fig.5 Local Delaunay meshes among the front inner particles

在波前推進的過程中，新顆粒的半徑是隨機生成（或按用戶指定尺寸分布）的，某些情況下，新生成的顆粒可能無法放置在這個波前單元上。當顆粒不能放置時，該位置可能就會產(chǎn)生較大孔隙，如果波前是一個有一定厚度的區(qū)域，就有可能在后繼的循環(huán)中充填這個孔隙，從而增加了充填密度。

3.3 邊界耦合算法

當顆粒體系隨著種子的膨脹，不斷產(chǎn)生的新顆粒將會逐漸與邊界重疊，此時的顆粒與邊界不能完全接觸，則會影響顆粒體系與邊界的耦合性。為了保證顆粒體系與邊界均勻、平穩(wěn)的接觸，此時可以采用容忍性填充方法：放棄當前的顆粒尺寸，由與兩個圓盤以及邊界相切條件計算出顆粒的半徑與位置，此時的顆粒半徑一般小于設定的尺寸范圍，如此處理可導致產(chǎn)生的顆粒與期望范圍不一致，但如此可以使得顆粒體系與邊界緊密耦合，克服顆粒與邊界較大的問題，如圖3 所示。

當充填區(qū)域為多邊形時，在式（1）基礎上，新生成的顆粒還需要滿足下式：

式中：nx、ny為邊界邊的切向量；xb、yb為邊界邊上的一點。

3.4 內(nèi)部重填算法

當顆粒體系填充完整個區(qū)域后，提高顆粒體系的密度與減少處于不平衡狀態(tài)顆粒的比例，對于用離散元法模擬密度較高材料時是非常必要的。

此時，可將所有顆粒劃分成m×n 個柵格區(qū)域，提取每個柵格區(qū)域內(nèi)的顆粒圓心，進行局部Delaunay 三角化，此時為了保證顆粒能夠填充到孔隙中，最小顆粒半徑可以取為遠小于設計最小半徑，三角化后，對任意滿足要求的邊進行新顆粒的生成，一旦判斷有解，則逐步增大顆粒半徑，使得新生成顆粒恰好與第3 個顆粒相切而不疊加。對當前體系中每個顆粒重復上述過程，直至沒有新顆粒生成為止，如圖6 所示。

圖6 內(nèi)部重填顆粒過程Fig.6 Refilling process for inner particles

3.5 距離空間控制流程

新生成顆粒是否位于填充區(qū)域內(nèi)，是否為邊界顆粒，均可借助于距離來控制。由于現(xiàn)實需要填充的區(qū)域可以是任何形式，可以是圓、橢圓、任意多邊形等，也可能是多種幾何體的合成。此時如果采用距離空間描述，則可以較好地控制新顆粒的生成與邊界顆粒的搜索。點與常見幾何域邊界的距離計算可參考相關文獻[13]。

在平面任意聯(lián)通域，一個給定點p 至邊界S 的距離可定義為點靠近邊界的遠近程度，如式（6）所示，函數(shù)值為正，表明點位于S 內(nèi)，否則在界外。

如果是簡單邊界的組合，則可以利用簡單邊界的并集、差集、交集來描述。即

當判斷對象是半徑為r 的圓形顆粒時，如果d＞r，則顆粒位于模型區(qū)域內(nèi)；如果d=r，則顆粒與模型邊界相切；d ＜0，則顆粒位于模型區(qū)域外，舍棄；如果 r＞d＞0，則表明該顆粒為邊界顆粒，需要降低半徑以與邊界耦合。

3.6 完整的算法流程

利用種子擴展為復雜模型的完整算法如圖7 所示。

在距離控制下，一旦判斷某個顆粒為邊界顆粒，此時隨機產(chǎn)生的顆粒并不能保證與邊界相切，通常需要降低顆粒半徑限制，以保證產(chǎn)生的邊界顆粒恰好與邊界相切，以保證顆粒體系與邊界緊密耦合。這一相切半徑搜索過程可以采用優(yōu)化算法進行，設邊界顆粒的半徑為 r0，則待尋半徑處于［0,r0］之間，由于半徑無窮小的顆粒往往導致顆粒數(shù)目巨大，較為實際的辦法是在顆粒體系最小半徑 rmin基礎上，容忍邊界顆粒的半徑適當降低 k0倍，即使得邊界顆粒生成的范圍為［k0rmin,r0］之間(k0∈(0,1))。通過兩個端點半徑計算圓心坐標，求出圓心后計算其與復雜邊界的距離。然后構(gòu)造如式（6）所示的最優(yōu)化函數(shù)：

式中：d′為點距離復雜邊界的距離；r′為嘗試半徑。

圖7 完整算法流程Fig.7 A complete algorithm process

采用二分法等優(yōu)化算法可求出式（6）最小值對應的半徑，利用該半徑得到的圓心及半徑即為最優(yōu)的邊界顆粒。設置一較小值作為閾值，如果f 最小值大于該值，則表明該邊界顆粒的半徑過小，舍棄。

4 模型質(zhì)量與應用探討

4.1 模型質(zhì)量

為了驗證該方法生成顆粒離散元模型的質(zhì)量，采用某邊坡地質(zhì)分層如圖8 所示，定義最大/最小顆粒半徑為尺寸比率，顆粒面積與模型區(qū)域面積的比值為堆積密度，根據(jù)圖7 所示建模流程，采用尺寸比率為2，尺寸范圍為1.0～2.0 m，前緣厚度為5，對模型區(qū)域進行填充，共生成顆粒4 037，堆積密度為0.807，得到如圖9 所示顆粒離散元模型。

圖8 邊坡地質(zhì)分層Fig.8 Geological stratification of a slope

圖9 尺寸比率為2 時的充填初始模型（堆積密度0.807）Fig.9 Initial model with a geometric ratio of 2(packing density 0.807)

采用重填算法后，在圖9 模型基礎上對重填顆粒的最小半徑取最小設計尺寸1.0 m 的0.5 倍進行重填，得如圖10 所示的重填模型，模型顆粒數(shù)目5 989 個，堆積密度上升為0.854，如在此基礎上再采用0.1 倍半徑顆粒重填，則顆粒數(shù)目增加為15 987個，堆積密度增加為0.934，此時所有顆粒的平均疊加量為1.3×10-6m（見圖11），完全可以滿足計算的要求。

為了能得到更加密實的堆積密度，一種辦法是先采用大顆粒填充作為骨架，然后再采用細顆粒重填；也可以采用提高尺寸比例的辦法來實現(xiàn)，圖11為采用圖9 所示同樣的平均顆粒半徑0.75 m，但尺寸比率為5，尺寸范圍為0.25～1.25 m，前緣厚度為5 進行的模型生成，共采用顆粒49 675 個。得到模型的堆積密度0.836，比尺寸比率2 時提高0.029，但所需的顆粒數(shù)目大幅上升了10 倍。

圖10 圖9 模型采用0.5 倍最小半徑重填模型（堆積密度0.854）Fig.10 Refilling model with a tolerance ratio of 0.5 on the basis of model shown in Fig.9(packing density 0.854)

由于圖9～11 所示模型僅在整體區(qū)域內(nèi)填充生成，如果模型區(qū)域內(nèi)需要劃分為多種材料分層，此時可以將每種材料分層的邊界設置為一個任意多邊形，通過判斷顆粒的圓心是否位于多邊形內(nèi)，按照射線法判斷規(guī)則，如從顆粒圓心向右側(cè)做一條水平射線[13]，如果與某個多邊形的交點數(shù)目為偶數(shù)則該顆粒圓心位于該多邊形外側(cè)，如果為奇數(shù)，則位于多邊形內(nèi)部，利用多邊形的屬性進行材料屬性的設置。對圖11 所示模型，采用圖8 所示的材料分層后，得到地質(zhì)模型如圖13 所示。將該模型按照計算格式輸出，即可得到可用的數(shù)值模型。

圖11 圖10 模型采用0.1 倍最小半徑重填模型（堆積密度0.934）Fig.11 Refilling model with a tolerance ratio of 0.1 on the basis of model shown in Fig.10(packing density 0.934)

圖12 尺寸比率為5 時的充填初始模型（堆積密度0.836）Fig.12 Initial model with a geometric ratio of 5(packing density 0.836)

圖13 分層后的顆粒地質(zhì)模型Fig.13 Particle geological model after being stratified

4.2 模型對比

采用常用的膨脹法[1-2]，設計孔隙率為0.15、尺寸比率為2、半徑為1～2 m 生成如圖14 所示地質(zhì)模型，共用顆粒4 006，結(jié)果表明，在模型不同位置會蓄積不等的應變能，表現(xiàn)在模型表面各處應力變化。如采用文獻[3]伺服釋放法，則面積需要變化0.24%方能將如圖14 所示邊界應力降低至0.1 MPa。而本文所述方法，顆粒相切，應變能幾乎忽略不計，數(shù)值檢測如圖14 所示3 個位置法向應力均小于1 kPa，同時不需改變模型的外幾何形態(tài)，適應性更強。

圖14 膨脹法建立顆粒模型Fig.14 Particle geological model with expanding method

4.3 模型應用探討

本文所述方法的實現(xiàn)相對簡單，其填充過程主要是利用圓形相切規(guī)則與距離控制實現(xiàn)。只要生成的顆粒尺寸合適，可適用于任意二維聯(lián)通域內(nèi)顆粒模型的構(gòu)建。

通常顆粒體系中，如果尺寸比率過大，容易導致顆粒數(shù)目急劇增加，計算效率下降。但尺寸比率過小，雖然計算速度快，但容易產(chǎn)生較大的孔隙率。為了解決二者間的矛盾，一般尺寸比率處于[1，4]間較為合理。

顆粒生成后，需要根據(jù)本構(gòu)模型對顆粒進行屬性的設置，因此，需要對顆粒進行材料區(qū)分。材料邊界的屬性可根據(jù)實際條件進行接觸設置。

基于本文所述方法建立的模型接觸精度高，可以滿足巖土工程復雜介質(zhì)變形破壞的模擬需要。同時相應原理也可方便地推廣至三維領域，在生成三維地質(zhì)模型與復雜顆粒構(gòu)成方面發(fā)揮作用。

5 結(jié) 論

（1）根據(jù)顆粒離散元Kelvin 接觸力計算模型，研究了合理數(shù)值模型應具備的條件。提出在復雜模型域內(nèi)，從一個種子逐步擴展填充到整個區(qū)域，借助局部Delaunay 三角化網(wǎng)格控制新顆粒的生成，采用復雜幾何體距離控制顆粒與模型邊界的相對位置，對靠近模型邊界的顆粒進行容忍性優(yōu)化填充，從而增加了模型顆粒與邊界的耦合性。

（2）對初始生成的顆粒離散元模型，在其孔隙內(nèi)進行再填充，保證了填充顆粒至少與3 個顆粒相切，提高了顆粒體系的耦合性，同時提高了模型的密度。

（3）采用任意多邊形控制材料邊界，將模型材料的設置簡化為判斷點是否在多邊形內(nèi)，從而對任意的顆粒進行材料屬性設置，簡化了復雜模型材料的設置過程。

（4）與膨脹顆粒生成法相比，該方法生成模型重疊量小、顆粒間及顆粒-邊界相互耦合、填充率高，因此，顆粒黏結(jié)力破壞后不會造成飛溢現(xiàn)象，可適用于任意聯(lián)通域模型的生成，能更好地實現(xiàn)復雜巖土細觀介質(zhì)變形破壞機制的模擬與研究，并可推廣至三維模型的構(gòu)建。

[1]CUNDALL P A.PFC2DUsers’ Manual(version3.1)[M].Minnesota:Itasca Consulting Group Inc.,2004.

[2]CUNDALL P A,STRACK O D L.A discrete numerical model for granular assemble[J].Geotechnique,1979,29(1):47－65.

[3]SHI Chong,WANG Hai-li,WANG Sheng-nian,et al.A construction method of complex discrete granular model[J].Journal of Theoretical and Applied Information Technology,2012,43(2):203－207.

[4]石崇,王盛年,劉琳,等.基于數(shù)字圖像分析的冰水堆積體結(jié)構(gòu)建模與力學參數(shù)研究[J].巖土力學,2012,33(11):3393－3398.SHI Chong,WANG Sheng-nian,LIU Lin,et al.Structure modeling and mechanical parameters research of outwash deposits based on digital image analysis[J].Rock and Soil Mechanics,2012,33(11):3393－3398.

[5]BENABBOU A,BOROUCHAKI H,LAUG P,et al.Geometrical modeling of granular structures in two and three dimensions—Application to nanostructures[J].International Journal for Numerical Methods in Engineering,2009,80(4):425－454.

[6]FENG Y T,HAN K,OWEN D R J.Filling domains with disks:An advancing front approach[J].International Journal for Numerical Methods in Engineering,2003,56(5):699－713.

[7]LOHNER R,ONATE E.A general advancing-front technique for filling space with arbitrary objects[J].International Journal for Numerical Methods in Engineering,2004,61(12):1977－1991.

[8]KANSAL A R,TORQUATO S,STILLINGER F H.Computer generation of dense polydisperse sphere packings[J].Journal of Chemical Physics,2002,117(18):8212－8218.

[9]HAN K,FENG Y T,OWEN D R J.Sphere packing with a geometric based compression algorithm[J].Powder Technology,2005,155:33－41.

[10]JERIER J F,IMBAULT D,RICHEFEU V,et al.Packing spherical discrete elements for large scale simulations[J].Computer Methods in Applied Mechanics and Engineering,2010,199:1668－1676.

[11]胡金星,馬照亭,吳煥萍,等.基于格網(wǎng)劃分的海量數(shù)據(jù)Delaunay 三角剖分[J].測繪學報,2004,33(2):163－167.HU Jin-xing,MA Zhao-ting,WU Huan-ping,et al.Massive data Delaunay triangulation based on grid partition method[J].Acta Geodaetica et Cartographica Sinica,2004,33(2):163－167.

[12]朱合華,吳江斌.基于Delaunay 構(gòu)網(wǎng)的地層2D、2.5D 建模[J].巖石力學與工程學報,2005,24(22):4073－4079.ZHU He-hua,WU Jiang-bin.2D and 2.5D modeling of strata based on Delaunay triangulation[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(22):4073－4079.

[13]SCHNEIDER P J,EBERLY D H.計算機圖形學幾何工具算法詳解[M].周長發(fā)譯.北京:電子工業(yè)出版社,2003.