潘洪武，王 偉，張丙印

(清華大學水沙科學與水利水電重點實驗室，北京 100084)

顆粒物質(zhì)是離散和相互作用的大量固體顆粒所組成的復雜體系。顆粒物質(zhì)廣泛存在于自然界，與人類的日常生活和生產(chǎn)密切相關(guān)。但是，顆粒體系的很多靜力學和動力學行為尚不能用一般的連續(xù)介質(zhì)力學理論很好地解釋?！禨cience》雜志于2005 年把顆粒物質(zhì)與湍流并列為100 個科學難題之一。

土體由離散的顆粒組成。顆粒間以強耗散的接觸摩擦為主，其宏觀力學特性主要取決于細觀層次的顆粒聯(lián)結(jié)和排列等。近年來，人們逐漸開始重視土顆粒體系細觀力學行為的研究。由于在物理試驗中存在對顆粒細觀行為測量技術(shù)的難題，因而各種數(shù)值模擬方法逐步成為主要的研究手段。對此，學者們已經(jīng)發(fā)展出了離散元法、非連續(xù)變形方法和連續(xù)-離散耦合分析方法等多種數(shù)值計算方法。

離散元法是Cundall[1]于1971 年提出的一種基于牛頓第二運動定律的非連續(xù)變形數(shù)值方法。該法將土體看成是離散單元的集合，通過給定顆粒間的接觸關(guān)系[2]，來反映顆粒材料的宏細觀力學響應(yīng)。離散元法具有很強的描述顆粒間不連續(xù)變形的能力，受到人們廣泛的青睞和重視，目前是進行顆粒材料細觀力學特性分析的最主要的方法，取得了豐碩的研究成果[3-8]。

由于計算方法的限制，離散元法對顆粒本身細觀力學特性的描述(如顆粒內(nèi)部應(yīng)力分布等)相對較為粗糙，其細觀數(shù)值模擬難以定量地再現(xiàn)顆粒材料的宏細觀力學特性。為此，Munjiza、馬剛、王永亮等[9-11]發(fā)展了連續(xù)-離散耦合分析方法(FDEM)。

有限元方法是目前應(yīng)用最廣泛的數(shù)值模擬計算方法。但一般認為，有限元法主要適合于連續(xù)介質(zhì)的模擬計算分析，其在不連續(xù)變形的描述方面存在缺陷。近年來，基于有限元方法的計算接觸力學算法[12-15]得到了迅速發(fā)展，使得有限元方法在模擬不連續(xù)變形和多體接觸問題方面的能力大幅提高。

土體是由離散的土顆粒組成的，其細觀力學性質(zhì)的研究為一個典型的多體接觸問題。從原理上講，計算接觸力學方法是適用于進行這種多體接觸分析的。基于上述的思路，本文開展了基于計算接觸力學的土顆粒材料細觀力學特性模擬的探索性研究工作。通過對粗顆粒土體材料進行三軸試驗的數(shù)值計算分析，探討了基于有限元方法的計算接觸力學算法對粗顆粒土體材料多體接觸問題的適用性。

1 計算接觸力學算法簡介

計算接觸力學方法來自固體力學領(lǐng)域，是非線性數(shù)值分析中最具挑戰(zhàn)性的課題之一，提出以來得到了來自計算數(shù)學和計算力學等領(lǐng)域諸多學者的研究。計算接觸力學研究內(nèi)容繁多，涉及范圍很廣。主要研究的問題包括接觸界面離散形式(點對點、點對面和面對面等)、接觸非線性求解、接觸本構(gòu)方程和接觸約束施加方法等。計算接觸力學方法被廣泛應(yīng)用在機械、航空和生物力學等領(lǐng)域。在巖土工程領(lǐng)域近年來才開始得到應(yīng)用[16-18]。

1.1 接觸問題的一般描述

直觀地講，接觸問題反映的是不同物體之間的相互作用。對于本文考慮的粗顆粒土體材料，不失一般性，考慮如圖1 所示的可變形顆粒 Ω1和Ω2的接觸問題。將顆粒 Ωi的邊界記作 ?Ωi(i=1,2)，每個邊界可分為Dirichlet 邊界、Neumann 邊界和接觸邊界，三者組成整體邊界且相互無重疊：

圖1 兩個顆粒的接觸問題示意圖Fig. 1 Schematic diagram of two particle contact problem

