劉廣煜， 徐文杰， 周 乾， 馮澤康

(清華大學(xué) 水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

1 引 言

自然界中的顆粒材料，如沙子和礫石，是由幾何上帶棱角的非球形顆粒組成的，而顆粒的幾何形狀對于顆粒系統(tǒng)的物理力學(xué)特性有非常重要的影響。近年來，非連續(xù)介質(zhì)力學(xué)方法的發(fā)展為顆粒材料的研究提供了更多手段。在這些方法中，最常用的是離散元法，由Cundall等[1]首次提出，在海洋工程、巖土工程及化學(xué)工程等領(lǐng)域有廣泛的應(yīng)用[2-5]。

根據(jù)顆粒的幾何形狀，離散元法可分為球體離散元法和塊體離散元法。球體離散元法由于接觸檢測和接觸力計(jì)算簡單，得到廣泛的發(fā)展和應(yīng)用。為更加真實(shí)地體現(xiàn)顆粒幾何形狀的影響，球體離散元法在顆粒碰撞時引入了彎矩和扭矩[6]。但是彎矩和扭矩的計(jì)算需要引入額外的細(xì)觀參數(shù)，這使參數(shù)反演過程更加復(fù)雜[7]。球體顆粒簇模型是離散元法中另一種研究顆粒幾何形狀的有效方法[8]，但是該方法的計(jì)算成本較高。塊體離散元法解決了球體顆粒難以表征顆粒形狀的問題，能夠模擬顆粒的復(fù)雜形狀，且無需引入復(fù)雜的細(xì)觀參數(shù)[9,10]。

在離散元法中，顆粒假定為剛體，并允許顆粒之間有微小的重疊?；谥丿B部分計(jì)算嵌入深度、法線方向、接觸面積和接觸點(diǎn)等接觸特性，可求得接觸力。球體離散元法計(jì)算碰撞球體間接觸特性相對簡單，而塊體離散元法計(jì)算碰撞塊體間接觸特性難度大且耗時。公共平面法是目前應(yīng)用較為廣泛的一種塊體接觸算法[10,11]，該方法對接觸關(guān)系進(jìn)行分類，如點(diǎn)-點(diǎn)、點(diǎn)-棱和點(diǎn)-面等，分別求解接觸特性及接觸力。公共平面法為了確保數(shù)值模擬的穩(wěn)定性，會針對某些特殊接觸類型(如棱-棱接觸)進(jìn)行特殊處理，導(dǎo)致接觸計(jì)算耗時較長。眾所周知，多面體在頂點(diǎn)處存在法線方向不連續(xù)的問題[9]，因而兩碰撞多面體在頂點(diǎn)處發(fā)生小的相對運(yùn)動時，公共平面法不能保證接觸法線方向的平滑演化，可能導(dǎo)致在數(shù)值模擬中引入額外的能量。為解決這類問題，馮云田等[9,12]提出了能量守恒接觸模型，證明了碰撞的多面體顆粒利用重疊多面體計(jì)算接觸力是符合能量守恒定律的。重疊多面體的計(jì)算目前是一個難點(diǎn)。Yade作為一種開源和擴(kuò)展性強(qiáng)的離散元法框架，基于CGAL算法庫計(jì)算重疊多面體[13]，但法線方向是通過最小二乘擬合得到，因而該方法耗時且模擬的塊體顆粒數(shù)有限。為提高計(jì)算效率，基于圖像處理器GPU(Graphic processing unit)，Govender等[14]發(fā)展了并行離散元計(jì)算軟件 Blaze -DEM，該軟件顯式計(jì)算得到重疊多面體，能夠?qū)崿F(xiàn)百萬級多面體顆粒模擬計(jì)算。但由于接觸計(jì)算精度較低且接觸模型少，難以在巖土工程得到應(yīng)用。

本文在耦合模擬器的塊體離散元框架(CoSim-DEM)下開發(fā)了多面體顆粒的接觸重疊算法，CoSim-DEM是基于GPU并行計(jì)算的面向?qū)ο蟮目蓴U(kuò)展性強(qiáng)的數(shù)值模擬仿真平臺。接觸重疊算法是基于幾何對偶理論[15]衍生的多面體接觸理論，是一種更精確的接觸算法?；谥丿B多面體，可以計(jì)算接觸點(diǎn)、接觸面積、法線方向和嵌入深度等接觸特性，進(jìn)而用于接觸力和力矩求解計(jì)算。

2 多面體接觸重疊算法

在塊體離散元中，顆粒的迭代計(jì)算主要包括鄰居搜索、接觸檢測、接觸力計(jì)算和插值計(jì)算四個步驟。其中，鄰居搜索就是通過均勻網(wǎng)格法[16]或邊界包圍盒層次樹法[17]等找到每個顆粒的鄰近顆粒列表。接觸檢測基于鄰居搜索得到的列表進(jìn)行重疊多面體計(jì)算，進(jìn)而求解接觸力。多面體接觸重疊計(jì)算算法流程如圖1所示。

