顆粒材料側(cè)向壓力分布特性離散元數(shù)值模擬分析

來志強(qiáng)，趙連軍，武彩萍，吳國英，張 磊

(1.黃河水利委員會黃河水利科學(xué)研究院 水利部黃河泥沙重點實驗室，鄭州 450003；2.中國科學(xué)院 地理科學(xué)與資源研究所，北京 100101；3.華北水利水電大學(xué) 水利學(xué)院，鄭州 450045)

0 引 言

水庫大壩工程中常見的顆粒材料有筑壩材料堆石體、土體、泥沙等。這些顆粒材料因自重或外荷載作用下會產(chǎn)生側(cè)向壓力，在高堆石壩、擋土墻、深埋隧道等工程建設(shè)時必須加以考慮。特別是有關(guān)土石側(cè)向壓力問題的研究已有200多年歷史，目前工程上常用的側(cè)向壓力計算理論是朗肯和庫倫土壓力理論。兩大經(jīng)典土壓力理論均假定側(cè)向壓力沿土石體深度呈線性分布。側(cè)向壓力線性分布的假定在國內(nèi)外擋土墻設(shè)計規(guī)范[1]、基坑結(jié)構(gòu)穩(wěn)定分析[2,3]及深埋隧道[4]等工程領(lǐng)域廣泛使用。然而土石等工程材料是非常復(fù)雜的顆粒材料, 傳統(tǒng)的連續(xù)介質(zhì)理論往往忽略了對其離散特性的考慮[5,6]?，F(xiàn)有研究表明土石等顆粒材料堆積后將導(dǎo)致拱結(jié)構(gòu)效應(yīng)的存在，這使得顆粒材料的側(cè)向壓力呈現(xiàn)非線性分布[7-9]，這與目前側(cè)向壓力線性分布假定相悖。因此，有必要研究顆粒材料側(cè)向壓力分布的非線性及其與經(jīng)典土壓力理論中線性分布假定之間的差異，對于目前工程設(shè)計具有十分重要的意義。

近年來，離散單元法(Discrete Element Method, DEM)的發(fā)展和計算機(jī)性能的飛速提高，使得通過細(xì)觀數(shù)值方法研究土石等顆粒材料的力學(xué)行為特性成為可能[10-14]。離散單元法模擬土石等顆粒材料時不需要定義宏觀本構(gòu)關(guān)系，其復(fù)雜的力學(xué)響應(yīng)將由其內(nèi)部顆粒間相互作用實現(xiàn)，同時不受連續(xù)介質(zhì)理論中連續(xù)性和小變形假定的約束。Jiang等[12]采用DEM分析了不同墻體運動模式下剛性擋墻側(cè)向壓力分布規(guī)律；孟云偉等[15]研究了堆石擋墻位移、應(yīng)力分布和破壞形式等力學(xué)行為；劉冬等[16]和萬勵等[17]分別探討了無黏性土體在平動模式下的破壞形式以及土壓力的分布特性；Nitka等[18]從細(xì)觀層面上揭示了擋墻土體剪切帶的破壞機(jī)制。以上數(shù)值研究表明采用離散單元法研究土石等離散顆粒材料的側(cè)向壓力分布特性具有獨特優(yōu)勢。

本文采用離散單元法DEM對土、石顆粒材料堆積過程進(jìn)行二維數(shù)值模擬，研究了顆粒粒徑(土和石)、顆粒摩擦系數(shù)對顆粒材料的側(cè)向壓力分布特性的影響規(guī)律，分析了DEM側(cè)向壓力分布特性與經(jīng)典土壓力理論解的差異，基于分形理論探明了不同顆粒粒徑和顆粒摩擦系數(shù)下側(cè)向壓力強(qiáng)度與作用位置分形分布特征，為進(jìn)一步完善擋墻、基坑等工程設(shè)計理論提供參考依據(jù)。

1 側(cè)向壓力分布DEM數(shù)值模擬分析

1.1 DEM數(shù)值計算模型

采用基于DEM理論的離散元軟件PFC2D模擬顆粒材料自然沉積過程：首先分別生成左、右、下三道剛性墻體，在擋墻內(nèi)隨機(jī)生成10 000個粒徑為d的圓形顆粒，并在自重作用下堆積至平衡狀態(tài)。顆粒最終的堆積高度為H0，堆積寬度為L0，如圖1所示。考慮到計算效率，因此采用能夠較合理反映土石力學(xué)性質(zhì)的單一粒徑分布的圓形顆粒。

圖1 DEM數(shù)值計算模型Fig.1 DEM numerical simulation model