接觸邊界上，兩物體被記作從面和主面。將物體表面的法向和切向單位向量分別記為 n、 τ1和τ2， 接觸區(qū)域的法向間隙 gn可表示為：

式中： xA為從面點的坐標； x?B為 點 xA在Γ上的投影點。接觸區(qū)域的物體接觸力 tC可記作：

式中： pn為 接觸力法向分量；tτ1、為接觸力的切向分量。

在計算接觸力學中，一般將上述接觸條件表

1.2 接觸問題的弱形式

1.3 接觸界面的空間離散

為了在接觸界面上計算接觸虛功，需要對接觸界面進行離散，常用的界面離散方法有點點方法、點面方法和面面方法等，其中，面面方法具有計算精度高、能通過接觸分片試驗等優(yōu)點，是國內(nèi)外學者們目前研究的熱點，也是本文數(shù)值研究采用的方法。

在接觸界面上，引入拉格朗日乘子，通過節(jié)點量插值記為：

1.4 接觸問題的非線性迭代

由于接觸區(qū)域和接觸狀態(tài)的雙重不確定性，接觸問題具有極強的非線性。計算接觸力學中常采用PDASS (primal dual active-set strategy)方法[21]對接觸問題進行非線性迭代求解。其基本思路為：首先將接觸界面上的從點劃分為脫開、滑移和粘結(jié)等不同狀態(tài)的點集，并分別引入相應(yīng)的接觸邊界條件；然后再根據(jù)相應(yīng)互補的接觸條件進行校核，判別從點的接觸狀態(tài)是否正確；對接觸狀態(tài)不正確的從點重新劃分其接觸狀態(tài)。如此迭代直至接觸界面從點的狀態(tài)不再發(fā)生改變?yōu)橹埂?/p>

1.5 接觸問題有限元方程的代數(shù)表達式

本文在界面離散上采用Dual mortar 元方法，其得到的D 矩陣為對角矩陣，可以高效地進行求逆操作。在激活節(jié)點上引入相對位移：

1.6 計算程序和驗證算例

基于課題組發(fā)展的多體接觸有限元計算程序Contana，針對顆粒集合體的特點，發(fā)展了相應(yīng)的土體顆粒多體接觸有限元計算程序。

采用四球堆積算例對所發(fā)展的有限元程序進行驗證。如圖2 所示，四個球體的半徑R 相等。球體B2、B3和B4位于平面W1上，其球心連線為邊長等于l=0.105 m 的等邊三角形。球體B1對稱放置在三個球體的頂部，四個球體的球心連線形成一個四面體。球體半徑R=0.05 m，密度ρ=7800 kg/m3，彈性模型E=100 GPa，泊松比ν=0.3。本算例中，采用六面體單元對球體進行建模，每個球體含有1734 個節(jié)點、1152 個單元。

當球體間以及球體和底面間的摩擦系數(shù)μ1和μ2足夠大時，四個球體處于穩(wěn)定狀態(tài)。記球體之間以及球體和底面之間接觸力的法向和切向分量分別為Fn1、Ft1、Fn2和Ft2。根據(jù)靜力平衡條件，可以求得不考慮球體變形條件下各接觸力的大小。由于在算例情況下球體變形非常小，下文稱該解答為理論解。

表1 給出了有限元接觸力計算結(jié)果與理論值的對比情況。其中，理論解扣除了有限元網(wǎng)格劃分所引起的質(zhì)量誤差。從表1 中可以看出，接觸力的計算值和相應(yīng)理論解的最大相對誤差為0.129%，數(shù)值解和理論解吻合較好。

圖2 四球堆積算例Fig. 2 Four sphere accumulation example

表1 接觸力計算結(jié)果的對比Table 1 Comparison between calculation and theoretical values

2 三軸數(shù)值試驗的計算模型

利用所發(fā)展的土體顆粒多體接觸有限元計算程序，進行了不同圍壓的常規(guī)三軸數(shù)值試驗，探討了計算接觸力學方法對土體顆粒系統(tǒng)細觀力學特性分析的適用性。

2.1 計算模型概述

三軸試樣由粒徑分別為21 mm、19 mm、17 mm、15 mm 和13 mm 的460 個球顆粒組成。上述不同粒徑的顆粒數(shù)目分別為70、95、130、95 和70。

