無黏結(jié)材料顆粒流模型的宏細(xì)觀參數(shù)關(guān)系研究

顏 敬，曾亞武，高 睿，杜 欣

(武漢大學(xué)土木建筑工程學(xué)院，武漢 430072)

1 研究背景

離散單元法(如PFC2D/3D，UDEC/3DEC)是近些年來興起的一種新的巖土工程數(shù)值計算方法。該方法一般按其用途，又可以分為細(xì)觀離散元和宏觀離散元。前者著重于處理數(shù)目眾多、具有不連續(xù)特征的接觸面或點的問題，例如破碎巖體中的破裂面、砂土中的接觸面(點)、散體材料中粒子之間的接觸面(點)等;后者則主要解決規(guī)模相對較大的不連續(xù)面，如斷層、節(jié)理、結(jié)構(gòu)與基礎(chǔ)之間的結(jié)合面等引起的一系列問題等［1－4］。PFC(Particle Flow Code)是在巖土界著名學(xué)者Peter Cundall院士主持下采用細(xì)觀離散元理論(又稱為粒子流或者顆粒流理論)開發(fā)的一套商業(yè)數(shù)值計算軟件，可以廣泛地應(yīng)用于研究細(xì)觀結(jié)構(gòu)控制的問題，具有很大的發(fā)展?jié)摿Α?/p>

ITASCA公司開發(fā)的PFC系列軟件，作為細(xì)觀離散元軟件平臺，具有2個最基本的特征［2－4］。

(1)允許顆粒發(fā)生有限的位移和轉(zhuǎn)動，顆粒之間允許完全脫離;

(2)程序在計算過程中能夠自動辨識新的接觸。該理論的核心思想是采用最基本的單元——粒子和最基本的力學(xué)關(guān)系——牛頓第二定律來描述介質(zhì)的復(fù)雜力學(xué)行為，而無需事先確定材料的宏觀本構(gòu)，模型介質(zhì)的本構(gòu)特征將由其內(nèi)部顆粒之間狀態(tài)特征的變化體現(xiàn)出來，顆粒間接觸的破壞與發(fā)展意味著介質(zhì)整體力學(xué)特性由線性向非線性轉(zhuǎn)化。

PFC的這種理念導(dǎo)致其在很大程度上不同于其他連續(xù)和非連續(xù)介質(zhì)力學(xué)理論方法與程序，其中的一個關(guān)鍵問題就是:介質(zhì)的宏觀基本物理力學(xué)特征不能通過直接賦值的形式實現(xiàn)，只有顆粒的幾何特性和顆粒間接觸的細(xì)觀力學(xué)參數(shù)可以賦值，并且這些參數(shù)大都無法測量，但介質(zhì)的宏觀力學(xué)特征又取決于顆粒的這些基本特性，一旦改變了這些基本特性，就意味著改變了介質(zhì)的宏觀力學(xué)性質(zhì)。

正是由于PFC具有這些特點，在投入工程計算之前，必須利用同樣的粒子幾何參數(shù)建立簡單的室內(nèi)試樣模型，并對這些試樣分別賦予不同的細(xì)觀力學(xué)參數(shù)進(jìn)行一系列的數(shù)值試驗，以獲得試樣的宏觀力學(xué)參數(shù)，并將這些測試結(jié)果與實際工程相應(yīng)的宏觀參數(shù)值進(jìn)行比較，選擇對應(yīng)的細(xì)觀力學(xué)參數(shù)作為實際模型計算用參數(shù)，即要進(jìn)行一個復(fù)雜的細(xì)觀力學(xué)參數(shù)標(biāo)定工作。這個過程與地質(zhì)力學(xué)材料模型試驗前期工作中選擇和配制模型材料的配合比非常類似［2］。但由于宏觀層次和細(xì)觀層次的參數(shù)眾多，盲目調(diào)試極有可能事倍功半，且極難尋找一個明確的數(shù)學(xué)模型來表征二者之間的關(guān)系，目前采用最多的方法就是變量控制法，即保持某些細(xì)觀參量不變，改變其他參數(shù)取值，探究介質(zhì)宏觀特性的演化規(guī)律，例如國內(nèi)外學(xué)者 S.Utili和R.Nova［5］，C.Wang，D.D.Tannant和P.A.Lilly［6］，Geng Yan，Yu Hai-sui和McDowell Glenn［7］，同 濟 大學(xué) 周 ?。?］、Chengbing WANG 和Hehua ZHU［9］，長江科學(xué)院李耀旭、王新和丁秀麗［10－11］等都開展了相關(guān)的研究工作。

