邵曉泉，遲世春*，張宗亮

(1.大連理工大學(xué) 海岸與近海工程國家重點實驗室，遼寧 大連 116024；2.大連理工大學(xué) 建設(shè)工程學(xué)部水利工程學(xué)院工程抗震研究所，遼寧 大連 116024；3.中國電建集團 昆明勘測設(shè)計研究院有限公司，云南 昆明 650051)

土石壩因其建造速度快、對基礎(chǔ)適應(yīng)能力強等優(yōu)點已經(jīng)成為應(yīng)用最廣的壩型之一。目前,在建或者擬建的200～300 m級土石壩多位于地震頻發(fā)的區(qū)域，其抗震安全問題至關(guān)重要。堆石料作為土石壩的主要支撐材料其動力特性研究得到了廣泛關(guān)注。

一般來說，堆石料的最大動模量隨著圍壓和固結(jié)應(yīng)力比的增加而增大[1-3]。相同循環(huán)荷載作用下，母巖的巖性越差，顆粒破碎率越大[4]。低圍壓下堆石料的動力特性受級配優(yōu)劣的影響并不顯著，但隨著圍壓增加，級配優(yōu)的試樣動力特性指標(biāo)更優(yōu)。對渾圓度較好的砂礫料，級配優(yōu)劣對顆粒破碎率的影響較小，相應(yīng)地，對動力變形特性的影響逐漸減弱[5]。施加動載前的顆粒破碎情況會影響后續(xù)動荷載作用過程中的顆粒破碎量，且顆粒破碎是引起試樣體積收縮的主要因素[6]。由此可知顆粒破碎是影響堆石料動力特性的重要方面。

離散單元法是從細觀角度分析巖土體靜動力特性的有效手段[7-9]。劉漢龍等[10]采用2維離散元方法對筑壩反濾料低圍壓下的模量阻尼特性進行模擬，提出Hertz-Mindlin非線性接觸模型可以較好地反映材料的非線性特征。楊貴等[11]根據(jù)中小尺寸堆石料動模量變化規(guī)律采用離散元模擬預(yù)測大尺寸試樣的模量值。通過引入破碎準(zhǔn)則研究了循環(huán)荷載作用下圍壓和固結(jié)應(yīng)力比對顆粒破碎的影響[12]。但是目前進行堆石料動力特性模擬時均采用圓盤或圓球模擬堆石顆粒，未考慮顆粒形狀的不規(guī)則性且考慮顆粒破碎對動模量影響的數(shù)值模擬研究較少。

以糯扎渡Ⅱ區(qū)粗堆石料為研究對象，采用多個小球組合模擬不規(guī)則的顆粒形狀，引入顆粒破碎準(zhǔn)則，通過擬合不同圍壓下動三軸骨架曲線確定一套細觀參數(shù)。在此基礎(chǔ)上研究了顆粒破碎和孔隙率對堆石料動彈性模量的影響規(guī)律。

1 數(shù)值試驗

1.1 試樣制備與加載

堆石顆粒采用多個直徑相等的小球組合而成。首先，生成長寬比范圍為1.2～1.8之間的橢球體，在橢球面上隨機選取50個點作為多面體的頂點[13]。隨后，按最密六方排列生成一系列小球，小球完全覆蓋多面體的區(qū)域。最后，通過判斷小球與多面體的位置關(guān)系，刪掉處于多面體外部的小球，剩下的小球組合得到不規(guī)則形狀的顆粒模板。部分顆粒模板如圖1（a）所示。

圖1 數(shù)值模擬的不規(guī)則顆粒與級配曲線Fig.1 Typical particles with irregular shape in simulations and grain size distribution

