基于統(tǒng)一標(biāo)定的勢接觸力計算

嚴(yán)成增，鄭 宏，葛修潤

（中國科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點實驗室，湖北 武漢 430071）

1 引 言

巖體中存在節(jié)理、斷層等大量不連續(xù)面，對這類問題，基于連續(xù)介質(zhì)力學(xué)的數(shù)值類方法，如有限元通常難于處理。為了處理這類非連續(xù)介質(zhì)力學(xué)問題，誕生了離散元類數(shù)值方法。其中最有代表性的為Cundall[1]提出離散元方法和石根華[2]提出非連續(xù)變形分析方法（DDA）。這兩種方法在計算接觸力方面的共同點是均引入了剛度參數(shù)，通過嵌入量計算接觸力，不同的是Cundall提出的離散元允許塊體間相互嵌入，而石根華[2]提出DDA則通過開閉迭代，使得塊體間不嵌入。鄔愛清等[3]對DDA中接觸力的計算過程予以詳細(xì)的介紹，指出通過目前的接觸狀態(tài)，按總體方程求解，當(dāng)系統(tǒng)所有接觸部位滿足不嵌入無張拉準(zhǔn)則時，此時計算出的接觸力分量則為真實接觸力分量，否則通過增刪彈簧繼續(xù)迭代求解。鄭宏等[4]為避免引入剛度參數(shù)和開-閉迭代，將互補理論引入到DDA中，將DDA中各潛在接觸對以互補形式來表現(xiàn)，通過求解非線性互補方程組來獲得問題的解答，該方法與傳統(tǒng)的接觸力計算方法有很大的不同。徐棟棟等[5]將Munjiza基于勢的接觸力計算方法引入到流形元中，用于處理流形元的接觸問題。以上計算接觸力的方法均存在一些問題，例如Cundall的離散元方法，不能處理角-角接觸的問題，在其開發(fā)的二維版本的離散元程序UDEC中通過圓角化來規(guī)避這一問題[6]。在三維情形下，雖提出了公共面法來求解接觸力，但仍存在公共面不惟一[7]，以及無法處理某些極端接觸問題[8]。石根華的DDA也不能很好地處理角-角接觸的問題，同時當(dāng)問題擴展到三維狀態(tài)時，接觸判斷變得相當(dāng)復(fù)雜，這也是三維DDA發(fā)展較為緩慢的一個重要原因。

與以上各種計算接觸力的方法不同，Munjiza在其提出的有限元-離散元（FDEM）耦合分析方法中，采用了一種非常新穎的基于勢的接觸力計算方法[9－10]。有關(guān)該方法的相關(guān)發(fā)展可參見文獻[11－12]。該方法計算的接觸力為分布力，不會產(chǎn)生應(yīng)力集中，對點-點接觸無需作特殊處理，而常規(guī)的離散元(DEM)為了處理點-點接觸及避免角點處產(chǎn)生應(yīng)力集中需要對角點進行圓角化?；趧莸慕佑|力計算法雖然規(guī)避了傳統(tǒng)DEM在處理點-點接觸上存在的困難，但仍存在勢的物理意義不夠明確，相同的嵌入量計算的接觸力大小不一致等問題。在已有研究的基礎(chǔ)上，本文提出了一種基于統(tǒng)一標(biāo)定的勢接觸力計算方法，該方法重新定義了勢函數(shù)，使得勢的物理意義更為明確，該方法中勢就是對嵌入的一種表征；同時，在相同的嵌入量時，計算的接觸力大小是一致的。在解決了原有勢接觸力計算方法物理意義不夠明確，對于相同嵌入量，其接觸力的大小不相同等問題的同時，保留了原有Munjiza勢接觸力方法的所有優(yōu)點，即接觸力為分布力，對點-點接觸無需特殊處理，接觸力所做的總功為0，系統(tǒng)的能量不會因為接觸力的存在而增加或者減少。

2 FDEM勢接觸力

