剪切波在堆石體中傳播的細觀數(shù)值模擬

楊尚慧 ，馬剛 ，周偉 ，陳遠 ，常曉林

(1. 武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢，430072；2. 武漢大學 水工巖石力學教育部重點實驗室，湖北 武漢 430072)

巖土顆粒材料的小應(yīng)變剪切模量是巖土工程設(shè)計與分析中的重要力學參數(shù)之一[1?2]，被廣泛應(yīng)用于地震、爆破等土動力問題，以及分析一些土工結(jié)構(gòu)物作用下地基土的靜力變形中。土體的小應(yīng)變剪切模量可通過現(xiàn)場試驗和室內(nèi)試驗來測定[3?5]，包括下孔法、跨孔法以及表面波法等。室內(nèi)試驗則以共振柱和彎曲元試驗為主，通過直接測定土體的剪切波速獲得。自SHIRLEY等[6]采用彎曲元法測試室內(nèi)制備高嶺土試樣的剪切波速和能量衰減特性，這項無損技術(shù)因原理明確、操作簡單、價格低廉等優(yōu)勢，被廣泛用于土體小應(yīng)變剪切模量的測定[7?8]。彎曲元通常是由緊密粘結(jié)的2片或多片壓電陶瓷晶體片組成，是一種可通過壓電效應(yīng)和逆壓電效應(yīng)實現(xiàn)機械能和電能之間轉(zhuǎn)換的機電傳感器[7]。進行土樣剪切波速測試時，固定在裝置一端的彎曲元元件作為激發(fā)元在脈沖電壓下產(chǎn)生振動，激發(fā)出剪切波，固定在裝置另一端的彎曲元元件作為接收元接收經(jīng)土體傳播而來的剪切波，并通過數(shù)據(jù)采集裝置將接收元產(chǎn)生的振動轉(zhuǎn)化為電信號，進而與激發(fā)元的電信號對比即可計算出土體的剪切波速。土體的剪切波速與眾多因素密切相關(guān)，如顆粒形狀、級配、孔隙比、應(yīng)力水平和應(yīng)力歷史等[9?11]。ZENG等[12]通過彎曲元法測試了砂土在無側(cè)向應(yīng)變條件下不同應(yīng)力面上的小應(yīng)變剪切模量，提出了不同應(yīng)力狀態(tài)下小應(yīng)變剪切模量的經(jīng)驗公式。吳宏偉等[13]利用彎曲元法研究了應(yīng)力狀態(tài)對上海原狀黏土小應(yīng)變剪切模量的影響，并進一步借助孔隙比函數(shù)建立起與土體應(yīng)力狀態(tài)的關(guān)系表達式。徐潔等[14]通過2組彎曲元試驗研究了吸力和干濕路徑對粉土小應(yīng)變剪切模量的作用，并結(jié)合試驗結(jié)果改進了小應(yīng)變剪切模量半經(jīng)驗公式。采用彎曲元技術(shù)進行的剪切波試驗時，如何準確地確定剪切波的傳播時間，一直存在較大的爭議。這是因為顆粒體系的內(nèi)在復(fù)雜性、信號干擾(如串音、壓縮波干預(yù)等)以及本身的組構(gòu)缺陷[15?18]，可能產(chǎn)生30%~50%的測量誤差[17]。目前，國內(nèi)外學者利用離散元法與室內(nèi)試驗相結(jié)合，針對顆粒材料的波動問題開展了大量研究，克服了物理試驗中的信號干擾問題，能夠得到質(zhì)量較高的剪切波接收信號[18?20]。徐曉敏等[18]針對無黏性顆粒材料進行了三軸試驗離散元模擬，并設(shè)置激發(fā)源和接收源實現(xiàn)了顆粒材料中剪切波傳播的模擬，系統(tǒng)分析了激發(fā)頻率等諸多因素對剪切波速的影響。O′DONOVAN等[19]模擬了均勻圓盤的理想六角堆積體的二維彎曲元試驗，比較了4種剪切波速的確定方法，并研究了離散元模型的細觀參數(shù)與剪切波速的關(guān)系。NING等[20]對由球形顆粒隨機堆積而成的圓柱體試樣進行了三維彎曲元試驗數(shù)值模擬，采用顆粒速度矢量場的演變表示波傳播過程中的顆粒運動情況。MA等[21?24]提出了隨機顆粒不連續(xù)變形方法，采用隨機模擬技術(shù)產(chǎn)生凸多邊形顆粒及其在空間中的分布，真實模擬了顆粒形態(tài)，將顆粒劃分成有限個單元，允許顆粒發(fā)生變形，顆粒與顆粒之間的相互作用力通過虛擬的法向和切向彈簧和阻尼器來傳遞。本文作者基于隨機顆粒不連續(xù)變形方法從堆石體顆粒位移演化以及頻域分析方面研究了剪切波傳播特性；采用3種方法確定剪切波速并與雙軸壓縮試驗結(jié)果進行比較，從而得到較優(yōu)的剪切波傳播時間確定方法；系統(tǒng)地分析激發(fā)頻率、激發(fā)幅值、顆粒接觸剛度以及圍壓對剪切波速的影響。

