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

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

0 引 言

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

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

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

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

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

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

圖1 DEM數(shù)值計(jì)算模型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é)果與朗肯土壓力理論解作對(duì)比；本文分別選取石塊、土體顆粒粒徑d為0.40 m、0.01 m，在保證其他數(shù)值計(jì)算參數(shù)不變的條件下，通過改變顆粒粒徑d和顆粒摩擦系數(shù)μ以研究?jī)烧邔?duì)顆粒材料的側(cè)向壓力分布特性的影響規(guī)律。本文最終數(shù)值模擬參數(shù)選取如表1所示。

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

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

為研究顆粒粒徑d、顆粒摩擦系數(shù)μ對(duì)顆粒材料的側(cè)向壓力分布特性的影響規(guī)律，分別模擬d為0.40 m，0.01 m和μ為0.30、0.60、0.90條件下顆粒堆積過程，其他計(jì)算參數(shù)均按表1所取。當(dāng)不同d、μ條件下顆粒堆積穩(wěn)定后，采用PFC2D中內(nèi)嵌FISH語言統(tǒng)計(jì)顆粒與左側(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)，對(duì)顆粒系統(tǒng)的穩(wěn)定起決定性作用。拱結(jié)構(gòu)效應(yīng)很好地解釋了DEM側(cè)向壓力分布呈鋸齒狀的原因。當(dāng)與墻體有接觸的顆粒為拱結(jié)構(gòu)組成顆粒時(shí)(如圖3中顆粒A、D和E)，強(qiáng)度較大的力將通過這些顆粒傳遞給墻體；當(dāng)顆粒沒有位于拱結(jié)構(gòu)時(shí)(如顆粒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ù)值計(jì)算的內(nèi)摩擦角略小于真實(shí)物理試驗(yàn)測(cè)得的。靜止土壓力合力Fr可表達(dá)為：

(3)

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

表2 不同顆粒粒徑d、顆粒摩擦系數(shù)μ下DEM和朗肯土壓力理論側(cè)向壓力對(duì)比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ù)值解對(duì)比

由表2可知，DEM側(cè)向壓力合力值小于朗肯土壓力理論解，當(dāng)顆粒粒徑d較大時(shí)，兩者的差異顯著減小；DEM側(cè)向壓力合力作用點(diǎn)位置與理論解差異很小。隨著顆粒摩擦系數(shù)μ的增大，顆粒之間由摩擦提供的咬合作用力增強(qiáng)，顆粒集合體內(nèi)部成拱率增加，導(dǎo)致DEM側(cè)向壓力合力增大；由圖2可知，DEM某些部分側(cè)向壓力值較周圍值驟增(圖中紅色圓圈所注)，這是朗肯土壓力理論無法計(jì)算得到的，工程設(shè)計(jì)時(shí)應(yīng)考慮這些危險(xiǎn)點(diǎn)的存在，例如采用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的尺子量取總長(zhǎng)度為L(zhǎng)的任意形狀直線(折線)時(shí)有如下關(guān)系式：

L=Nr

(4)

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

同時(shí)，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。

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

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

顆粒材料的離散特性使得DEM側(cè)向壓力作用位置分布具有非連續(xù)性。采用分形理論研究不同d、μ條件下DEM側(cè)向壓力作用位置分布特性。圖6給出了分析方法：圖中黑色直線代表長(zhǎng)度為20 m剛性墻體，藍(lán)色短線代表剛性墻體上側(cè)向壓力作用位置；間距為r的紅色短線將剛性墻體平均分成數(shù)目為20/r(20/r若為分?jǐn)?shù)，則四舍五入取整)個(gè)子組，N為子組中包含作用點(diǎn)(藍(lán)色短線)的子組數(shù)目。DEM側(cè)向壓力作用位置分布分形維數(shù)D值在0～1之間，0表示作用位置分布幾乎集中在某一點(diǎn)處，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較大時(shí)，D較小，表明側(cè)向壓力作用位置分布較集中；當(dāng)顆粒摩擦系數(shù)μ增大時(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ì)，計(jì)算得到側(cè)向壓力分布隨堆積深度均勻增大且連續(xù)分布。由于顆粒材料堆積后內(nèi)部呈現(xiàn)的拱結(jié)構(gòu)效應(yīng)，DEM模擬得到的側(cè)向壓力分布呈鋸齒狀且具有不連續(xù)性，部分側(cè)向壓力值較周圍值驟增，在工程設(shè)計(jì)時(shí)應(yīng)結(jié)合相同工況下DEM數(shù)值模擬分析得到的驟增側(cè)向壓力值及作用位置，在擋墻相應(yīng)位置處進(jìn)行錨固、灌漿等加固處理。

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

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

□