本文擬在前人研究的基礎(chǔ)之上，以某無黏結(jié)顆粒材料為例，利用正交試驗方法設(shè)計試樣并在不同側(cè)壓下進(jìn)行雙軸試驗，以探求細(xì)觀參數(shù)不同組合對介質(zhì)宏觀特性的影響，從而避免控制變量法固定某些參數(shù)的局限，更加科學(xué)地分析細(xì)觀參量對介質(zhì)宏觀特性影響的敏感程度，并進(jìn)一步地利用人工神經(jīng)網(wǎng)絡(luò)實現(xiàn)兩個層次參數(shù)間的互演，以供PFC模型進(jìn)行實際巖土工程計算時參考。

2 正交試驗設(shè)計

在實際生產(chǎn)和科學(xué)研究中，往往需要通過一定的試驗來獲取一些試驗數(shù)據(jù)，并對這些數(shù)據(jù)進(jìn)行科學(xué)分析和數(shù)學(xué)處理，以幫助人們找出問題的主要矛盾及其間的相互關(guān)系，明確問題的內(nèi)在規(guī)律，從而尋求問題的解決辦法。對于多因素影響的試驗，如果采用完全組合將會導(dǎo)致試驗次數(shù)急劇增加且試驗結(jié)果整理困難，而正交試驗可以用最少的試驗次數(shù)反映比較全面的情況，并帶有充足的變異信息，可以反映不同因素水平組合對結(jié)果的影響［12］。

本文以無黏結(jié)材料中顆粒的法向剛度kn、剛度比kn/ks(法向與切向剛度之比)、摩擦系數(shù)f、平均半徑Rm這4個細(xì)觀參數(shù)進(jìn)行3水平正交設(shè)計(依正交表L934，文獻(xiàn)［12］)，形成9 類試樣，具體見表1。

表1 試樣設(shè)計Table 1 Orthogonal designs for the samples

3 PFC雙軸試驗

PFC2D軟件用4個墻體限制試樣形狀為矩形，頂墻和底墻模擬加荷板，左墻和右墻模擬試樣側(cè)限。在數(shù)值試驗過程中，通過指定頂板和底板的速度施加軸壓，采用平面應(yīng)變控制式加載，而左右墻體的速度均由數(shù)值伺服器自動控制以保持側(cè)壓恒定(圖1)，具體試驗原理見文獻(xiàn)［3］。本文選取試樣大小為300mm×600mm，孔隙率取為0.1。

試樣的應(yīng)力狀態(tài)通過邊墻的平均接觸應(yīng)力表示(見圖2)，即

式中:nc表示與墻體接觸的顆粒個數(shù);Fi表示接觸集中力;ld是墻體面積，l為墻體有效長度，d表示沿縱向的厚度(一般取單位厚度)。

圖1 試驗裝置Fig.1 Testing device

圖2 顆粒-墻體接觸Fig.2 Particle-wall contact

依據(jù)PFC使用手冊，試樣的應(yīng)變狀態(tài)依照下式進(jìn)行計算(大應(yīng)變模型)，即

式中:L為變形后的長度;L0為變形前的長度。

體積應(yīng)變?yōu)閤，y方向應(yīng)變之和，即

若側(cè)壓為σc，軸向主應(yīng)力和偏差應(yīng)力分別為σa和σd，則有彈性模量和泊松比為:

本文側(cè)壓分別取10，15，20 kPa，記錄每一個試樣在3種側(cè)壓下的偏應(yīng)力－軸應(yīng)變曲線、體應(yīng)變－軸應(yīng)變曲線、破壞偏差應(yīng)力，測試結(jié)果見表2。

4 數(shù)值試驗結(jié)果整理與分析

4.1 變形參數(shù)E，v