1 彎曲元數(shù)值試驗

1.1 數(shù)值試樣

生成二維多邊形顆粒及其在空間中的分布，再壓縮成指定尺寸的試樣，試樣級配曲線如圖1所示，孔隙率為0.135。試樣長×寬為12 m×6 m，最大粒徑dmax為 80 mm，dmax/dmin=4(其中，dmin為最小粒徑)。定義顆粒的外接橢圓的長短徑之比為顆粒的形狀指標，長短徑比在 1.4~1.6之間呈均勻分布。試樣中共生成 10 217個顆粒，采用二階四邊形單元離散得到64 838個單元和243 520個節(jié)點。為了模擬室內(nèi)彎曲元剪切波速測試裝置，試樣上下端為剛性板，四周用橡膠模包裹住，橡膠模上下端與剛性板分離，圍壓施加在橡膠模上。橡膠模的厚度取 2 mm，采用超彈性膜單元模擬橡膠模。

圖1 數(shù)值試樣級配曲線圖Fig. 1 Particle size distribution of numerical sample

1.2 剪切波產(chǎn)生與傳播

由于正弦波適用土樣范圍較廣，而方波僅適用于剛度較小的土樣[15]，本文采用正弦波作為激發(fā)信號。在試樣頂部選取直徑為0.75 m、厚度為0.20 m的薄層顆粒作為激發(fā)源，施加水平方向單周期正弦位移產(chǎn)生剪切波，與由單個顆粒作為激發(fā)源相比，薄層波源能夠有效地減少物理試驗中由壓縮波和剪切波相互耦合形成的近場效應(yīng)的干擾[15,20]。同時，沿試樣軸線選取3個顆粒(R1，R2，R3)作為接收源，距離試樣頂部的距離分別為6，9和12 m，并輸出顆粒位移時程曲線，顆粒位移取顆粒上所有節(jié)點位移的平均值。數(shù)值試樣如圖2所示，模型主要計算參數(shù)如表1所示，所選參數(shù)可以保證基本消除彌散特性(λ/d50≥10，λ為波長；d50為顆粒平均粒徑)與近場效應(yīng)的影響(L/λ＞2，L為傳播距離)[25?26]。

數(shù)值模擬開始時，先對試樣施加圍壓進行等向固結(jié)，固結(jié)完成后對試樣進行全約束繼續(xù)運行一定時步，使得顆粒速度與激發(fā)幅值相比足夠小，實現(xiàn)“除噪”。當試樣中顆粒趨于平衡，對薄層激發(fā)源顆粒施加水平方向單周期正弦位移產(chǎn)生剪切波，側(cè)向薄膜施加恒定圍壓，同時監(jiān)測接收源顆粒水平方位移和豎直向坐標隨時間變化動態(tài)。

圖2 數(shù)值試樣Fig. 2 Numerical sample

表1 計算參數(shù)Table 1 Parameters for calculation

2 數(shù)值模擬結(jié)果及分析

2.1 位移演化

圖3所示為時間t為0和2T(T為施加的正弦位移的周期)時刻堆石體顆粒水平方向位移云圖。從圖3可以看出：在初始時刻，僅薄層顆粒產(chǎn)生位移，隨著時間的推進，薄層顆粒的運動向下傳遞，且運動幅度逐漸降減弱，顯示了波在堆石體中的傳播過程。

圖3 不同時刻堆石體顆粒位移云圖Fig. 3 Contours of particle displacement for different time

