基于計算接觸力學的粗顆粒土體材料細觀性質模擬

潘洪武，王 偉，張丙印

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

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

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

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

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

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

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

1 計算接觸力學算法簡介

計算接觸力學方法來自固體力學領域，是非線性數值分析中最具挑戰性的課題之一，提出以來得到了來自計算數學和計算力學等領域諸多學者的研究。計算接觸力學研究內容繁多，涉及范圍很廣。主要研究的問題包括接觸界面離散形式(點對點、點對面和面對面等)、接觸非線性求解、接觸本構方程和接觸約束施加方法等。計算接觸力學方法被廣泛應用在機械、航空和生物力學等領域。在巖土工程領域近年來才開始得到應用[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， 接觸區域的法向間隙 gn可表示為：

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

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

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

1.2 接觸問題的弱形式

1.3 接觸界面的空間離散

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

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

1.4 接觸問題的非線性迭代

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

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

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

1.6 計算程序和驗證算例

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

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

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

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

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

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

2 三軸數值試驗的計算模型

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

2.1 計算模型概述

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

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

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

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

2.2 數值試驗計算過程

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

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

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

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

3 三軸數值試驗計算結果分析

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

3.1 試樣的宏觀力學特性

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

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

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

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

由圖5(b)所示的εv～ε1曲線可以看出，所得體變曲線總體符合粗粒土三軸試驗中的體變特性。對一個圍壓的數值試驗，試樣在初始先發生體縮。當軸向應變較大時，體積變形逐步轉變為剪脹。對不同圍壓的數值試驗，圍壓較小時，試驗初始段的體縮較小，后段的剪脹較為顯著。隨著圍壓的提高，試樣的剪縮性總體增大，剪脹性逐步減小。根據試驗結果繪制莫爾圓可得試樣內摩擦角約為44°。

3.2 顆粒細部應力分析

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

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

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

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

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

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

圖8 顆粒單元大主應力統計分布Fig. 8 Distribution of large principle stress in a particle

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

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

圖9 大主應力剖視圖 /MPa Fig. 9 Section view of large principal stress

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

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

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

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

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

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

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

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

5 結論

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

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

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

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