如圖3 所示，計算采用10 cm×10 cm×20 cm的柱體試樣。試樣的四個側(cè)面和兩個頂面共設(shè)置6 塊加載板。加載板使用六面體單元離散。在加載板之間不設(shè)任何的接觸關(guān)系，它們在計算中不會發(fā)生相互作用，可獨立施加所需的應(yīng)力或位移。每個球顆粒均被剖分為510 個四面體單元、133個節(jié)點。整個計算模型共包含有78748 個節(jié)點、245564 個單元。

圖3 計算模型示意圖Fig. 3 Schematic diagram of simulation model

計算中，顆粒和加載板均采用線彈性模型模擬，彈性模量分別取E=1 GPa 和25 GPa，泊松比均取ν=0.3。不考慮重力影響。顆粒-顆粒之間以及顆粒-加載板之間均采用庫侖摩擦定律模擬，摩擦系數(shù)分別取0.4 和0.2。

2.2 數(shù)值試驗計算過程

為了能與土體三軸試驗的結(jié)果進行定性對比，參照土工三軸物理試驗的流程，將本文數(shù)值試驗的計算流程也劃分為制樣、固結(jié)和加載三個過程。數(shù)值試驗的試樣參照目前離散元方法所采用的三軸試樣生成方法進行生成[5]。

采用生成的同一個數(shù)值試驗試樣，分別進行了圍壓σ3=200 kPa、400 kPa 和600 kPa 的常規(guī)三軸數(shù)值壓縮試驗。計算時，首先模擬施加圍壓的固結(jié)過程。采用Δσ3=50 kPa 分級施加圍壓力σ3。

對每一個圍壓的試驗，固結(jié)完成后，保持側(cè)向圍壓力σ3大小不變，采用應(yīng)變控制的加載方式進行剪切加載。具體方式為，每個加載步使加載頂板向下位移0.01 cm，對應(yīng)試樣發(fā)生0.05%的軸向應(yīng)變。

計算在個人計算機在上進行。計算機CPU 的型號為Intel? CoreTMi7-2600K(3.40 GHz)，每個圍壓試驗的計算耗時約為40 h。

3 三軸數(shù)值試驗計算結(jié)果分析

在本文的計算結(jié)果中，應(yīng)力和應(yīng)變均以壓為正。如前所述，采用計算接觸力學的方法，可以對顆粒本身劃分有限元網(wǎng)格。因此，所得數(shù)值試驗的結(jié)果不僅能反映試樣的整體宏觀力學特性，并且也能反映顆粒本身的細部力學性質(zhì)。

3.1 試樣的宏觀力學特性

圖4 給出了σ3=600 kPa 數(shù)值試驗過程中試樣的側(cè)視圖。可以看出，在試驗過程中，顆粒本身不發(fā)生大的變形，試樣的變形主要由顆粒在空間相對位置的變化所致，即主要由顆粒間的滑移和顆粒轉(zhuǎn)動等所導致的顆粒集合體排列方式的改變所產(chǎn)生。

圖4 試驗過程中的試樣側(cè)視圖(σ3=600 kPa)Fig. 4 Side view of specimen during experiment(σ3=600 kPa)

圖5 給出了不同圍壓的應(yīng)力-應(yīng)變曲線和體變計算曲線。由圖5(a)所示的(σ1-σ3)-ε1曲線可以看到，各圍壓的應(yīng)力-應(yīng)變曲線均表現(xiàn)出一定的應(yīng)變軟化特征。偏差應(yīng)力隨軸向應(yīng)變的增加逐步增大，在開始階段增加相對較快，之后逐步變得較為平緩，并達到峰值。峰值后隨軸向應(yīng)變增加，偏差應(yīng)力會發(fā)生一定的下降。對不同圍壓的試驗，對應(yīng)相同的軸向應(yīng)變，圍壓大對應(yīng)的偏差應(yīng)力也越高，表明試樣的剛度或強度隨圍壓逐步增大，表現(xiàn)出了顆粒材料壓硬性的特點。這些特點均符合土體常規(guī)三軸壓縮試驗應(yīng)力-應(yīng)變曲線的一般特征。

圖5 不同圍壓數(shù)值試驗的計算結(jié)果Fig. 5 Computation results of numerical experiment at different confining pressures