圖4所示為數(shù)值試樣中激發(fā)源與接收源顆粒水平方向位移時程曲線，其中，S代表激發(fā)源顆粒，R1，R2和R3代表接收源顆粒。從圖4可以看出：接收源顆粒按照距離波源的遠近依次做出響應(yīng)，且波動幅度逐漸衰減；接收源顆粒 R1初始出現(xiàn)與激發(fā)源相位相反的位移，這是由于近場效應(yīng)壓縮波(P波)傳播速度大于剪切波(S波)速度，先期到達；隨著傳播距離的增加，近場效應(yīng)減弱，R2和R3顆粒初始反相位位移現(xiàn)象逐漸消失。對于單個接收源顆粒而言，運動經(jīng)歷了初期衰減隨后在0.04~0.05 s間再次增強的過程，這是由底部剛性板邊界的反射造成的。綜合考慮剪切波傳播過程中的近場效應(yīng)以及由于阻尼等因素引起的振幅衰減所造成的影響，選取接收源顆粒 R1位移隨時間變化曲線作為確定剪切波速的依據(jù)。

2.2 頻域分析

頻域分析方法又被稱為快速傅里葉變換(fast Fourier transform, FFT)分析法，通過對獲取的波信號數(shù)據(jù)進行快速傅里葉變換，將時域信號轉(zhuǎn)變?yōu)轭l域信號，并對其時域和頻域特性進行分析。為了更加精確地比較發(fā)射信號和接收信號的各頻率組分，分別對其進行快速傅里葉變換，由于接收信號振幅衰減明顯，與發(fā)射信號相差較大，不便于比較，因此，將接收信號振幅擴大20倍，結(jié)果如圖5所示。

圖4 輸入接收信號波形圖Fig. 4 Signals transmitted and received by bender element

圖5 0.8 kHz發(fā)射信號與接收信號頻域?qū)Ρ菷ig. 5 Comparison between input signals and received signals in frequency domain for 0.8 kHz

由圖5可知：堆石體顆粒對高頻率剪切波信號具有明顯的過濾作用。對于發(fā)射信號而言，單周期正弦波具有較寬的頻率范圍，最高頻率達到1.6 kHz附近，峰值振幅對應(yīng)的頻率為0.8 kHz，而接收信號所包含的頻率范圍較窄，主要集中在1 kHz以下，最高頻率大約在1.1 kHz附近；當頻率大于0.5 kHz時，振幅迅速衰減，峰值振幅對應(yīng) 3組頻率分別為 0.12，0.15和0.18 kHZ。改變發(fā)射信號的頻率，分別進行正弦位移頻率為0.4，0.6和1 kHz的數(shù)值試驗，并對發(fā)射信號和接收信號進行快速傅里葉變換，結(jié)果如圖6所示。

由圖5和圖6可知：隨著發(fā)射信號頻率的增加，其所包含的頻率范圍也隨之增大，峰值振幅對應(yīng)的頻率分別為0.4，0.6和1.0 kHz，而接收信號的頻率范圍與發(fā)射信號的頻率無關(guān)，主要集中在1.0 kHz以下部分，并且當接收信號的頻率大于0.5 kHz時，其對應(yīng)的振幅將迅速減小。定義位移振幅趨于0 mm時對應(yīng)的頻率為堆石體振幅衰減的臨界頻率，計算得到不同發(fā)射信號頻率下振幅衰減的臨界頻率的方差為3.2×10?5，說明其隨發(fā)射信號頻率變化不明顯，證明臨界頻率是與試驗參數(shù)無關(guān)的堆石體自身的特性參數(shù)。

圖6 發(fā)射信號與接收信號頻域?qū)Ρ菷ig. 6 Comparison between input signals and received signals in frequency domain for different input signals

3 剪切波傳播時間的確定

剪切波速由發(fā)射源和激發(fā)源的距離和剪切波的傳遞時間來確定，目前有2類確定剪切波傳播時間的方法：直接確定法和計算分析法。直接確定法根據(jù)輸入波形和輸出波形的特征點確定傳播時間，常見的有初達波法(S-S法)和峰值-峰值法(P-P法)，計算法通過互相關(guān)分析確定傳播時間。分別根據(jù)S-S法、P-P法和互相關(guān)分析法薄層顆粒S和R1接收源顆粒位移波形曲線(見圖7)確定顆粒體的剪切波速，并與雙軸壓縮試驗結(jié)果進行對比，從而選擇出較優(yōu)的剪切波傳播時間確定方法。