參照DEM模擬土、石的相關(guān)文獻(xiàn)[15,16]，顆粒密度取為2 600 kg/m3，顆粒法向和切向剛度取為1.0×108N/m，墻體剛度與顆粒剛度保持一致，墻體摩擦系數(shù)選取為0以便于數(shù)值結(jié)果與朗肯土壓力理論解作對比；本文分別選取石塊、土體顆粒粒徑d為0.40 m、0.01 m，在保證其他數(shù)值計算參數(shù)不變的條件下，通過改變顆粒粒徑d和顆粒摩擦系數(shù)μ以研究兩者對顆粒材料的側(cè)向壓力分布特性的影響規(guī)律。本文最終數(shù)值模擬參數(shù)選取如表1所示。

表1 DEM數(shù)值計算參數(shù)[15,16]Tab.1 DEM numerical simulation parameters

1.2 DEM側(cè)向壓力分布

為研究顆粒粒徑d、顆粒摩擦系數(shù)μ對顆粒材料的側(cè)向壓力分布特性的影響規(guī)律，分別模擬d為0.40 m，0.01 m和μ為0.30、0.60、0.90條件下顆粒堆積過程，其他計算參數(shù)均按表1所取。當(dāng)不同d、μ條件下顆粒堆積穩(wěn)定后，采用PFC2D中內(nèi)嵌FISH語言統(tǒng)計顆粒與左側(cè)墻體接觸力鏈的強(qiáng)度及其位置，得到DEM側(cè)向壓力分布如圖2所示。可以看出，不同d、μ條件下DEM側(cè)向壓力沿深度均呈鋸齒狀分布，這與經(jīng)典朗肯土壓力理論中側(cè)向壓力沿深度呈線性分布的假定不一致。

圖2 不同顆粒粒徑d、顆粒摩擦系數(shù)μ下DEM側(cè)向壓力分布Fig.2 Lateral force distributions obtained from DEM with different particle size d and friction coefficient of particles μ

圖3 DEM力鏈分布中的拱結(jié)構(gòu)效應(yīng)Fig.3 Arch structure effect in the force chains of DEM

顆粒在重力作用下相互接觸形成力鏈，圖3給出了d=0.01 m、μ= 0.30條件下顆粒堆積體力鏈分布圖(右側(cè))和拱結(jié)構(gòu)效應(yīng)圖(左側(cè))，圖中力鏈的粗細(xì)代表了接觸力強(qiáng)度的大小。強(qiáng)度較大的力鏈構(gòu)成了顆粒集合體承載外力的拱結(jié)構(gòu)，對顆粒系統(tǒng)的穩(wěn)定起決定性作用。拱結(jié)構(gòu)效應(yīng)很好地解釋了DEM側(cè)向壓力分布呈鋸齒狀的原因。當(dāng)與墻體有接觸的顆粒為拱結(jié)構(gòu)組成顆粒時(如圖3中顆粒A、D和E)，強(qiáng)度較大的力將通過這些顆粒傳遞給墻體；當(dāng)顆粒沒有位于拱結(jié)構(gòu)時(如顆粒B和C)，強(qiáng)度較小甚至為0(如顆粒C)的力將傳遞給墻體，因此顆粒材料側(cè)向壓力在拱結(jié)構(gòu)效應(yīng)的影響下呈現(xiàn)鋸齒狀分布。

2 不同顆粒粒徑和摩擦系數(shù)下側(cè)向壓力理論解和DEM數(shù)值解

2.1 朗肯土壓力理論

經(jīng)典朗肯土壓力理論認(rèn)為靜止?fàn)顟B(tài)下墻后填土側(cè)向壓力分布沿深度呈線性分布，可表示為：

f=K0ρgH

(1)

式中：f為側(cè)向應(yīng)力分布強(qiáng)度；K0為靜止土壓力系數(shù)，H為堆積深度。

靜止土壓力系數(shù)K0可依據(jù)文獻(xiàn)[19]定義如下:

K0=1-sinφ

(2)

式中：φ為顆粒集合體內(nèi)摩擦角。

根據(jù)文獻(xiàn)[20]可知，顆粒集合體內(nèi)摩擦角φ約等于其休止角度。本文通過DEM分別模擬了μ為0.30、0.60、0.90條件下顆粒集合體在無側(cè)限約束下的堆積過程，得到了相應(yīng)的休止角度(見表2)。由于DEM采用圓形顆粒，沒有考慮顆粒形狀咬合作用，使得表2中數(shù)值計算的內(nèi)摩擦角略小于真實物理試驗測得的。靜止土壓力合力Fr可表達(dá)為：

(3)