由圖5(b)所示的εv～ε1曲線可以看出，所得體變曲線總體符合粗粒土三軸試驗中的體變特性。對一個圍壓的數(shù)值試驗，試樣在初始先發(fā)生體縮。當軸向應(yīng)變較大時，體積變形逐步轉(zhuǎn)變?yōu)榧裘洝Σ煌瑖鷫旱臄?shù)值試驗，圍壓較小時，試驗初始段的體縮較小，后段的剪脹較為顯著。隨著圍壓的提高，試樣的剪縮性總體增大，剪脹性逐步減小。根據(jù)試驗結(jié)果繪制莫爾圓可得試樣內(nèi)摩擦角約為44°。

3.2 顆粒細部應(yīng)力分析

圖6 給出圍壓σ3=600 kPa 試驗試樣內(nèi)部的力鏈分布情況。由圖6 可見，在試驗過程中，水平方向接觸力小于豎直方向接觸力，荷載主要沿豎直方向傳遞。對比不同軸向應(yīng)變下試樣力鏈分布可知，顆粒間接觸力隨著軸向應(yīng)變增加而增大。

圖6 試樣力鏈分布Fig. 6 Distribution of the force chain

圖7 進一步給出了圍壓σ3=600 kPa 試驗中顆粒間法向接觸力大小的分布情況。由圖7 可以看出，當試樣固結(jié)完成尚未進行剪切時，顆粒間法向接觸力總體相對較小，絕大多數(shù)分布在小于4 kN以下的范圍。開始加載剪切后，顆粒間法向接觸力逐步增大。隨軸向應(yīng)變的增大，小法向力接觸對數(shù)目逐漸降低，大法向力接觸對數(shù)目逐漸增加。當軸向應(yīng)變達到14%時，一些顆粒的接觸應(yīng)力可達16 kN 以上。

如前所述，采用計算接觸力學方法可以對顆粒本身劃分所需的有限元的單元，從而可以計算得到其詳細的細觀應(yīng)力和變形狀態(tài)。這是計算接觸力學方法相比離散元方法的一個優(yōu)勢。在本次計算中，每個球顆粒均包含133 個節(jié)點、510 個單元，因而可以得到任意時刻所有顆粒詳細的應(yīng)力狀態(tài)。

本文的計算模型中，顆粒采用四面體常應(yīng)變單元進行網(wǎng)格剖分，且單元尺寸相近，因此可根據(jù)顆粒內(nèi)部單元應(yīng)力的統(tǒng)計分布來評估顆粒的總體應(yīng)力狀態(tài)。圖8 給出了σ3=600 kPa 圍壓數(shù)值試驗所得某個顆粒中單元大主應(yīng)力大小的分布情況。當施加圍壓力之后但尚未進行剪切之前，試樣處于等向壓縮應(yīng)力狀態(tài)。此時，假如試樣是連續(xù)的均質(zhì)材料，則試樣中所有單元的三個主應(yīng)力都應(yīng)為600 kPa。但由于試樣由土體顆粒組成，因而計算所得單元應(yīng)力的大小并不相等。其中，少數(shù)單元處于受拉狀態(tài)，計算所得大主應(yīng)力的最大值約為1.5 MPa。

圖7 顆粒間法向接觸力大小統(tǒng)計分布Fig. 7 Distribution of normal contact forces between particles

圖8 顆粒單元大主應(yīng)力統(tǒng)計分布Fig. 8 Distribution of large principle stress in a particle

當進行剪切后，由圖8 可見，計算所得顆粒單元的大主應(yīng)力總體增大，承受高應(yīng)力的單元數(shù)目顯著增加。少數(shù)單元仍會處于受拉狀態(tài)。當軸向應(yīng)變達到14%時，一些顆粒單元計算所得大主應(yīng)力的最大值增大約為3.0 MPa。

圖9 為σ3=600 kPa 圍壓的三軸數(shù)值試驗，當軸向應(yīng)變ε1=10%時，某典型顆粒中的大主應(yīng)力分布的計算結(jié)果。由圖9 可以看出，計算所得單元的大主應(yīng)力在顆粒中的分布非常的不均勻。其中，在顆粒中的大部分區(qū)域，計算所得到的應(yīng)力值相對均較低。但在顆粒接觸點附近的局部小區(qū)域，則會計算得到相對很高的應(yīng)力分布。對于顆粒材料，其所受到的作用力是通過顆粒接觸點進行傳遞的。對于本文所采用的球形顆粒，在接觸點部位會作用集中的粒間力的作用。計算所得上述的應(yīng)力分布特征符合赫茲接觸問題的一般規(guī)律。