圖7 剪切波傳播時間的確定Fig. 7 Determination of travel time of shear wave

3.1 初達波法(S-S法)

初達波法利用剪切波的首次到達來確定傳遞時間，通常是以接收信號的首個偏轉(zhuǎn)點作為剪切波的到達位置。發(fā)射彎曲元在激發(fā)剪切波的同時也產(chǎn)生壓縮波，由于壓縮波在土體中的傳播速度比剪切波的大，將會引起接收信號產(chǎn)生初始位移而隨后到達的剪切波將會改變波信號的方向，因此，將A′點視為剪切波到達接收源顆粒R1的初始位置，對應(yīng)的時間為 18.288 ms，此時顆粒與薄層底部的距離為5.780 63 m，計算得到剪切波速為316.09 m/s。

3.2 特征點法(P-P法)

特征點法分別選取發(fā)射信號和接收信號的首個峰值點，通過峰值點間的時間和位移差來計算剪切波速，由圖7可知：接收源顆粒R1位移波形的第1峰值B′對應(yīng)的時間為20.603 ms，此時顆粒與薄層底部的距離為5.780 59 m，計算得到剪切波速為280.57 m/s。與初達波法相比，特征點法中接收信號的峰值點由于不受噪音干擾能夠被清楚識別，但對選擇首個峰值點還是最大峰值點問題仍存在一定爭議[27]。

3.3 互相關(guān)法

初達波法和特征點法中均需要人為判斷剪切波到達和傳播時間，具有一定的主觀性，而互相關(guān)法是一種數(shù)學方法，通過互相關(guān)函數(shù)即CXY來衡量2個信號之間的相關(guān)性，其數(shù)學表達式為

式中：T1為信號總測試時間；ξ為2個信號之間的時間差。

對發(fā)射信號和接收信號進行快速傅里葉變換得到如圖8所示的互相關(guān)函數(shù)，取CXY最大值對應(yīng)的時間tc作為剪切波的傳播時間，得到tc=19.015 ms，計算得到剪切波速為299.91 m/s。

圖8 互相關(guān)法確定剪切波到達時間Fig. 8 Determination of arrive time using cross-correlation method

3.4 剪切波速對比

進行雙軸壓縮試驗得到小應(yīng)變剪切模量Gmax，并進一步根據(jù)Gmax=ρvs2(其中，ρ為顆粒密度；vs為剪切波速)，計算剪切波速vs作為參照。雙軸壓縮試驗與彎曲元試驗采用相同的邊界條件，均保持恒定圍壓，并將頂部剛性板按照1 kHz頻率施加6級往復(fù)荷載，根據(jù)由6個循環(huán)荷載級的軸向動應(yīng)力幅值σd和軸向動應(yīng)變幅值εd繪制的σd?εd曲線的初始斜率計算得到數(shù)值試樣的動彈性模量Emax=448.38 MPa，根據(jù)公式Gmax=Emax/2(1+v)最終得到vs的參考值為292.93 m/s。

將3種方法計算得到的剪切波速與雙軸壓縮試驗參考值進行比較，如表2所示。從表2可以看出：采用互相關(guān)法計算得到的剪切波速與雙軸試驗結(jié)果最接近，說明基于隨機顆粒不連續(xù)變形SGDD方法的堆石體二維彎曲元試驗的數(shù)值模擬中，采用互相關(guān)法確定的剪切波傳播時間精確度最高，與AHN等[28]基于離散元方法得到的結(jié)論相一致，因此，本文選擇互相關(guān)分析作為計算剪切波傳播時間的方法。

表2 不同計算方法剪切波速對比Table 2 Comparison of values of vs calculated by different methods

4 剪切波傳播影響因素分析

堆石體中剪切波的傳播受諸多因素的影響，主要包括激發(fā)幅值、激發(fā)頻率、圍壓以及細觀層次上的顆粒間接觸剛度。

4.1 激發(fā)幅值