圖1 多面體重疊接觸算法流程

本文僅考慮凸多面體的接觸算法，因?yàn)閷τ谌我庑螤畹膹?fù)雜塊體，可以由若干個凸多面體拼接而成，因此凸多面體的算法同樣適用于復(fù)雜塊體的碰撞?；趲缀螌ε祭碚撗苌鰜淼亩嗝骟w接觸理論，將位于初始坐標(biāo)系的兩多面體轉(zhuǎn)換到對偶空間中，此時兩多面體的所有面轉(zhuǎn)換成了對偶空間中的點(diǎn)，在對偶空間中計(jì)算所有點(diǎn)的凸包，再將該凸包轉(zhuǎn)換到初始坐標(biāo)系下生成一個新的多面體，此時初始坐標(biāo)系下的這個多面體即為兩碰撞多面體的重疊體。多面體接觸重疊算法基于該理論，融合了GJK算法和快速凸包算法，基于重疊體計(jì)算接觸面積、接觸法線方向、嵌入深度和接觸點(diǎn)等接觸特性，可以提高計(jì)算效率，保證計(jì)算結(jié)果精確且高效。該算法主要分為GJK算法、幾何對偶算法、快速凸包算法和接觸特性計(jì)算四步。

2.1 GJK算法

對于兩個塊體A和B，其接觸關(guān)系在歐氏空間中可以分為碰撞、觸碰和分離三種類型(圖2)。

圖2 兩多面體間接觸狀態(tài)

GJK算法計(jì)算成本低且通用性高，通過迭代，能夠檢查兩個塊體是否發(fā)生碰撞。該算法利用多面體A和B各自的頂點(diǎn)點(diǎn)集a和b，通過尋找Minkowski差分A-B(A-B={a-b,a∈A,b∈B})內(nèi)部或外部邊界上的點(diǎn)V，找到距離原點(diǎn)O最小的點(diǎn)，如圖2所示。在確定兩塊體的接觸狀態(tài)為碰撞或觸碰時，執(zhí)行仿射變換將點(diǎn)V轉(zhuǎn)換為v′(v′∈A，v′∈B)，該點(diǎn)將用于后續(xù)對偶空間變換算法。

2.2 對偶空間變換算法

GJK算法篩除了分離狀態(tài)的多面體對，而對于存在接觸的多面體對，需要對其進(jìn)行幾何對偶化，將其轉(zhuǎn)化到對偶空間中，用于后續(xù)凸包的計(jì)算。幾何對偶化具體操作是將初始坐標(biāo)系的多面體的面變換到對偶空間的點(diǎn)。對偶空間的點(diǎn)vdual,j可表示為

d=(ci-v)·ni，vdual,j=ni/d

(1，2)

式中ci和ni分別為多面體每個面的面心和法線方向，v為重疊體內(nèi)部的任意一點(diǎn)。GJK算法計(jì)算出的點(diǎn)v′不能直接作為重疊多面體的內(nèi)部點(diǎn)，因?yàn)関′有可能在重疊體的面上(圖3)。因此，Koziara[18]提出了向重疊多面體內(nèi)部延伸點(diǎn)v′的迭代算法，如圖3所示，最終延伸方向nt可以是若干個多面體的面法線反方向(圖3的n1，n2和n3)的組合，進(jìn)而得到重疊體內(nèi)部點(diǎn)v。延伸的距離由自定的閾值決定，可根據(jù)精度要求進(jìn)行設(shè)置。根據(jù)式(1，2)，基于重疊體內(nèi)部的該點(diǎn)v，多面體A和B的面將轉(zhuǎn)換為對偶空間中的散點(diǎn)vdual,j。

圖3 重疊體內(nèi)部點(diǎn)v計(jì)算算法

2.3 快速凸包算法

在對偶空間中，通過快速凸包算法，搭建包含所有散點(diǎn)vdual,j的最小凸多面體。圖4顯示了基于給定點(diǎn)集構(gòu)建最終凸包的整個迭代過程。如圖4(a)所示，創(chuàng)建一個初始四面體作為初始凸包，給定散點(diǎn)集由此分成兩組，第一組是初始四面體四個不共面的頂點(diǎn)(圖4(a)的VT k)；另一組是位于四面體面外(圖4(a)的fk)的點(diǎn)(圖4(a)的VF i,j)。初始四面體搭建完成后，需要迭代計(jì)算建立整個凸包。