試驗材料為糯扎渡Ⅱ區(qū)粗堆石料，巖性為花崗巖[14]。為保證數(shù)值模擬的可行性和計算效率，研究者多采用截斷級配的方式[15-16]。張宜等[17]指出堆石體的強度和變形特性隨截斷粒徑的增大而單調(diào)變化，且截斷粒徑取15 mm時其影響可以接受。截斷后的模擬級配與室內(nèi)試驗級配如圖1（b）所示。采用PFC內(nèi)置的Hertz接觸模型來模擬堆石料的非線性特征[10]。制樣過程中通過給定不同的摩擦系數(shù)可以得到不同孔隙率的試樣[18]。模擬時不考慮墻體與顆粒之間的摩擦，制樣過程中不考慮破碎且顆粒間的摩擦系數(shù)設(shè)為0.15。制樣時,首先在給定空間范圍內(nèi)生成滿足級配要求的球形顆粒。然后將球形顆粒替換為不規(guī)則顆粒，替換完成后緩慢移動墻體將試樣壓縮至30 kPa圍壓下固結(jié)穩(wěn)定。最后將摩擦系數(shù)設(shè)為其真實值。生成的數(shù)值試樣尺寸為300 mm×300 mm×300 mm，由3 519個顆粒共84 310個小球組成，初始孔隙率為0.55。該孔隙率大于室內(nèi)試樣的孔隙率主要是因為模擬時對級配進行截斷處理導(dǎo)致沒有足夠的小顆粒填充孔隙[19]。Huang等[20]模擬循環(huán)荷載下砂土的力學(xué)特性時設(shè)置局部阻尼系數(shù)為0.1，并通過敏感性分析發(fā)現(xiàn)振動角頻率ω小于10π rad/s時模擬結(jié)果幾乎不受影響。局部阻尼系數(shù)和振動角頻率均參考其取值。

動三軸試驗圍壓為900、1 500和2 200 kPa，固結(jié)比Kc=σ1/σ3為1.5，加載工況如表1所示。數(shù)值模擬時先將試樣固結(jié)到指定圍壓，穩(wěn)定后施加偏應(yīng)力直到試樣的固結(jié)比為1.5。動力試驗采用應(yīng)力控制，加載波型為正弦波。一個試樣施加一個動應(yīng)力σd，每一級動應(yīng)力條件下振動3圈，取第3圈計算動彈性模量。不排水循環(huán)加載條件的模擬與劉漢龍等學(xué)者的處理方法一致[10,20-21]。循環(huán)加載過程中固定上側(cè)墻體，每循環(huán)一步均對墻體的速度進行調(diào)整。先根據(jù)當(dāng)前時間和時步大小計算下一循環(huán)步結(jié)束對應(yīng)的時間t，計算得到動應(yīng)力的目標(biāo)值為σdsin(ωt)。通過伺服機制指定底部墻體的運動速度使試樣動應(yīng)力滿足要求。再根據(jù)底部墻體的速度計算周圍四面墻體的運動方向與速度，確保整個振動過程中試樣的體積不變。

表1 模擬的試驗方案Tab. 1 Test schemes for dynamic simulation

1.2 破碎準(zhǔn)則

當(dāng)顆粒受到的最大接觸力或等效應(yīng)力達到極限狀態(tài)時顆粒會發(fā)生破碎[22-23]。Russell等[24]提出遭受多個外荷載的顆粒，其破碎主要受最大接觸力控制。在其工作的基礎(chǔ)上，Ciantia等[25]提出一個通用高效的破壞準(zhǔn)則。如圖2所示，當(dāng)最大法向接觸力滿足如下關(guān)系式時，顆粒發(fā)生破碎。

圖2 接觸力作用于球體表面Fig.2 Contact force acts to sphere surface

式中：Fmax、d分別為作用在顆粒上的最大接觸力和顆粒直徑；θ0決定了接觸面積的大小，根據(jù)赫茲接觸理論確定[25]；σlim為極限強度，與材料的泊松比、單軸抗壓和抗拉強度相關(guān)。

極限強度σlim假設(shè)符合Weibull分布，且極限強度是粒徑相關(guān)的：

式中：m為Weibull模量，表征顆粒強度分布的離散程度；d0為參考粒徑，本文選為38 mm；σlim,0為粒徑d0的顆粒殘存概率37%對應(yīng)的極限強度；η為顆粒極限強度分布的尺寸效應(yīng)參數(shù)。

圖3為顆粒在加載過程中可能發(fā)生的破壞過程示意圖。

圖3 顆粒的可能破碎過程示意圖Fig.3 Possible crushing process of a particle

對任意顆粒，當(dāng)其受力滿足公式（1）時顆粒發(fā)生破碎，該顆粒沿最大接觸力方向劈裂為2個碎片。當(dāng)碎片顆粒再次達到破碎條件時劈裂為2個新的碎片，直到碎片顆粒僅由1個基本球組成時停止。新生成的碎片顆粒屬性與原顆粒一致，強度根據(jù)碎片顆粒的等效直徑重新指定[25-26]。

2 細觀參數(shù)的確定