由于堆石體的剪切剛度是隨應(yīng)變水平變化的函數(shù)，因此剪切波速在不同的應(yīng)變水平下有所不同。彎曲元試驗中試樣產(chǎn)生的剪切應(yīng)變則是通過激發(fā)幅值來控制。分別進行A/d50=1×10?6~1×10?2(其中，A代表激發(fā)幅值)的6組平行試驗，計算得到6種情況下的剪切波速如圖9所示，同時提取出剪切波傳播過程中接觸失效率隨時間變化曲線如圖10所示。由圖9可知：隨著A/d50減少，剪切波速略有增加，當激發(fā)幅值小于10?3時，剪切波速趨于收斂，基本穩(wěn)定。對應(yīng)圖10中的接觸失效率變化曲線，當A/d50小于 10?3時，顆粒接觸失效率的變化曲線近似重合，而當A/d50大于10?3時，接觸失效率曲線上升，對應(yīng)的剪切波速略有降低。

圖9 激發(fā)幅值對剪切波速的影響Fig. 9 Influence of excitation amplitude on shear wave velocity

圖10 不同激發(fā)幅值下接觸失效率Fig. 10 Percentage of sliding contacts for different excitation amplitudes

4.2 激發(fā)頻率

圖11所示為激發(fā)頻率分別為0.4，0.6，0.8和1.0 kHz時接收源顆粒R1輸出的位移波形圖。從圖11可以看出：當激發(fā)頻率小于0.6 kHz時，接收源顆粒的位移波形較微弱；當激發(fā)頻率達到0.8 kHz后，信號強度增強并且波形相似，但1.0 kHz下波形振幅略有下降，即當激發(fā)頻率為0.8 kHz時，接收信號振幅達到峰值，對應(yīng)堆石體顆粒以最大振幅做振動。又由圖5和圖 6可知：接收信號包含的頻率主要集中在0.8 kHz以下，因此，可近似認為0.8 kHz為堆石體的共振頻率。另外，從圖11還可以直接觀察到不同頻率下剪切波的初達時間，從而確定剪切波速，進一步得出當激發(fā)頻率從0.4 kHz增加到0.8 kHz時，剪切波速略有增大；當激發(fā)頻率達到共振頻率后，剪切波速將保持穩(wěn)定。

堆石體的強度部分來源于堆石體顆粒間摩擦的激勵，為研究堆石體顆粒接觸處的摩擦特性，定義堆石體顆粒摩擦激勵指標為Im=|ft|/(fntanφ)(其中，ft和fn分別表示堆石體顆粒間切向接觸力和法向接觸力)。圖12所示為不同激發(fā)頻率下堆石體平均摩擦激勵指標Im的演化曲線(Im為相應(yīng)積分點對應(yīng)的堆石體中所有接觸處的平均摩擦激勵)。從圖12可以看出，初始階段0.4 kHz和0.6 kHz激發(fā)頻率下堆石體的摩擦激勵系數(shù)演化曲線位于0.8 kHz和1 kHz曲線下方，當激發(fā)頻率大于0.8 kHz時，曲線近似重合，對應(yīng)圖11中當激發(fā)頻率小于共振頻率時，堆石體響應(yīng)較微弱，當激發(fā)頻率達到共振頻率后，響應(yīng)增強并且波形相似。同時，隨著波的傳播，不同激發(fā)頻率下堆石體的摩擦激勵系數(shù)最終穩(wěn)定在同一數(shù)值附近，對應(yīng)堆石體的剪切波速隨激發(fā)頻率波動較小，穩(wěn)定在1個定值附近。

圖11 不同激發(fā)頻率下接收波形Fig. 11 Waveform of a receiver in the time domain,excited by different frequencies

圖12 不同激發(fā)頻率下摩擦激勵系數(shù)演化過程Fig. 12 Evolution of friction incentive index, excited by different frequencies

4.3 圍壓

堆石體的應(yīng)力狀態(tài)通過影響堆石體顆粒間的剛度從而使剪切波速發(fā)生改變，分別計算堆石體顆粒在0.4，0.5，0.6，0.8和1.0 MPa圍壓固結(jié)下的剪切波速，結(jié)果如圖13所示。從圖13可以看出：隨著圍壓的增加，剪切波速逐漸增大，且低圍壓下的增加速度略大于高圍壓下的增加速度，即隨著體系壓強的增大，顆粒之間的接觸更加緊密，對應(yīng)的平均接觸數(shù)增加，顆粒體系內(nèi)部結(jié)構(gòu)更不容易發(fā)生調(diào)整，相應(yīng)的剪切波速增加速度降低。這與徐小敏等[18]得到的剪切波速隨圍壓變化的基本規(guī)律一致，但剪切波速增長速度的變化仍存在差異，造成的原因可能是堆積方式以及填充率不同。