為說明迭代算法，迭代計(jì)算凸包的其中一個迭代步如圖4(b)所示。VF 1,1是f1外部且距離f1最遠(yuǎn)的點(diǎn)。算法中，f1將標(biāo)記且刪除，因?yàn)樵撁娌皇峭拱倪吔?。遍歷f1的所有相鄰面(圖4(b)的f2，f3和f4)，滿足式(3)的面也會標(biāo)記刪除。然后基于點(diǎn)VF 1,1和已刪除面(f1)的邊，創(chuàng)建新的面fN 1，fN 2和fN 3。此時更新后的凸包由原始的面(圖4(b)的f2，f3和f4)和新創(chuàng)建的面(圖4(b)的fN 1，fN 2和fN 3)構(gòu)成。最終，當(dāng)凸包外部沒有點(diǎn)時，迭代終止，此時整體凸包構(gòu)造完成(圖4(d))。

圖4 快速凸包算法

nf·(VF i,j-cf)>0

(3)

式中cf和nf分別為凸包表面的面心和法線方向。

2.4 接觸特性計(jì)算

圖5 兩碰撞多面體的重疊體

(4)

(5)

kn=(EA+EB)·a/(RA+RB)

(6)

(7)

ks=(EA·vA+EB·vB)·a/(RA+RB)

(8)

3 多面體接觸重疊算法驗(yàn)證

為驗(yàn)證接觸重疊算法的適用性和準(zhǔn)確性，模擬分析了兩個算例，分別是顆粒碰撞的基準(zhǔn)測試和砌體結(jié)構(gòu)破壞過程，對比了模擬結(jié)果與試驗(yàn)結(jié)果。

3.1 塊體碰撞測試

構(gòu)造了三種塊體碰撞類型的基準(zhǔn)測試，分別是點(diǎn)-點(diǎn)、點(diǎn)-面和棱-棱碰撞，如圖6所示。在這些基準(zhǔn)測試中，具有相同質(zhì)量的顆粒在重力作用下從相同高度(1 m)落下，并與下方固定的顆粒相撞，此時回彈系數(shù)設(shè)置為0.1。此外，顆粒的幾何形狀均采用軸對稱多面體，包括八面體、六面體和三棱柱。為研究回彈系數(shù)的影響，針對點(diǎn)-面碰撞類型采用了三個回彈系數(shù)進(jìn)行模擬，并對比了不同回彈系數(shù)下碰撞力的演化過程。表1是模擬中使用的計(jì)算參數(shù)。

圖6 顆粒碰撞基準(zhǔn)測試

表1 顆粒碰撞的數(shù)值計(jì)算參數(shù)

圖7展示了三種不同接觸類型在顆粒碰撞時接觸力的演化。如圖7(a,b)所示，在點(diǎn)-點(diǎn)接觸的情況下，盡管嵌入深度已經(jīng)大于八面體顆粒尺寸的15%，接觸力僅作用在垂直方向(Z方向)。在三種碰撞類型下(點(diǎn)-點(diǎn)、點(diǎn)-面和棱-棱接觸)，X和Y方向的接觸力都相同且都等于0，且由于黏滯力消耗能量，多次碰撞后，上方顆粒均可以靜止于下方顆粒上方。由于顆粒形狀均為軸對稱圖形，碰撞力在X和Y方向分量為0是合理的。進(jìn)而對比了三種碰撞類型在垂直方向的接觸力的演化，如圖7(c)所示，由于顆粒碰撞中黏壺的作用，部分能量得到釋放，最終接觸力可以收斂到上部顆粒的重力。不同碰撞類型接觸力收斂速度不同，由式(6)可發(fā)現(xiàn)，接觸面積越大會導(dǎo)致接觸剛度越大，點(diǎn)-面接觸由于接觸面積最大，因此收斂速度最快。回彈系數(shù)對接觸力的收斂速度有更為顯著的影響，較小的回彈系數(shù)可以明顯加快接觸力收斂速度，如圖7(d)所示。在回彈系數(shù)為1.0時，發(fā)生多次碰撞后，上部顆粒都能夠回到初始位置，每次碰撞力的演化均一致，由于顆粒系統(tǒng)中無摩擦和黏滯力導(dǎo)致的能量損失，因此碰撞過程符合能量守恒，驗(yàn)證了碰撞過程中沒有額外的能量增加。該算例充分驗(yàn)證了多面體接觸重疊算法的適用性，可廣泛用于多體接觸模擬計(jì)算。

圖7 碰撞測試下不同接觸類型的接觸力和顆粒位置的演化

3.2 砌體結(jié)構(gòu)破壞過程

砌體結(jié)構(gòu)在結(jié)構(gòu)工程中十分常見。為進(jìn)一步驗(yàn)證算法在計(jì)算顆粒系統(tǒng)的旋轉(zhuǎn)、位移和接觸力等方面的準(zhǔn)確性，采用了文獻(xiàn)[19]的砌體結(jié)構(gòu)破壞試驗(yàn)進(jìn)行模擬對比。