細觀模擬需要確定的參數(shù)主要分為接觸模型相關(guān)參數(shù)和顆粒強度參數(shù)兩類。接觸模型參數(shù)為剪切模量G、泊松比v、摩擦系數(shù)μ。顆粒強度參數(shù)包括Weibull模量m、尺寸效應(yīng)參數(shù)η和極限強度σlim,0。細觀參數(shù)的確定流程是在參考前人細觀參數(shù)選取規(guī)則和成果的基礎(chǔ)上，確定細觀參數(shù)的大致范圍；然后，對不同的參數(shù)進行敏感性分析，了解各個參數(shù)對動力特性的影響；最后，通過微調(diào)得到一套可以反映研究對象變形特性的參數(shù)。根據(jù)楊貴等[11]進行的參數(shù)敏感性分析可知，細觀泊松比v在小于0.3的情況下，對試樣的動模量影響不大，參考其取值，細觀泊松比v取0.3。文獻[27]確定花崗巖顆粒強度相關(guān)的分布參數(shù)m和η。這樣需要確定的細觀參數(shù)為剪切模量G、摩擦系數(shù)μ和極限強度σlim,0。

首先，對剪切模量G、摩擦系數(shù)μ和極限強度σlim,0進行敏感性分析，3個參數(shù)對骨架曲線的影響如圖4所示。

圖4 參數(shù)敏感性分析Fig.4 Parameter sensitivity analysis

剪切模量直接影響試樣的變形能力。由圖4（a）可知：剪切模量越大，相同動應(yīng)變對應(yīng)的動應(yīng)力越大；但是隨著剪切模量增加，骨架曲線隨動應(yīng)變的增長趨勢變緩。圖4（b）為摩擦系數(shù)對骨架曲線的影響。由圖4（b）可知：摩擦系數(shù)越小，顆粒之間越容易發(fā)生相對滑動，在相同荷載作用的條件下試樣產(chǎn)生的變形越大；動應(yīng)變相同時，動應(yīng)力隨著摩擦系數(shù)的增加而增大。圖4（c）為極限強度對骨架曲線的影響。由圖4（c）可知：動應(yīng)變較小時，不同極限強度對應(yīng)的動應(yīng)力差異不大；隨著動應(yīng)變增加，極限強度越大的試樣其動應(yīng)力也越大。

根據(jù)參數(shù)敏感性分析結(jié)果，通過與2 200 kPa圍壓下的骨架曲線對比，可以大致確定細觀參數(shù)的取值范圍。通過進一步試算，并與2 200 kPa圍壓下的骨架曲線對比，可以得到一套合理的細觀參數(shù)，如表2所示。

表2 細觀模擬參數(shù)Tab. 2 Micro-parameters for simulation

為進一步驗證參數(shù)的合理性，采用表2中的細觀參數(shù)，模擬900和1 500 kPa圍壓下的骨架曲線，如圖5所示。模擬的動應(yīng)力-動應(yīng)變骨架曲線與室內(nèi)試驗結(jié)果具有較好的一致性，說明該模型可以基本再現(xiàn)糯扎渡Ⅱ區(qū)粗堆石料的動力變形特性。

圖5 模擬和試驗骨架曲線對比Fig.5 Comparison of simulation and experimental skeleton curves

3 顆粒破碎的影響

圖6為考慮顆粒破碎和不考慮顆粒破碎條件下，動應(yīng)變隨循環(huán)振次的變化規(guī)律。由圖6可知：試樣在循環(huán)振動過程中，均發(fā)生了比較明顯的塑性變形，這主要是因為初始偏應(yīng)力的存在使得動應(yīng)變曲線偏向壓縮側(cè)，累積塑性變形均為壓應(yīng)變，同一圍壓下，動應(yīng)變幅值和塑性變形均隨動應(yīng)力的增加而增大；動荷載相同時，圍壓大的試樣對應(yīng)的動應(yīng)變幅值小；相同應(yīng)力狀態(tài)下，考慮顆粒破碎的試樣動應(yīng)變幅值較不考慮顆粒破碎的試樣大，相應(yīng)的塑性變形也越大，說明顆粒破碎導(dǎo)致試樣整體剛度降低，使試樣產(chǎn)生更多不可恢復(fù)的變形。

圖6 動應(yīng)變隨循環(huán)振次的變化關(guān)系Fig.6 Curves of dynamic strain-cycle number

顆粒破碎程度采用Marsal破碎率衡量。整理振動過程中不同工況下試樣的顆粒破碎率，如圖7所示。由圖7可知：圍壓越大，固結(jié)階段顆粒破碎越嚴(yán)重；顆粒破碎率在初始加載的1/2周期內(nèi)顯著增加，隨后增長減緩，說明循環(huán)荷載作用下顆粒破碎主要發(fā)生在第1圈的前1/2周期內(nèi)；圍壓相同時，動應(yīng)力越大，顆粒破碎率也越大。

