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

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

1 研究背景

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

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

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

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

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

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

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

2 正交試驗設計

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

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

表1 試樣設計Table 1 Orthogonal designs for the samples

3 PFC雙軸試驗

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

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

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

圖1 試驗裝置Fig.1 Testing device

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

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

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

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

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

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

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

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

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

表2 破壞偏差應力Table 2 Maximum deviatoric stresses kPa

圖3 應力-應變曲線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

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

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

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

圖5 細觀參數(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可以得出結論:

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

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

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

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

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

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

4.2 強度參數(shù)φ

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5.1 正演計算

建立列向量之間的映射關系:

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

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

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

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

5.2 反演計算

建立列向量之間的映射關系:

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

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

表9 反演計算測試結果Table 9 Test results of inversion

表10 反演計算結果校驗Table 10 Verification results of inversion

6 結論

本文以無黏結顆粒材料為例，通過一組正交設計的試驗，討論了PFC顆粒細觀參數(shù)取值對介質(zhì)宏觀特性的影響，并分析了各個細觀參量對宏觀參量的敏感度，提出了數(shù)值實驗時宏細觀參數(shù)匹配調(diào)整的基本原則，并利用神經(jīng)網(wǎng)絡實現(xiàn)了宏細觀參數(shù)間的智能互演，為PFC系列軟件在實際巖土工程計算中提供了參考。黏結顆粒材料例如黏土、巖石、混凝土等亦可按照本文的研究思路進行分析。但是本文在分析試樣破壞大小應力關系時，采用了摩爾庫侖強度理論，認為二者為線性關系，實際上隨著側壓的提高，主應力并未嚴格地按線性比例增大，而是呈現(xiàn)非線性變化，σ3－σ1關系曲線出現(xiàn)逐漸平緩的趨勢(這與莫爾理論極為相似)，故可以考慮引入其他的宏觀本構模型并研究細觀參量對本構參數(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及其在礦山崩落開采研究中的應用［J］.巖石力學與工程學報，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］周 健，王家全，曾 遠，等.顆粒流強度折減法和重力增加法的邊坡安全系數(shù)研究［J］.巖土力學，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］李耀旭.顆粒流方法在土石混合體力學特性研究中的應用［D］.武漢:長江科學院，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］王 新.土石混合體力學特性影響因素及破環(huán)機制研究［D］.武漢:長江科學院，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ù)據(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)絡仿真與應用［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))