圖13 圍壓對剪切波速的影響Fig. 13 Influence of confining pressure on shear wave velocity

4.4 顆粒接觸剛度

堆石體的細觀數(shù)值試驗采用基于罰函數(shù)的線性接觸模型，即通過法向和切向彈簧來傳遞堆石體顆粒間的法向和切向相互作用力，因此，需要堆石體顆粒間接觸剛度作為細觀參數(shù)。而在該接觸模型中，接觸剛度并沒有明確的物理意義，也很難通過常規(guī)試驗獲取，因此，本文計算堆石體的法向接觸剛度Kn分別取值1×105，5×105，1×106，5×106和 1×107kN/m3時的剪切波速，對其進行敏感性分析，如圖14所示。由圖14可知：剪切波速隨法向接觸剛度的增大而增加，并與法向接觸剛度的開方近似成直線關(guān)系。DONOVAN等[19]采用離散元法模擬二維彎曲元試驗也得到了相似的剪切波速與顆粒接觸剛度之間的關(guān)系。

圖14 顆粒接觸剛度對剪切波速的影響Fig. 14 Influence of inter-particle contact spring stiffness on shear wave velocity

5 結(jié)論

1)堆石體顆粒對高頻率的剪切波信號具有過濾作用，接收信號相對于發(fā)射信號而言包含的頻率范圍變窄，主要集中在1 kHz以下，并且振幅衰減的臨界頻率是與試驗參數(shù)無關(guān)的堆石體自身的特性參數(shù)。

2)與初達波法和特征點法相比，互相關(guān)法計算得到的剪切波速與雙軸壓縮試驗結(jié)果最接近，進一步驗證了將其作為確定剪切波傳播時間方法的合理性。

3)堆石體顆粒中剪切波的傳播受諸多因素的影響，其中，隨著激發(fā)幅值減少，剪切波速略有增加，當A/d50小于10?3時，剪切波速趨于收斂，基本穩(wěn)定；當激發(fā)頻率小于共振頻率時，接收源顆粒的位移波形較微弱，且剪切波速隨頻率略有增加，當激發(fā)頻率達到共振頻率后，信號強度增強并且波形相似，同時剪切波速將保持穩(wěn)定；隨圍壓增加，剪切波速逐漸增大；剪切波速隨法向接觸剛度的增大而增加，并與法向接觸剛度的開方近似成直線關(guān)系。

[1]CLAYTON C R I. Stiffness at small strain: research and practice[J]. Géotechnique, 2011, 61(1): 5?37.

[2]SANTAMARINA J C, KLEIN A, FAM M A. Soils and waves:Particulate materials behavior, characterization and process monitoring[J]. Journal of Soils and Sediments, 2001,1(2): 130?130.

[3]KOKUSHO T. Cyclic triaxial test of dynamic soil properties for wide strain range[J]. Soils & Foundations, 1980, 20(2): 45?60.

[4]BENEDETTO H D, EZAOUI A. Experimental measurements of the global anisotropic elastic behaviour of dry Hostun sand during triaxial tests, and effect of sample preparation[J].Géotechnique, 2009, 59(7): 621?635.

[5]CHAUDHARY S K. Measurement of quasi-elastic stiffness parameters of dense Toyoura sand in hollow cylinder apparatus and triaxial apparatus with bender elements[J]. Geotechnical Testing Journal, 2004, 27(1): 1?13.

[6]SHIRLEY D J, HAMPTON L D. Shear-wave measurements in laboratory sediments[J]. The Journal of the Acoustical Society of America, 1978, 63(2): 607?613.