根據(jù)朗肯土壓力理論可知合力作用點與左側(cè)墻底距離為Lr=H0/3。依據(jù)上式可得不同d、μ條件下朗肯土壓力理論解，具體如表2所示。

表2 不同顆粒粒徑d、顆粒摩擦系數(shù)μ下DEM和朗肯土壓力理論側(cè)向壓力對比Tab.2 Comparison of lateral force distributions obtained from DEM and Rankine theory with different particle size d and friction coefficient of particles μ

2.2 理論解和DEM數(shù)值解對比

由表2可知，DEM側(cè)向壓力合力值小于朗肯土壓力理論解，當(dāng)顆粒粒徑d較大時，兩者的差異顯著減?。籇EM側(cè)向壓力合力作用點位置與理論解差異很小。隨著顆粒摩擦系數(shù)μ的增大，顆粒之間由摩擦提供的咬合作用力增強(qiáng)，顆粒集合體內(nèi)部成拱率增加，導(dǎo)致DEM側(cè)向壓力合力增大；由圖2可知，DEM某些部分側(cè)向壓力值較周圍值驟增(圖中紅色圓圈所注)，這是朗肯土壓力理論無法計算得到的，工程設(shè)計時應(yīng)考慮這些危險點的存在，例如采用DEM數(shù)值模擬相同工況以得到驟增側(cè)向壓力值及作用位置，在擋墻相應(yīng)位置處進(jìn)行錨固、灌漿等加固處理。

3 DEM側(cè)向壓力分布分形特性研究

分形理論目前被廣泛地應(yīng)用于巖土工程領(lǐng)域中以定量描述表面復(fù)雜無序而又內(nèi)在具有規(guī)律的系統(tǒng)。傳統(tǒng)的歐幾里得幾何學(xué)定義任意形狀直線的維度為一維，平面的維度為二維。分形理論則采用非整數(shù)的形式來描述任意形狀直線(折線)的分布狀態(tài)。本文采用分形理論量化DEM側(cè)向壓力的分布特性，研究不同d、μ條件下DEM側(cè)向壓力強(qiáng)度、作用位置分布內(nèi)在規(guī)律。

3.1 DEM側(cè)向壓力強(qiáng)度分布分形特性研究

根據(jù)分形理論的定義可知[21,22]，當(dāng)采用尺度為r的尺子量取總長度為L的任意形狀直線(折線)時有如下關(guān)系式：

L=Nr

(4)

式中：N為該直線長度所包含尺度r尺子的數(shù)目，如圖4所示。

同時，r、L和該直線的分形維數(shù)D之間具有如下關(guān)系式[20,21]：

L=Cr(1-D)

(5)

式中：C為常數(shù)。

由式(4)及式(5)可得：

N=Cr-D

(6)

通過擬合logr與logN之間的線性關(guān)系得到斜率即分形維數(shù)D。

參照圖4中量取任意形狀直線(折線)的方法，分別統(tǒng)計出不同d、μ條件下DEM側(cè)向壓力強(qiáng)度分布(圖2中的折線)所包含不同尺度r的數(shù)目N(見表3)，并對其分形維數(shù)D進(jìn)行擬合，如圖5所示。分形維數(shù)D越大，表明DEM側(cè)向壓力強(qiáng)度分布的波動性越大。由表3和圖5分析結(jié)果可知，當(dāng)d、μ發(fā)生改變時，DEM側(cè)向壓力分布強(qiáng)度分形維數(shù)D變化較小，分布在1.17～1.29范圍。這表明DEM側(cè)向壓力分布的波動性受d、μ的影響很小，可能是顆粒材料的本質(zhì)特性。

3.2 DEM側(cè)向壓力作用位置分布分形特性研究

顆粒材料的離散特性使得DEM側(cè)向壓力作用位置分布具有非連續(xù)性。采用分形理論研究不同d、μ條件下DEM側(cè)向壓力作用位置分布特性。圖6給出了分析方法：圖中黑色直線代表長度為20 m剛性墻體，藍(lán)色短線代表剛性墻體上側(cè)向壓力作用位置；間距為r的紅色短線將剛性墻體平均分成數(shù)目為20/r(20/r若為分?jǐn)?shù)，則四舍五入取整)個子組，N為子組中包含作用點(藍(lán)色短線)的子組數(shù)目。DEM側(cè)向壓力作用位置分布分形維數(shù)D值在0～1之間，0表示作用位置分布幾乎集中在某一點處，1則表示作用位置分布為線性連續(xù)分布。D越大表明側(cè)向壓力作用位置分布越分散均勻。

圖4 不同尺度r下任意形狀直線(折線)分形分析方法示意圖Fig.4 Schematic diagram of fractal analysis method of arbitrary shaped lines with different ruler r