圖7 循環(huán)加載過程中的顆粒破碎率Fig.7 Particle crushing during cyclic loading process

對任意顆粒，配位數(shù)為與該顆粒接觸的周圍顆粒數(shù)量。Thornton[28]指出，配位數(shù)小于2的懸浮顆粒對集合體的力學(xué)特性貢獻較小，可以忽略。試樣的有效配位數(shù)定義為：

式中，Ntotal、Neff、N0和N1分別為總的顆粒數(shù)量、去除懸浮顆粒后的有效顆粒數(shù)、配位數(shù)為0和配位數(shù)為1的顆粒數(shù)，C為顆粒接觸總數(shù)。

圖8（a）為循環(huán)振動過程中不同動應(yīng)力條件下，有效配位數(shù)的變化規(guī)律。由圖8（a）可知：不考慮顆粒破碎條件下，動應(yīng)力相同時圍壓越大，試樣的有效配位數(shù)越多；圍壓相同時，動應(yīng)力越大，試樣的有效配位數(shù)越?。谎h(huán)過程中有效配位數(shù)會緩慢減小，且顆粒破碎也會導(dǎo)致有效配位數(shù)的減少；配位數(shù)的頻率為具有相同配位數(shù)的顆粒數(shù)量與有效顆粒數(shù)之比。整理圍壓2 200 kPa、動應(yīng)力為2 000 kPa時，振動前后的配位數(shù)的頻率分布，如圖8（b）所示，圖8（b）中振動前指試樣固結(jié)完成施加動荷載前的狀態(tài)。由圖8（b）可知，考慮破碎和不考慮破碎條件下配位數(shù)的頻率分布圖相似；振動前，考慮破碎和不考慮破碎，試樣的峰值配位數(shù)均為5.0；振動后，位于峰值配位數(shù)左側(cè)的頻率增加，右側(cè)降低，且發(fā)生顆粒破碎的試樣峰值配位數(shù)降低到4.0。說明振動和顆粒破碎均會導(dǎo)致試樣中力學(xué)不穩(wěn)定顆粒的增加，使試樣內(nèi)部的結(jié)構(gòu)性變差，且顆粒破碎量越大，對試樣的結(jié)構(gòu)性影響越顯著，最終導(dǎo)致試樣的剛度和承載能力降低。

圖8 有效配位數(shù)的演化曲線與配位數(shù)頻率分布Fig.8 Evolution curves of effective coordination number and frequency distribution of coordination number of samples before and after cyclic loading

4 孔隙率對動力參數(shù)的影響

為了研究孔隙率對動力特性的影響，將制樣時顆粒的摩擦系數(shù)設(shè)置為0，采用同樣的制樣方法制得最密實的試樣，孔隙率為0.52。施加與前文孔隙率為0.55的試樣相同的振動工況，記錄振動過程中的顆粒破碎率和配位數(shù)等參數(shù)。

圖9為孔隙率對動應(yīng)變、顆粒破碎率和配位數(shù)分布的影響。

圖9 孔隙率對動應(yīng)變、顆粒破碎和配位數(shù)分布的影響Fig.9 Influence of porosity on dynamic strain, particle crushing and the distribution of coordination number

由圖9（a）可知，孔隙率對動力變形特性影響顯著。在同樣的圍壓和動應(yīng)力條件下，孔隙率大的試樣具有更大的動應(yīng)變。振動過程中的顆粒破碎發(fā)展如圖9（b）所示，孔隙率小的試樣在相同應(yīng)力狀態(tài)下，顆粒破碎率較孔隙率大的試樣小。

整理動應(yīng)力為2 000 kPa振動前后配位數(shù)的頻率分布圖如圖9（c）所示。由圖9（c）可知：孔隙率為0.55和0.52的試樣頻率峰值對應(yīng)的配位數(shù)均為5；孔隙率為0.55的試樣峰值左側(cè)的頻率均高于孔隙率為0.52的試樣，峰值右側(cè)頻率均比孔隙率為0.52的試樣低；振動后頻率分布曲線的峰值有所下降，峰值左側(cè)上抬，右側(cè)下降，且孔隙率為0.55的試樣變化趨勢更顯著。說明同一級配的同種材料，密實試樣內(nèi)部的平均接觸數(shù)更多，在相同受力狀態(tài)下不易發(fā)生顆粒破碎，力學(xué)性能更好。

5 討 論

同一圍壓下的動應(yīng)變與動彈性模量可以近似用直線關(guān)系式（4）擬合[29]：

式中，a和b分別為擬合直線的截距和斜率，1/a為最大動彈性模量Edmax。斜率b越大，表明動彈性模量隨著動應(yīng)變的增加衰減越快。

