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

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

1 引 言

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

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

2 FDEM勢(shì)接觸力

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

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

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

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

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

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

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

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

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

式中：P為靶單元中長度為L 的邊上的任意一點(diǎn)。

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

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

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

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

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

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

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

勢(shì)接觸力為

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

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

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

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

3.1 基本思路

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

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

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

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

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

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

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

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

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

3.2 實(shí)施方案

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

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

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

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

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

4 算 例

4.1 算例1

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

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

4.2 算例2

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

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

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

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

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

4.3 算例3

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

（5）基于統(tǒng)一標(biāo)定的勢(shì)接觸力計(jì)算方法，計(jì)算出的接觸力大小分布更為均勻，不會(huì)在局部產(chǎn)生過大的接觸力，使得所有單元的接觸力計(jì)算都基于一個(gè)統(tǒng)一的標(biāo)準(zhǔn)，不容易產(chǎn)生應(yīng)力集中，計(jì)算更為穩(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)定性驗(yàn)證及其在巖質(zhì)邊坡穩(wěn)定性分析中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào)，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]鄭宏，江巍.基于互補(bǔ)理論的非連續(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ì)算方法中的若干問題探討[J].巖石力學(xué)與工程學(xué)報(bào)，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中摩擦力的實(shí)施及驗(yàn)證[J].巖石力學(xué)與工程學(xué)報(bào)，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.