如圖1所示，將接觸塊體分別稱之為目標(biāo)塊體cβ 和靶塊體βt，目標(biāo)塊體嵌入到靶塊體的面積dA所引起的接觸力為[10]

式中：Pt為重疊區(qū)域內(nèi)位于靶塊體 βt內(nèi)的一點；φt(Pt)分別為 Pt點在靶塊體 βt內(nèi)的勢，grad為梯度；pn為法向罰參數(shù)。

圖1 接觸力計算示意圖Fig.1 Potential based contact force between two elements

同理，靶塊體βt嵌入到目標(biāo)塊體βc中的面積dA 引起的接觸力為

特別說明，式（2）計算的接觸力為目標(biāo)塊體βc提供勢對靶塊體βt的作用力，于是靶塊體βt對目標(biāo)塊體βc的反作用力大小與式（2）相同，方向相反，即為

于是，重疊區(qū)域dA 引起的總接觸力（靶塊體對目標(biāo)塊體）為[10]

于是，整個重疊區(qū)域所引起總的接觸力可通過對上式積分得到

對式（5）應(yīng)用格林公式，可將面積分轉(zhuǎn)化為對重疊區(qū)域邊界的線積分：

Munjiza等[9－10]指出，勢函數(shù)φ 在離散的單元邊界上應(yīng)為常量，只有這樣才能保證接觸力做的總功為0，這也是基于勢的接觸力對勢函數(shù)定義的惟一要求。據(jù)此，Munjiza[10]定義了如下勢函數(shù)：

式中：P為靶單元中長度為L 的邊上的任意一點。

如圖2所示，Ai(i=1，2，3)為P 點與三角形3條邊所構(gòu)成的子三角形面積，A為三角形的面積。這樣在三角形中心點的勢等于1，而在邊界上的勢等于0。為了便于理解，以三角形的中心點C 將三角形分為3個子三角形，要據(jù)式（7）求解勢，只需先判斷點P 在那個子三角形內(nèi)，例如點P 在C12子三角形內(nèi)，則該點的勢據(jù)式（7）可進一步改寫為

式中：hP-12表示點P 到12邊的距離；h0-12表示0點到12邊的距離，亦即三角形12邊的高。

圖2 三角形單元內(nèi)一點P 的勢Fig.2 Potential at point P in a triangle element

這種勢函數(shù)的定義，雖然滿足了在邊界上勢為常量的要求，但存在物理意義不夠明確，相同的嵌入量計算出的接觸力大小不等的問題。

對該勢函數(shù)稍作分析，便可發(fā)現(xiàn)其中的問題，如圖3所示，左右兩組接觸單元，除了靶單元的大小不同外，兩組接觸單元的嵌入量完全相同，其中目標(biāo)單元標(biāo)號為L 的邊與靶單元的12、1′ 2 ′邊平行，且該邊到12、1′ 2 ′邊的距離相等，亦即嵌入量相同，L 邊的外法線向量為 nL。以圖3靶單元中長度為L的邊為例，據(jù)式（6）、（8）對該邊積分得到的勢接觸力，其中圖3中右邊的長度為L 的邊積分得到的勢接觸力為

式中：pn為法向罰參數(shù)。

同理，圖3中左邊的長度為L 的邊積分得到的

勢接觸力為

圖3 嵌入量相同的兩組接觸Fig.3 Embedding the same amount of two group contacts

此處僅僅是對重疊區(qū)域邊界的一條邊作了分析，總的接觸力需要對整個重疊區(qū)域的邊界進行積分，與上述分析類似，對圖3所示的情況，原有程序計算出目標(biāo)單元受到的總接觸力如表1所示（均已分配到單元的節(jié)點上），采用的罰參數(shù)pn=67.3 GPa。表1表明，在嵌入量完全相同的情況下，目標(biāo)單元所受的總接觸力并不相同，而這是與事實不相符合的。也就是說，此處的勢函數(shù)的定義物理意義是不明確的，無法對嵌入量這一指標(biāo)基于統(tǒng)一的標(biāo)準(zhǔn)進行表征。為此，在接下來的部分，提出了基于統(tǒng)一標(biāo)定的勢接觸力計算方法，采用了新的勢函數(shù)定義，使得勢可以基于統(tǒng)一標(biāo)準(zhǔn)對嵌入量進行表征，在相同嵌入量的情形下，計算出的接觸力大小相等。