表3 不同顆粒粒徑d、顆粒摩擦系數(shù)μ下DEM側(cè)向壓力強(qiáng)度分布分形分析
Tab.3 Fractal analysis of lateral force intensities obtained from DEMwith different particle sizedand friction coefficient of particlesμ

顆粒粒徑、摩擦系數(shù)N (r=0.1)N (r=0.2)N (r=0.3)N (r=0.4)分形維數(shù) Dd= 0.01 m、 μ= 0.30139.1160.1230.0526.781.23d= 0.01 m、 μ= 0.60192.2692.9549.0232.531.29d= 0.01 m、 μ= 0.90220.5199.4060.1939.551.23顆粒粒徑、摩擦系數(shù)N (r=2.0)N (r=8.0)N (r=16.0)N (r=32.0 )分形維數(shù) Dd= 0.40 m、 μ= 0.30510.08107.6442.1019.561.19d= 0.40 m、 μ= 0.60362.9767.6024.5010.961.28d= 0.40 m、 μ= 0.90423.4985.3736.9620.131.17

圖5 不同顆粒粒徑d、顆粒摩擦系數(shù)μ下DEM側(cè)向壓力強(qiáng)度分布log r- log N線性擬合Fig.5 Linear fitting relationship between log r and log N for lateral force intensities obtained from DEM with different particle size d and friction coefficient of particles μ

圖6 不同尺度r下側(cè)向力作用位置分布分形分析方法示意圖Fig.6 Schematic diagram of fractal analysis method of lateral force positions with different ruler r

參照上述方法，表4給出了不同d、μ條件下DEM側(cè)向壓力作用位置分布分形特性分析結(jié)果，通過擬合表4中l(wèi)ogr與logN的線性關(guān)系(見圖7)，根據(jù)斜率即可求出D。由表4和圖7可知，DEM側(cè)向壓力作用位置分布具有良好的分形特性；當(dāng)顆粒粒徑d較大時，D較小，表明側(cè)向壓力作用位置分布較集中；當(dāng)顆粒摩擦系數(shù)μ增大時，D減小，這表明側(cè)向壓力作用位置分布整體傾向于集中密集，剛性擋墻的某些部位將承載更多的側(cè)向壓力，這將易導(dǎo)致?lián)鯄植慨a(chǎn)生破壞，不利于其安全穩(wěn)定。

表4 不同顆粒粒徑d、顆粒摩擦系數(shù)μ下DEM側(cè)向壓力作用位置分布分形分析Tab.4 Fractal analysis of lateral force positions obtained from DEM with different particle size d and friction coefficient of particles μ

圖7 不同顆粒粒徑d和顆粒摩擦系數(shù)μ下DEM側(cè)向壓力作用位置分布log r- log N線性擬合Fig.7 Linear fitting relationship between log r and log N for lateral force positions obtained from DEM with different particle size d and friction coefficient of particles μ

4 結(jié) 論

(1)經(jīng)典朗肯土壓力理論假定土石等顆粒材料為連續(xù)介質(zhì)，有悖于其離散特性的本質(zhì)，計算得到側(cè)向壓力分布隨堆積深度均勻增大且連續(xù)分布。由于顆粒材料堆積后內(nèi)部呈現(xiàn)的拱結(jié)構(gòu)效應(yīng)，DEM模擬得到的側(cè)向壓力分布呈鋸齒狀且具有不連續(xù)性，部分側(cè)向壓力值較周圍值驟增，在工程設(shè)計時應(yīng)結(jié)合相同工況下DEM數(shù)值模擬分析得到的驟增側(cè)向壓力值及作用位置，在擋墻相應(yīng)位置處進(jìn)行錨固、灌漿等加固處理。

(2)DEM側(cè)向壓力合力值小于朗肯土壓力理論解，當(dāng)顆粒粒徑增大時，顆粒側(cè)向壓力合力與朗肯土壓力理論解差異顯著減小。兩者計算得到的合力作用點位置幾乎沒有差異。

(3)采用分形理論分析了不同顆粒粒徑和顆粒摩擦系數(shù)條件下DEM側(cè)向壓力分布特性。側(cè)向壓力強(qiáng)度分布分形維數(shù)受顆粒粒徑和顆粒摩擦系數(shù)影響不大，分布范圍在1.17～1.29。側(cè)向壓力分布位置受顆粒粒徑和顆粒摩擦系數(shù)影響較大；顆粒粒徑和顆粒摩擦系數(shù)越大，側(cè)向壓力作用位置分形維數(shù)越小，作用力分布位置越密集，越不利于剛性擋墻的安全穩(wěn)定。

□