最大動彈性模量與試樣的應(yīng)力狀態(tài)相關(guān)，可表示為：

式中：σm為平均有效固結(jié)主應(yīng)力，σm=（σ1+2σ3）/3；Pa為大氣壓力；k和n為試驗常數(shù)。

整理不同圍壓下的動彈性模量與動應(yīng)變值，繪制動彈性模量的倒數(shù)與動應(yīng)變的關(guān)系，采用式（4）擬合。圖10匯總了不同情況下，斜率b的分布規(guī)律。由圖10可知：斜率b均隨著平均有效固結(jié)主應(yīng)力的增加而減小，說明動彈性模量的衰減速率隨有效固結(jié)主應(yīng)力的增加變緩；平均有效固結(jié)主應(yīng)力相同時，考慮顆粒破碎的試樣較不考慮顆粒破碎的試樣具有更大的b值，動彈性模量衰減速率增加。不同孔隙率試樣對應(yīng)的斜率如圖10（b）所示。由圖10（b）可知，相同圍壓條件下小孔隙試樣的斜率均小于大孔隙試樣，說明疏松試樣的動彈性模量隨動應(yīng)變的增加衰減更快，這與不同密實程度卵石的動模量衰減規(guī)律一致[30]。

圖10 顆粒破碎和孔隙率對斜率b的影響Fig.10 Influence of particle breakage and porosity on the value of b

不同孔隙率試樣的最大動彈性模量與平均有效固結(jié)主應(yīng)力之間的關(guān)系如圖11（a）所示。小孔隙率試樣的最大動彈性模量比大孔隙率試樣大。平均有效固結(jié)主應(yīng)力相同時，孔隙率為0.52的試樣最大動彈性模量約為孔隙率0.55試樣的1.2倍。采用式（5）擬合，k值對應(yīng)于平均有效圍壓為Pa時的最大動彈性模量。孔隙率為0.55的試樣對應(yīng)的k值為（5 615），明顯小于孔隙率為0.52試樣的k值（7 028）。n值的差異不顯著。

Hardin等[31-32]認為材料的變形特性主要與孔隙比和應(yīng)力狀態(tài)相關(guān)，提出可以用經(jīng)驗公式（6）來描述變形模量與孔隙比和平均有效應(yīng)力之間的關(guān)系：

式中，A為材料常數(shù)，e為孔隙比，f(e)為孔隙比的函數(shù)。不同的研究者針對不同的巖土材料提出不同的表達形式用來消除孔隙率的影響，常用的函數(shù)形式為f(e)=(2.973-e)2/(1+e)。

采用函數(shù)f（e）對最大動彈性模量進行歸一化，如圖11（b）所示。由圖11（b）可知，不同孔隙不同圍壓下的最大動彈性模量可以采用統(tǒng)一的公式描述：Edmax=4085f(e)Pa(σm/Pa)0.34。

圖11 孔隙率對最大動彈性模量的影響Fig.11 Influence of porosity on the maximum dynamic elastic modulus

整體而言，顆粒破碎主要影響最大動彈性模量的衰減速率?？紫堵市〉脑嚇邮芰π阅芨?，同樣應(yīng)力條件下顆粒不易發(fā)生破碎，模量隨動應(yīng)變的衰減速率也相應(yīng)較慢。這對土石壩的抗震具有較好的借鑒意義，施工中可通過提高填筑密度來改善壩體的抗震性能。

6 結(jié) 論

采用PFC3D離散元程序?qū)ε丛散騾^(qū)粗堆石料循環(huán)荷載作用下的變形特性進行模擬。研究了小應(yīng)變條件下，顆粒破碎和孔隙率對動彈性模量的影響。主要結(jié)論如下：

1）不同圍壓下，數(shù)值模擬的骨架曲線與室內(nèi)堆石料的試驗結(jié)果吻合較好。圍壓越大，相同動應(yīng)變條件下試樣的動彈性模量越大。

2）循環(huán)振動和顆粒破碎均會降低集合體內(nèi)部的有效配位數(shù)。

3）孔隙率小的試樣動力特性要優(yōu)于孔隙率大的試樣。相同圍壓和動應(yīng)力條件下，孔隙率小的試樣有效配位數(shù)大，顆粒破碎量小。

4）顆粒破碎和孔隙率均影響動彈性模量的衰減速率。最大動彈性模量主要受平均有效應(yīng)力和孔隙率的影響，可以用Hardin等提出的經(jīng)驗公式描述孔隙率對最大動彈性模量的影響。