表1 基于原有勢函數(shù)計算的接觸力Table 1 Results of contact force with original potential function

3 統(tǒng)一標(biāo)定的勢接觸力

3.1 基本思路

為了解決原有勢接觸力計算方法存在的問題，對勢函數(shù)進行了重新定義，為了使得勢能夠?qū)η度肓炕诮y(tǒng)一標(biāo)準(zhǔn)進行表征，三角形單元內(nèi)一點的勢與該點到3條邊的最短距離成正比，于是可以定義如下勢函數(shù)：

式中：k為比例常數(shù)，hP-01、hP-12、hP-20分別表示點P 到三角形01邊、12邊、20邊的垂直距離，如圖4所示。

圖4 三角形單元內(nèi)一點P 的新勢函數(shù)Fig.4 A new potential function at point P in a triangle element

為了使得整個離散系統(tǒng)的單元的勢基于統(tǒng)一標(biāo)準(zhǔn)對嵌入量進行表征，同時保證勢φ為無量綱量，于是定義一個標(biāo)尺長度H，認(rèn)為當(dāng)一點的嵌入量（一點到三角形3條邊的最大距離即為該點的嵌入量）為H 時，該點的勢為1。于是據(jù)式（11）有k=1/H，這樣勢函數(shù)的表達式可進一步寫為

此處需要特別說明，H為統(tǒng)一標(biāo)定嵌入量，可以根據(jù)需要設(shè)定H 的值。為了使得法向罰參數(shù) pn具有與傳統(tǒng)DEM中法向剛度類似的意義，取H=1 m，這樣法向罰參數(shù)就表示一點嵌入為1 m時，提供與法向罰參數(shù)大小相同的壓力，H 也可設(shè)置為所有三角形單元中最短的高，來防止小尺寸的單元過大地嵌入。

為了便于理解，對式（12）作進一步的分析，依據(jù)三角形3條角平分線的交點（因為角平分線上的點到兩邊的距離相等），該交點即為三角形的內(nèi)心I，將三角形分為3個子三角形，如圖4所示。這樣要求一點的勢，只需先判斷該點位于那個子三角型內(nèi)，比如點P 位于I12子三角形內(nèi)，則P 點勢據(jù)式（12），可寫為

據(jù)以上分析可知，三角形單元位于內(nèi)心的點的勢最大，但該點處的勢不一定為1，同時不同三角形單元內(nèi)心點處的勢的大小不同，這與Munjiza定義的勢函數(shù)在三角形內(nèi)的中心處勢最大，且不同三角形在中心處的勢均為1是不同的。比較式（8）和式（13）可以發(fā)現(xiàn)，兩者的主要不同在于，一個分母是 h0-12，一個分母是H，由于 h0-12的大小會隨三角形的不同而不同，而采用H，則因為所有三角形單元都基于統(tǒng)一的標(biāo)準(zhǔn)，這樣的勢函數(shù)定義基于統(tǒng)一的標(biāo)準(zhǔn)表征了嵌入量，物理意義明確。同時，新定義的勢函數(shù)在三角形的邊界上勢為0，完全滿足勢函數(shù)在邊界上為常量這一要求。

采用新的勢函數(shù)的定義，可以使得在嵌入量相同的情況下，計算得到的接觸力大小相同。仍然以圖3為例，據(jù)式（6）和式（13）對目標(biāo)單元中長為L 的邊積分得到的接觸力為