通過27次雙軸試驗，發(fā)現(xiàn)偏應(yīng)力－軸應(yīng)變曲線均表現(xiàn)為圖3的特征，在原點附近加載曲線可以近似為一條直線，介質(zhì)呈現(xiàn)線彈性響應(yīng)，之后進(jìn)入非線性階段，達(dá)到峰值強度后破壞。依照每一次試驗的記錄，計算試樣在線性階段的彈性模量和泊松比(初始彈性模量和泊松比，即變形參數(shù)，見表3)，并進(jìn)行統(tǒng)計分析(見表4)。

表2 破壞偏差應(yīng)力Table 2 Maximum deviatoric stresses kPa

圖3 應(yīng)力-應(yīng)變曲線Fig.3 Stress-strain curve

表3 彈性模量與泊松比Table 3 Elastic moduli and poisson’s ratios

表4 彈性模量和泊松比的統(tǒng)計分析Table 4 Statistic analysis of elastic moduli and poisson’s ratios

觀察試驗結(jié)果可以發(fā)現(xiàn):在一定范圍內(nèi)，側(cè)壓應(yīng)力不同，試樣的彈性模量和泊松比略顯差別，但是起伏不大，極差與平均值的比值均在10%以內(nèi)，即可以用均值表示它們的大小，而無需嚴(yán)格考慮側(cè)壓應(yīng)力的變化。

按照正交試驗的結(jié)果分析方法［12］，可以分析各個細(xì)觀參數(shù)取值水平對該種材料彈模和泊松比的影響，形成圖4和圖5，并以彈模為例給出數(shù)據(jù)的處理過程，見表5和表6，泊松比和下文摩擦角的處理與此類似。

圖4 細(xì)觀參數(shù)取值對彈性模量的影響Fig.4 The influence of particle mesoscopic parameter on elastic modulus

圖5 細(xì)觀參數(shù)取值對泊松比的影響Fig.5 The influence of particle mesoscopic parameter on poisson’s ratio

表5 各因素水平試樣的彈性模量Table 5 Analysis of elastic moduli of the level of each factor

分析圖4、圖5、表5和表6可以得出結(jié)論:

(1)法向剛度與彈性模量呈現(xiàn)顯著的正相關(guān)關(guān)系，剛度比增大即切向剛度減小可以降低彈性模量。顆粒的摩擦系數(shù)和平均粒徑對彈模的影響較小。

(2)極差R的大小，反映了試驗中各因素作用的大小，極差大表明該因素對指標(biāo)的影響大。對于彈性模量而言，各因素作用由強到弱排序為:法向剛度、剛度比、平均粒徑、摩擦系數(shù)。

(3)顆粒剛度在法向和切向的大小差異(剛度比)是影響泊松比的最主要因素，剛度比越大，介質(zhì)體在宏觀上的泊松比就越大。

(4)顆粒法向剛度的增加和粒徑的減少可以導(dǎo)致泊松比減小，顆粒摩擦系數(shù)對泊松比幾乎沒有影響。

(5)對于泊松比而言，各因素作用由強到弱排序為:剛度比、法向剛度、平均粒徑、摩擦系數(shù)。

表6 彈性模量統(tǒng)計分析Table 6 Statistical analysis of elastic moduli kPa

4.2 強度參數(shù)φ

對于無黏結(jié)的顆粒材料，根據(jù)摩爾庫侖強度理論，達(dá)到破壞極限時，大小主應(yīng)力應(yīng)滿足其中內(nèi)摩擦角度σ3－σ1平面上破壞點應(yīng)該在一條過原點的直線上，將試驗結(jié)果繪成曲線，見圖6，可以觀察到試驗結(jié)果基本滿足這一規(guī)律，故可認(rèn)為在一定應(yīng)力水平下摩爾庫侖強度理論成立且用方程σ1=λσ3進(jìn)行最小二乘回歸可以得到內(nèi)摩擦角，通過進(jìn)一步分析可以得出細(xì)觀參數(shù)對內(nèi)摩擦角的影響(見圖7)。

圖6 試樣破壞時的σ1-σ3關(guān)系Fig.6 Relation between σ1 and σ3 when samples fail

(1)內(nèi)摩擦角與顆粒摩擦系數(shù)呈現(xiàn)顯著的正相關(guān)關(guān)系，法向剛度的增加可以適當(dāng)增加內(nèi)摩擦角。