[7]陳云敏, 周燕國, 黃博. 利用彎曲元測試砂土剪切模量的國際平行試驗[J]. 巖土工程學報, 2006, 28(7): 874?880.CHEN Yunmin, ZHOU Yanguo, HUANG Bo. International parallel test on the measurement of shear modulus of sand using bender elements[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(7): 874?880.

[8]FONSECA A V D, FERREIRA C, FAHEY M. A framework interpreting bender element tests, combining time-domain and frequency-domain methods[J]. Geotechnical Testing Journal,2009, 32(2): 91?107.

[9]周燕國. 土結(jié)構(gòu)性的剪切波速表征及對動力特性的影響[D].杭州: 浙江大學建筑工程學院, 2007: 95?116.ZHOU Yanguo. Shear wave velocity-based characterization of soil structure and its effects on dynamic behavior[D]. Hangzhou:Zhejiang University. College of Civil Engineering and Architecture, 2007: 95?116.

[10]HARDIN B O, RICHART F E. Elastic wave velocities in granular soils[J]. Journal of the Soil Mechanics & Foundations Division, 1963, 89: 33?66.

[11]SUITS L D, SHEAHAN T C, PATEL A, et al. An empirical relationship for determining shear wave velocity in granular materials accounting for grain morphology[J]. Geotechnical Testing Journal, 2009, 32(1): 1?10.

[12]ZENG X, NI B. Stress-induced anisotropic G max of sands and its measurement[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1999, 125(9): 741?749.

[13]吳宏偉, 李青, 劉國彬. 利用彎曲元測量上海原狀軟黏土各向異性剪切模量的試驗研究[J]. 巖土工程學報, 2013, 35(1):150?156.WU Hongwei, LI Qing, LIU Guobin. Measurements of small-strain inherent stiffness anisotropy of intact Shanghai soft clay using bender elements[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(1): 150?156.

[14]徐潔, 周超. 干濕路徑影響粉土小應(yīng)變剪切模量的試驗研究[J]. 巖土力學, 2015: 36(1): 577?381.XU Jie, ZHOU Chao. Experimental study of effect of wetting-drying path on small-strain shear modulus of silt[J].Rock and Soil Mechanics, 2015: 36(1): 577?381.

[15]LEE J S, SANTAMARINA J C. Bender elements: performance and signal interpretation[J]. Journal of geotechnical and geoenvironmental engineering, 2005, 131(9): 1063?1070.

[16]ALVARADO G, COOP M R. On the performance of bender elements in triaxial tests[J]. Géotechnique, 2012, 62(1): 1?17.

[17]ARROYO M, MUIR WOOD D, GREENING P D. Source near-field effects and pulse tests in soil samples[J]. Géotechnique,2003, 53(3): 337?345.

[18]徐小敏, 凌道盛, 黃博, 等. 離散元模擬中顆粒材料剪切波速的剪切振動確定方法[J]. 巖土工程學報, 2011, 33(9):1462?1468.XU Xiaomin, LING Daosheng, HUANG Bo, et al.Determination of shear wave velocity in granular materials by shear vibration within discrete element simulation[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(9): 1462?1468.

[19]O’DONOVAN J, O’SULLIVAN C, MARKETOS G.Two-dimensional discrete element modelling of bender element tests on an idealised granular material[J]. Granular Matter, 2012,14(6): 733?747.

[20]NING Z, KHOUBANI A, EVANS T M. Shear wave propagation in granular assemblies[J]. Computers & Geotechnics, 2015, 69:615?626.

[21]MA G, ZHOU W, CHANG X L, et al. Combined fem/dem modeling of triaxial compression tests for rockfills with polyhedral particles[J]. International Journal of Geomechanics,2014, 14(4): 04014014.

[22]MA G, ZHOU W, CHANG X L. Modeling the particle breakage of rockfill materials with the cohesive crack model[J].Computers & Geotechnics, 2014, 61(61): 132?143.

[23]MA G, ZHOU W, NG T T, et al. Microscopic modeling of the creep behavior of rockfills with a delayed particle breakage

model[J]. Acta Geotechnica, 2015, 10(4): 481?496.

[24]MA G, ZHOU W, CHANG X L, et al. Formation of shear bands in crushable and irregularly shaped granular materials and the associated microstructural evolution[J]. Powder Technology,2016, 301: 118?130.

[25]TOOMEY A, BEAN C J. Numerical simulation of seismic waves using a discrete particle scheme[J]. Geophysical Journal International, 2000, 141(3): 595?604.

[26]CHANEY R C, DEMARS K R, PENNINGTON D S, et al.Horizontally mounted bender elements for measuring anisotropic shear moduli in triaxial clay specimens[J]. Geotechnical Testing Journal, 2001, 24(2): 133?144.

[27]LEONG E C. Measuring shear wave velocity using bender elements[J]. Astm Geotechnical Testing Journal, 2005, 28(5):488?498.

[28]AHN J, BISCONTIN G, ROESSET J M. Wave propagation in nonlinear one-dimensional soil model[J]. International Journal for Numerical & Analytical Methods in Geomechanics, 2010,33(33): 487?509.