由式（14）可知，對圖3所示的左右兩種接觸情況，因為嵌入量相等，即hP-12相同，這樣兩種情況計算得到的接觸力相同，此處僅對重疊區(qū)域的一條邊界作了分析，其余邊界也可類似分析，對圖3左右兩種接觸情況，基于新的勢函數(shù)定義計算得到的總的接觸力相同。這一點可在第4節(jié)中算例1中得到進一步證實。

3.2 實施方案

由3.1節(jié)論述可知，本文僅僅是定義了新的勢函數(shù)，再無其他更改，以圖5所示，要求靶單元對目標(biāo)單元AB 邊的作用力，由于勢函數(shù)的值在AB邊上為分段線性函數(shù)，基于原有勢函數(shù)定義，I 點為三角形的中心點，I2為三角形01邊的中線，即勢在I2邊與AB 邊的交點處發(fā)生轉(zhuǎn)折，采用新的勢函數(shù)的定義，只需將三角形的中心點I 的坐標(biāo)替換為三角形的內(nèi)心點的坐標(biāo)，來求I2與AB 邊的交點，并采用式（13）計算勢即可。整個過程非常容易實施，只需對原有程序作非常小的改動即可。

圖5 靶單元對目標(biāo)單元AB 邊的接觸力計算Fig.5 Contact force of target element to edge AB of contactor element

如圖5所示，在求得AB 邊與角平分線I2的交點P1以及AB與三角形01邊的交點P2后，依據(jù)式（5）計算靶單元對目標(biāo)單元AB 邊的接觸力為

圖5中，φt在AB 邊上是分段線性分布的，故為圖5中所示區(qū)域的面積，于是可進一步寫為

φt(A)、φt(P1)、φt(P2)由式（14）可得

4 算 例

4.1 算例1

對圖3所示的算例，采用基于統(tǒng)一標(biāo)定的勢接觸力計算方法進行計算，采用的罰參數(shù)仍為pn=67.3 GPa 。對圖3左右兩組接觸對計算得到對目標(biāo)單元的總接觸力（均已分配到節(jié)點上）如表2所示。由表可知，采用本文的方法計算，對左右兩種接觸情形，計算的接觸力完全一致。

表2 基于新勢函數(shù)計算的接觸力Table 2 Results of contact force with new potential function

4.2 算例2

圖3所示的算例給出的是一種較為特殊情況，該算例主要是便于進行理論分析，不能完全體現(xiàn)出基于統(tǒng)一標(biāo)定的勢接觸力計算方法的全部優(yōu)勢，同時也無法完全暴露原有勢函數(shù)定義所存在的不合理之處。這里給出一個更為通用的算例，如圖6所示，采用Munjiza原有的勢函數(shù)定義方法，在目標(biāo)單元和靶單元重疊區(qū)域保持不變的條件下，靶單元的AB 邊變?yōu)锳′ B 邊，或者CD 邊變?yōu)镃D′邊，計算出的總的接觸力都會發(fā)生變化，因為AB 邊變?yōu)锳′ B邊，或者CD 邊變?yōu)镃D′邊，使得三角形的高發(fā)生了變化，必然導(dǎo)致重疊區(qū)域的勢發(fā)生變化，從而使得在重疊區(qū)域完全一致的情形下，計算出的總接觸力并不相同，而這違背了基本的常識，是讓人無法理解和不可接受的。正確的情況是，只要接觸的局部區(qū)域相同，那么計算出的接觸力大小就應(yīng)該一致。如圖7所示，只要重疊區(qū)域不變，不管相互接觸三角形單元的另外一條邊如何，或者說三角形單元的大小、整體形狀如何，計算出的接觸力大小都應(yīng)該相等，本文提出的基于統(tǒng)一標(biāo)定的勢接觸力計算方法，對圖7所示的情況，計算的總接觸力是完全相同的，因為重疊區(qū)域的勢只與該點到三角形單元的兩條邊的距離有關(guān)，因此，重疊區(qū)域的勢是固定的，因而計算出的總接觸力完全一致。這說明新定義的勢函數(shù)，計算出的接觸力具有很好的局部特性，勢基于統(tǒng)一標(biāo)準(zhǔn)對嵌入量進行了表征，解決了原有勢函數(shù)定義存在的問題，新的勢函數(shù)定義使得基于勢的接觸力計算有了堅實的物理基礎(chǔ)，意義重大。