(2)對于內(nèi)摩擦角而言，各因素作用由強到弱排序為:摩擦系數(shù)、法向剛度、剛度比和平均粒徑。

圖7 細(xì)觀參數(shù)取值對內(nèi)摩擦角的影響Fig.7 The influence of particle mesoscopic parameter on internal friction angle

4.3 宏細(xì)觀參數(shù)匹配調(diào)整原則

在進(jìn)行實際巖土工程計算時，要求實際物理試驗測試出的介質(zhì)宏觀參數(shù)與PFC數(shù)值試驗測試出的介質(zhì)宏觀特性參數(shù)盡可能的一致，這就需要合理調(diào)整顆粒流模型的細(xì)觀參數(shù)，如果盲目調(diào)試則極有可能事倍功半。根據(jù)上面的分析，針對本文討論的無黏結(jié)顆粒材料，可以得到一些基本的調(diào)整原則。

(1)僅增大(減小)顆粒摩擦系數(shù)可以使內(nèi)摩擦角增大(減小)，而保持介質(zhì)彈性模量和泊松比基本不變;

(2)增加法向剛度的同時，增加剛度比，可以提高介質(zhì)體的彈性模量，且保證內(nèi)摩擦角基本不變;

(3)保持法向剛度不變，增加剛度比，可以增大介質(zhì)體的泊松比，且保持彈性模量和內(nèi)摩擦角度基本不變;

(4)在一定粒徑范圍內(nèi)，改變顆粒的大小對宏觀參數(shù)的影響十分有限;

(5)在處理問題時，應(yīng)該抓住主要矛盾，將逼近主要宏觀參數(shù)作為調(diào)整目標(biāo)，而適當(dāng)放寬對其他參數(shù)的模擬要求。

5 宏細(xì)觀參數(shù)間的智能互演

利用人工神經(jīng)網(wǎng)絡(luò)［13］的非線性、動態(tài)性與魯棒性的特點，以正交試驗結(jié)果為訓(xùn)練樣本(9個)，建立宏觀細(xì)觀參數(shù)之間的相互映射，即互演機制，以實現(xiàn)2個層次參數(shù)間的相互預(yù)測。在Matlab的平臺上，自編程序構(gòu)建3層BP網(wǎng)絡(luò)系統(tǒng)。

5.1 正演計算

建立列向量之間的映射關(guān)系:

訓(xùn)練形成網(wǎng)絡(luò)net1，并測試得出(試樣粒度分布為3～7mm)

將細(xì)觀參數(shù)輸入PFC雙軸試驗程序以檢驗預(yù)測結(jié)果的正確性，得到表7。將破壞時的(σ3，σ1)標(biāo)在主應(yīng)力坐標(biāo)平面，用方程σ1=λσ3進(jìn)行回歸，得到 λ=2.53，進(jìn)一步計算 φ=25.75°，可以發(fā)現(xiàn)預(yù)測結(jié)果和數(shù)值試驗結(jié)果比較接近(見表8)。

表7 正演計算測試結(jié)果Table 7 Test results of forward computation

表8 正演計算結(jié)果校驗Table 8 Verification results of forward computation

5.2 反演計算

建立列向量之間的映射關(guān)系:

訓(xùn)練形成網(wǎng)絡(luò) net2，并測試得出(試樣粒度分布3.22～7.22mm)

將細(xì)觀參數(shù)輸入PFC雙軸試驗程序以檢驗預(yù)測結(jié)果的正確性，得到表8。同理可得λ=2.54，φ=25.92°，可見預(yù)測結(jié)果和數(shù)值試驗結(jié)果比較接近，見表9和表10。

表9 反演計算測試結(jié)果Table 9 Test results of inversion

表10 反演計算結(jié)果校驗Table 10 Verification results of inversion

6 結(jié)論

