潘洪武，王 偉，張丙印

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

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

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

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

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

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

土體是由離散的土顆粒組成的，其細(xì)觀力學(xué)性質(zhì)的研究為一個典型的多體接觸問題。從原理上講，計算接觸力學(xué)方法是適用于進(jìn)行這種多體接觸分析的?；谏鲜龅乃悸罚疚拈_展了基于計算接觸力學(xué)的土顆粒材料細(xì)觀力學(xué)特性模擬的探索性研究工作。通過對粗顆粒土體材料進(jìn)行三軸試驗(yàn)的數(shù)值計算分析，探討了基于有限元方法的計算接觸力學(xué)算法對粗顆粒土體材料多體接觸問題的適用性。

1 計算接觸力學(xué)算法簡介

計算接觸力學(xué)方法來自固體力學(xué)領(lǐng)域，是非線性數(shù)值分析中最具挑戰(zhàn)性的課題之一，提出以來得到了來自計算數(shù)學(xué)和計算力學(xué)等領(lǐng)域諸多學(xué)者的研究。計算接觸力學(xué)研究內(nèi)容繁多，涉及范圍很廣。主要研究的問題包括接觸界面離散形式(點(diǎn)對點(diǎn)、點(diǎn)對面和面對面等)、接觸非線性求解、接觸本構(gòu)方程和接觸約束施加方法等。計算接觸力學(xué)方法被廣泛應(yīng)用在機(jī)械、航空和生物力學(xué)等領(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為從面點(diǎn)的坐標(biāo)； x?B為 點(diǎn) xA在Γ上的投影點(diǎn)。接觸區(qū)域的物體接觸力 tC可記作：

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

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

1.2 接觸問題的弱形式

1.3 接觸界面的空間離散

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

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

1.4 接觸問題的非線性迭代

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

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

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

1.6 計算程序和驗(yàn)證算例

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

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

當(dāng)球體間以及球體和底面間的摩擦系數(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ù)值試驗(yàn)的計算模型

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

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ā)生相互作用，可獨(dú)立施加所需的應(yīng)力或位移。每個球顆粒均被剖分為510 個四面體單元、133個節(jié)點(diǎn)。整個計算模型共包含有78748 個節(jié)點(diǎn)、245564 個單元。

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

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

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

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

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

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

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

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

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

3.1 試樣的宏觀力學(xué)特性

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

圖4 試驗(yàn)過程中的試樣側(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)變的增加逐步增大，在開始階段增加相對較快，之后逐步變得較為平緩，并達(dá)到峰值。峰值后隨軸向應(yīng)變增加，偏差應(yīng)力會發(fā)生一定的下降。對不同圍壓的試驗(yàn)，對應(yīng)相同的軸向應(yīng)變，圍壓大對應(yīng)的偏差應(yīng)力也越高，表明試樣的剛度或強(qiáng)度隨圍壓逐步增大，表現(xiàn)出了顆粒材料壓硬性的特點(diǎn)。這些特點(diǎn)均符合土體常規(guī)三軸壓縮試驗(yàn)應(yīng)力-應(yīng)變曲線的一般特征。

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

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

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

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

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

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

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

本文的計算模型中，顆粒采用四面體常應(yīng)變單元進(jìn)行網(wǎng)格剖分，且單元尺寸相近，因此可根據(jù)顆粒內(nèi)部單元應(yīng)力的統(tǒng)計分布來評估顆粒的總體應(yīng)力狀態(tài)。圖8 給出了σ3=600 kPa 圍壓數(shù)值試驗(yàn)所得某個顆粒中單元大主應(yīng)力大小的分布情況。當(dāng)施加圍壓力之后但尚未進(jìn)行剪切之前，試樣處于等向壓縮應(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

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

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

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

4 計算接觸力學(xué)法和離散元法對比

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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