圖6 嵌入量相同的4組接觸Fig.6 Four group contacts with same embed amount

圖7 相同嵌入量的一般化Fig.7 Generalization of same embed amount

為了進一步驗證上述論述，對圖6所示的4種接觸情況，即ABE 與CDF 接觸對、ABE 與CD F′ 接觸對，即A′ BE 與CDF 接觸對、A′ BE 與CD F′ 接觸對，對這4組接觸對的總接觸力采用原有勢接觸力計算方法和本文的方法進行計算，由于三角形單元在變化，分配到各個節(jié)點上的力也在變化，為此僅僅比較作用在目標(biāo)單元上總接觸力，計算結(jié)果如表3所示。由表可知，原有勢函數(shù)定義計算的接觸力大小不等，而新定義的勢函數(shù)計算的接觸力完全一致。

表3 基于新、舊勢函數(shù)計算的接觸力比較Table 3 Comparison of contact forces between new and old potential functions

4.3 算例3

算例1和算例2均是對兩種勢接觸力計算方法的驗證和對比，此處給出一個應(yīng)用性的算例，來充分展現(xiàn)基于統(tǒng)一標(biāo)定的勢接觸力計算法的優(yōu)點，由于采用統(tǒng)一標(biāo)定，所有單元的接觸力計算都基于同一標(biāo)準(zhǔn)，這樣便不會出現(xiàn)局部某個接觸對單元的接觸力偏大或者偏小，而導(dǎo)致系統(tǒng)算出的接觸力分布不均勻。由于采用統(tǒng)一標(biāo)定的勢函數(shù)定義，不會因為單元形狀和尺寸的差異而導(dǎo)致接觸力計算相差懸殊，對求解域的網(wǎng)格劃分不再有嚴(yán)格的要求，而對原有方法則一般要求采用尺寸幾乎相同的網(wǎng)格，否則因為尺寸的差異，導(dǎo)致局部即便嵌入相等，但接觸力大小不一致。同時，單元的尺寸相差較大時，對小尺寸單元，相當(dāng)于變相增加其剛度值，使得剛度參數(shù)失去了原有的物理意義。實際上，只要材料確定了，其剛度是不會改變的，當(dāng)然不應(yīng)隨單元的尺寸和形狀而改變。而新的勢函數(shù)則不存在這一問題，剛度參數(shù)為材料的本質(zhì)屬性，不隨單元尺寸和形狀而改變。采用一個巴西圓盤劈裂的算例來加以說明，巴西圓盤的網(wǎng)格如圖8所示，其網(wǎng)格尺寸和形狀并不均一，用兩種勢接觸力計算方法對該算例進行計算，由于兩種方法采用不同的勢函數(shù)定義，為了使得程序有相同參數(shù)輸入，兩者需要采用不同的法向罰參數(shù)，亦即采用原有程序計算時，剛度參數(shù)需要按照系統(tǒng)內(nèi)所有三角形單元的高的平均值進行折算，已知所劃分的三角形單元高平均值為0.24 m，原有程序輸入的法向罰參數(shù)值，使得兩者采用的罰參數(shù)不一致，是由于原有勢函數(shù)定義方法存在變相增加剛度的問題，這樣處理使得兩種方法輸入的實際等效剛度接近，其余的力學(xué)參數(shù)完全一致，只有這樣才能對兩者的計算結(jié)果進行比較。這種方法的模擬結(jié)果如圖8～16所示。