本文以無黏結(jié)顆粒材料為例，通過一組正交設(shè)計的試驗，討論了PFC顆粒細(xì)觀參數(shù)取值對介質(zhì)宏觀特性的影響，并分析了各個細(xì)觀參量對宏觀參量的敏感度，提出了數(shù)值實驗時宏細(xì)觀參數(shù)匹配調(diào)整的基本原則，并利用神經(jīng)網(wǎng)絡(luò)實現(xiàn)了宏細(xì)觀參數(shù)間的智能互演，為PFC系列軟件在實際巖土工程計算中提供了參考。黏結(jié)顆粒材料例如黏土、巖石、混凝土等亦可按照本文的研究思路進(jìn)行分析。但是本文在分析試樣破壞大小應(yīng)力關(guān)系時，采用了摩爾庫侖強度理論，認(rèn)為二者為線性關(guān)系，實際上隨著側(cè)壓的提高，主應(yīng)力并未嚴(yán)格地按線性比例增大，而是呈現(xiàn)非線性變化，σ3－σ1關(guān)系曲線出現(xiàn)逐漸平緩的趨勢(這與莫爾理論極為相似)，故可以考慮引入其他的宏觀本構(gòu)模型并研究細(xì)觀參量對本構(gòu)參數(shù)的影響。

［1］CUNDALL P A，STACK O D L.A Discrete Numerical Model for Granular Assemblies［J］.Geotechique，1979，29(1):47－65.

［2］朱煥春.PFC及其在礦山崩落開采研究中的應(yīng)用［J］.巖石力學(xué)與工程學(xué)報，2006，25(9):1927－1931.(ZHU Huan-chun.PFC and Application Case of Caving Study［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25(9):1927－1931.(in Chinese))

［3］Itasca Consulting Group，Inc.，PFC2D(Particle Flow Code in 2 Dimensions)，Version1.1［K］.Minneapolis，MN:ICG，1999.

［4］Itasca Consulting Group，PFC3D(Particle Flow Code in 3 Dimensions)，Version1.1［K］.Minneapolis，MN:ICG，1999.

［5］UTILI S，NOVA R.DEM Analysis of Bonded Granular Geomaterials［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2008，32(17):1997－2031.

［6］WANG C，TANANT D D，LILLY PA.Numerical Analysis of the Stability of Heavily Jointed Rock Slopes Using PFC2D［J］.International Journal of Rock Mechanics＆Mining Sciences，2003，40(3):415－424.

［7］GENG Y，YU H S，GLENN M.Simulation of Granular Material Behaviour Using DEM［J］.Procedia Earth and Planetary Science，2009，1(1):598－605.

［8］周 健，王家全，曾 遠(yuǎn)，等.顆粒流強度折減法和重力增加法的邊坡安全系數(shù)研究［J］.巖土力學(xué)，2009，30(6):1549－1554.(ZHOU Jian，WANG Jia-quan，ZENG Yuan，et al.Slope Safety Factor by Methods of Particle Flow Code Strength Reduction and Gravity Increase［J］.Chinese Journal of Rock and Soil Mechanics，2009，30(6):1549－1554.(in Chinese))

［9］WANG C B，ZHU H H.Study on the Relationship Between Particle Micro-parameters and Soil Mechanical Properties［C］∥ Proceedings of the 1st International Conference on Information Technology in Geo-Engineering(ICITG).Shanghai，China，Sep.16－17，2010:599－606.

［10］李耀旭.顆粒流方法在土石混合體力學(xué)特性研究中的應(yīng)用［D］.武漢:長江科學(xué)院，2009.(LI Yao-xu.The Application of Particle Flow Code in Studying of Mechanical Characteristics of Soil-Rock Mixture［D］.Wuhan:Yangtze River Scientific Research Institute，2009.(in Chinese))

［11］王 新.土石混合體力學(xué)特性影響因素及破環(huán)機制研究［D］.武漢:長江科學(xué)院，2010.(WANG Xin.Research on Influence Factors of Mechanics Characteristics and Failure Mechanism of Soil-Rock Mixture［D］.Wuhan:Yangtze River Scientific Research Institute，2010.(in Chinese))

［12］吳貴生.試驗設(shè)計與數(shù)據(jù)處理［M］.北京:冶金工業(yè)出版社，1997.(WU Gui-sheng.Test Design and Data Processing［M］.Beijing:Metallurgy Industry Press，1997.(in Chinese))

［13］張德豐.Matlab神經(jīng)網(wǎng)絡(luò)仿真與應(yīng)用［M］.北京:電子工業(yè)出版社，2009.(ZHANG De-feng.Simulation and Application of Neural Network Based on Matlab［M］.Beijing:Publishing House of Electronics Industry，2009.(in Chinese))