圖9 大主應(yīng)力剖視圖 /MPa Fig. 9 Section view of large principal stress

4 計算接觸力學法和離散元法對比

由本文前面的分析可知，計算接觸力學方法是基于有限元方法的一種進行多體接觸問題的數(shù)值求解方法。與離散元方法相比，在基本原理、求解方法、顆粒本身的模擬以及接觸條件的模擬等方面都存在顯著的不同。表2 總結(jié)列出了兩種方法在模擬顆粒材料特性方面的主要差別。

表2 兩種模擬方法對比Table 2 Comparison between two simulation methods

總體而言，將計算接觸力學方法引入土體顆粒集合體細觀應(yīng)力變形特性的模擬分析，可能會具有如下的優(yōu)勢：

1)在基本原理和求解方面，計算接觸力學方法是一種基于變分原理的隱式數(shù)值求解方法，具有嚴格的理論基礎(chǔ)。離散元方法的基本原理是牛頓第二定律，根據(jù)顆粒受到的合力計算顆粒在各時間步的位移，是一種顯示求解方法。自提出至今，離散元及其相關(guān)方法的研究已有近50 年歷史，眾多學者發(fā)表了大量學術(shù)論文。但是，離散元方法缺乏理論的嚴密性。劉凱欣和高凌天[22]曾評述“研究離散元的理論和算法的文章卻很少”“自誕生的那天起就帶有缺乏理論嚴密性的先天不足”。

2)在顆粒本身的模擬方面，采用計算接觸力學方法可以對顆粒本身劃分所需的有限單元，從而可以得到其詳細的細觀應(yīng)力和變形狀態(tài)，可為顆粒破碎等的模擬計算提供基礎(chǔ)。傳統(tǒng)離散元方法的模型一般由剛性塊體組成，此時無法考慮顆粒體本身變形，也無法計算得到顆粒體內(nèi)部應(yīng)力的分布。

3)在顆粒間的接觸特性模擬方面，計算接觸力學方法通過拉格朗日乘子法引入KKT 接觸條件，來模擬顆粒接觸界面的接觸特性，是一種高精度的接觸問題求解方法。拉格朗日乘子的物理意義為接觸應(yīng)力。該方法無須引入接觸彈簧，可避免罰模量選取對人工經(jīng)驗的依賴等。離散元方法采用罰函數(shù)法，通過在接觸點引入法向和切向彈簧模擬顆粒間的接觸行為，存在彈簧剛度以及罰模量選取對人工經(jīng)驗的依賴等。

但是，大規(guī)模的土體顆粒集合體是一個大型強非連續(xù)顆粒系統(tǒng)，將海量的顆粒間接觸條件引入有限元方程后，會造成系統(tǒng)剛度矩陣的高度病態(tài)。為此，研發(fā)針對大型復雜顆粒集合體的高效接觸搜索算法、高度病態(tài)矩陣預處理技術(shù)和具有較強魯棒性的迭代算法等，是該種方法能否成功應(yīng)用于土體顆粒集合體細觀力學行為模擬分析的關(guān)鍵。這也是該法亟待在未來研究中解決的關(guān)鍵技術(shù)難題。

5 結(jié)論

本文開展了基于計算接觸力學的粗顆粒土體材料細觀力學行為模擬分析的探索性研究工作，主要取得如下的結(jié)論：

(1)提出了一種基于計算接觸力學的粗顆粒土體材料細觀力學行為模擬分析新方法。該法將顆粒剖分成一定數(shù)量的有限元單元，通過計算接觸力學方法模擬顆粒間的接觸行為。

(2)和離散元法方法相比，本文提出的模擬計算方法是一種基于變分原理的隱式數(shù)值求解方法，具有嚴格的理論基礎(chǔ)。此外，該方法在描述顆粒本身的力學特性方面具有優(yōu)勢，既可以計算分析粗顆粒土體材料的宏細觀力學特性，也可以計算得到顆粒內(nèi)部的應(yīng)力分布情況，可為顆粒破碎等的細觀行為的模擬計算提供依據(jù)。

(3)利用所發(fā)展的土體顆粒多體接觸有限元計算程序，進行了不同圍壓的常規(guī)三軸數(shù)值試驗，模擬計算結(jié)果符合土體三軸試驗的一般規(guī)律，初步驗證了所提出的計算方法的適用性。