由圖8～16的計算結(jié)果可知，從應(yīng)力分布的均勻程度上看，新方法獲得應(yīng)力分布更為均勻，應(yīng)力云圖中不同顏色的交界面更為光滑，同時在同一種顏色區(qū)域內(nèi)部較少出現(xiàn)其他顏色的空洞，亦即代表應(yīng)力值的顏色分布較均勻，如圖9、10中的時步為10 000、15 000。從裂紋的擴展形態(tài)上看，新方法獲得裂紋形態(tài)基本沿著圓盤的豎向直徑擴展，與理論分析較為接近；而原有勢接觸力計算方法獲得的裂紋出現(xiàn)了分支，圓盤下端出現(xiàn)了偏離豎向直徑的裂紋，如圖14、15中的時步為35 000、40 000所示。從兩種方法計算得到荷載-位移曲線上看（如圖16所示），兩種方法獲得的峰值強度基本一致，但新方法獲得的荷載-位移曲線更為平滑，震蕩幅度較小，這一點在曲線的峰后部分表現(xiàn)得尤為明顯，這說明基于統(tǒng)一標(biāo)定的勢接觸力計算方法要比原有勢接觸力方法計算更穩(wěn)定。綜上所述可知，基于新的勢函數(shù)定義獲得的應(yīng)力分布也更均勻，較少出現(xiàn)應(yīng)力集中，裂紋擴展形態(tài)與理論解更為接近，同時計算過程更為穩(wěn)定，對網(wǎng)格劃分的魯棒性更高。

圖8 兩種方法計算結(jié)果的比較（時步為5 000）Fig.8 Comparison of two methods of calculation (step is 5 000)

圖9 兩種方法計算結(jié)果的比較（時步為10 000）Fig.9 Comparison of two methods of calculation (step is 10 000)

圖10 兩種方法計算結(jié)果的比較（時步為15 000）Fig.10 Comparison of two methods of calculation (step is 15 000)

圖11 兩種方法計算結(jié)果的比較（時步為20 000）Fig.11 Comparison of two methods of calculation (step is 20 000)

圖12 兩種方法計算結(jié)果的比較（時步為25 000）Fig.12 Comparison of two methods of calculation (step is 25 000)

圖13 兩種方法計算結(jié)果的比較（時步為30 000）Fig.13 Comparison of two methods of calculation (step is 30 000)

圖14 兩種方法計算結(jié)果的比較（時步為35 000）Fig.14 Comparison of two methods of calculation (step is 35 000)

圖15 兩種方法計算結(jié)果的比較（時步為40 000）Fig.15 Comparison of two methods of calculation (step is 40 000)

圖16 巴西圓盤劈裂試驗的荷載-位移曲線Fig.16 Load-displacement curves of the Brazilian disc tests

5 結(jié) 論

（1）指出了原有勢函數(shù)定義，存在物理意義不明確，計算出的接觸力與物理直觀不符，同時勢無法對嵌入量進行統(tǒng)一表征，接觸力的計算沒有建立在一個統(tǒng)一的標(biāo)準(zhǔn)之上等問題。

（2）提出了一種新的勢函數(shù)定義方法，三角形單元內(nèi)一點的勢與該點到三角形單元3條邊的最短距離成正比，三角單元內(nèi)位于內(nèi)心點處的勢最大，但并不一定為1，且每個三角形內(nèi)位于內(nèi)心點處的勢的值不一定相同。

（3）新提出的勢函數(shù)，對嵌入量進行了統(tǒng)一的表征，物理意義明確、直觀，使得基于勢的接觸力計算建立在堅實的物理基礎(chǔ)之上，意義重大。

（4）新定義的勢函數(shù)，解決原有勢函數(shù)所存在的重大問題，只要重疊區(qū)域不變，計算的總的接觸力將是不變的，接觸力計算具有局部特性，而與劃分的三角形單元的整體形狀和大小無關(guān)，即新的接觸力計算方法，具有更好的網(wǎng)格魯棒性。