Portioli等[19]使用的模型是干接縫凝塊巖塊制成的砌體面板(圖8)。面板由兩種砌塊組成，尺寸分別為100 mm×200 mm×50 mm和100 mm×100 mm×50 mm。面板的尺寸是高400 mm、寬1100 mm和厚100 mm。底部左側(cè)設(shè)置長度為 600 mm 的固定支撐臺，支撐臺上方的磚塊僅可以沿水平方向移動。底部右側(cè)設(shè)置長度為500 mm的升降臺，可在水平和垂直方向自由移動。為監(jiān)測升降臺上的承載力和位移，設(shè)置機(jī)械千分表和臺秤分別測量升降臺的垂直位移和承載力。試驗(yàn)過程是升降臺每次緩慢降低1 mm，待面板運(yùn)動穩(wěn)定后，記錄升降臺的位移變化，同時記錄升降臺的承載力變化。根據(jù)模型試驗(yàn)，相應(yīng)的數(shù)值模型如圖8(c)所示，塊體之間沒有設(shè)置黏聚力，其中固定支撐結(jié)構(gòu)和升降臺分別由固定塊體和可移動塊體表示。在模擬過程中，可移動塊體以10-4m/s的速度向下移動?？梢苿訅K體與上方塊體之間的摩擦角設(shè)置為零，以確保在水平方向上有自由位移，塊體顆粒之間的摩擦力設(shè)置為20°，計(jì)算過程中監(jiān)測可移動塊體受到的外力變化，進(jìn)而與試驗(yàn)結(jié)果對比。其他參數(shù)列入表2。

圖8 砌體結(jié)構(gòu)模型試驗(yàn)及離散元模型[19]

表2 砌體結(jié)構(gòu)的數(shù)值計(jì)算參數(shù)

數(shù)值結(jié)果和模型試驗(yàn)結(jié)果在面板破壞時的比較如圖9所示。可以看出，數(shù)值結(jié)果與模型試驗(yàn)結(jié)果吻合較好。根據(jù)變形特征，砌體結(jié)構(gòu)中部發(fā)生滑動和轉(zhuǎn)動破壞(圖9(b)的第2部分)，表明該算法可以準(zhǔn)確模擬塊體的旋轉(zhuǎn)和平移。左側(cè)部分(圖9(b)的第1部分)表現(xiàn)得像不可變形體，位移非常??；右側(cè)部分(圖9(b)的第3部分)也類似不可變形體，但通過階梯狀裂縫與中間部分分開。圖10顯示了作用在升降臺的承載力隨垂直位移的演化。從數(shù)值結(jié)果與試驗(yàn)結(jié)果的對比來看，兩者吻合較好。而且，砌體結(jié)構(gòu)發(fā)生破壞時數(shù)值模擬計(jì)算得到的承載力約為104 N，這也與試驗(yàn)結(jié)果相一致。該算例充分驗(yàn)證了算法的精確性。

圖9 可移動塊體豎直位移20 mm的模型試驗(yàn)與模擬結(jié)果對比[19]

圖10 作用在升降臺的承載力隨位移的演化

4 結(jié) 論

離散元法已廣泛用于顆粒材料的物理和力學(xué)研究。在塊體離散元法的框架CoSim-DEM下，開發(fā)了多面體接觸重疊算法，該算法能夠顯式計(jì)算得到精確的重疊多面體。根據(jù)計(jì)算流程依次分為四個子步驟，分別是GJK算法、對偶空間變換算法、快速凸包算法和接觸特性計(jì)算。基于精確的重疊多面體，可以直接計(jì)算法線方向、嵌入深度、接觸點(diǎn)和接觸面積等接觸特征，更合理地計(jì)算接觸力。

采用了顆粒碰撞的基準(zhǔn)測試驗(yàn)證該算法的適用性和準(zhǔn)確性。通過計(jì)算結(jié)果可以發(fā)現(xiàn)，該算法能夠準(zhǔn)確計(jì)算各種接觸類型的接觸力，如點(diǎn)-點(diǎn)接觸、點(diǎn)-面接觸和棱-棱接觸，并發(fā)現(xiàn)回彈系數(shù)越小，接觸力收斂越快。進(jìn)而對砌體結(jié)構(gòu)的破壞過程進(jìn)行了模擬，并與模型試驗(yàn)進(jìn)行了對比，表明數(shù)值模擬對破壞機(jī)理的揭示和再現(xiàn)承載力的演化非常有效。本文開發(fā)的多面體接觸重疊算法避免了多面體之間接觸搜索的繁瑣幾何判斷，其計(jì)算結(jié)果具有很好的可靠性，可應(yīng)用于非規(guī)則顆粒材料的離散元模擬。