（5）基于統(tǒng)一標(biāo)定的勢接觸力計算方法，計算出的接觸力大小分布更為均勻，不會在局部產(chǎn)生過大的接觸力，使得所有單元的接觸力計算都基于一個統(tǒng)一的標(biāo)準(zhǔn)，不容易產(chǎn)生應(yīng)力集中，計算更為穩(wěn)定。

[1]CUNDALL P A.A computer model for simulating progressive large scale movement in block rock system[C]//Nancy:Symposium ISRM，1971，2:129－136.

[2]石根華.數(shù)值流形方法與非連續(xù)變形分析[M].裴覺民譯.北京:清華大學(xué)出版社，1997:192－194.

[3]鄔愛清，丁秀麗，盧波，等.DDA方法塊體穩(wěn)定性驗證及其在巖質(zhì)邊坡穩(wěn)定性分析中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報，2008，27(4):664－664.WU Ai-qing，DING Xiu-li，LU Bo，et al.Validation for rock block stability and its application to rock slope stability evaluation using DDA method[J].Chinese Journal of Rock Mechanics and Engineering，2008，27(4):664－664.

[4]鄭宏，江巍.基于互補理論的非連續(xù)變形分析方法[J].中國科學(xué)(E輯:技術(shù)科學(xué))，2009，39(10):1702－1708.ZHENG Hong，JIANG Wei.Discontinuous deformation analysis based on complementary theory[J].Science China (Series E:Technological Sciences)，2009，39(10):1702－1708.

[5]徐棟棟，孫冠華，鄭宏.接觸處理Munjiza方法在數(shù)值流形法中的應(yīng)用[J].巖土力學(xué)，2013，34(2):526－534.XU Dong-dong，SUN Guan-hua，ZHENG Hong.Application of Munjiza’s algorithm to contact treatment in numerical manifold method[J].Rock and Soil Mechanics，2013，34(2):526－534.

[6]Itasca Consulting Group Ltd.User manual of UDEC code[S].Minneapolis:Itasca Consulting Group Ltd.，2011.

[7]陳文勝，王桂堯，劉輝，等.巖石力學(xué)離散單元計算方法中的若干問題探討[J].巖石力學(xué)與工程學(xué)報，2005，24(10):1639－1644.CHEN Wen-sheng，WANG Gui-yao，LIU Hui，et al.Insight into some aspects of discrete element numerical methods for rock mass[J].Chinese Journal of Rock Mechanics and Engineering，2005，24(10):1639－1644.

[8]CUNDALL P A.Formulation of a three-dimensional distinct element model-part I:A scheme to detect and represent contact in composed of many polyhedral blocks[J].International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts，1988，25(3):107－116.

[9]MUNJIZA A，ANDREWS K R F.Penalty function method for combined finite-discrete element systems comprising large number of separate bodies[J].International Journal for Numerical Methods in Engineering，2000，49(11):1377－1369.

[10]MUNJIZA A.The combined finite-discrete element method[M].London:John Wiley &Sons，Ltd.，2004:29－32.

[11]嚴(yán)成增，孫冠華，鄭宏，等.基于局部單元劈裂的FEM/DEM自適應(yīng)分析方法[J].巖土力學(xué)，2014，35(7):2064－2070.YAN Cheng-zeng，SUN Guan-hua，ZHENG Hong，et al.Adaptive FEM/DEM analysis method based on the local splitting elements[J].Rock and Soil Mechanics，2014，35(7):2064－2070.

[12]嚴(yán)成增，孫冠華，鄭宏，等.三維FEM/DEM中摩擦力的實施及驗證[J].巖石力學(xué)與工程學(xué)報，2014，33(6):1248－1256.YAN Cheng-zeng，SUN Guan-hua，ZHENG Hong，et al.Implementation of friction to 3D FEM/DEM and verification[J].Chinese Journal of Rock Mechanics and Engineering，2014，33(6):1248－1256.