劉曉宇 張國(guó)華 孫其誠(chéng) 趙雪丹 劉尚

1)(北京科技大學(xué)物理系,北京 100083)

2)(清華大學(xué),水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100084)

二維圓盤(pán)顆粒體系聲學(xué)行為的數(shù)值研究?

劉曉宇1)張國(guó)華1)?孫其誠(chéng)2)?趙雪丹1)劉尚1)

1)(北京科技大學(xué)物理系,北京 100083)

2)(清華大學(xué),水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100084)

顆粒物質(zhì),聲速,非線性,聲衰減

1 引 言

在非熱條件下壓縮顆粒系統(tǒng)可使系統(tǒng)經(jīng)歷堵塞(jamming)轉(zhuǎn)變[1?3].根據(jù)O’Hern的分析,在jamming點(diǎn)(J點(diǎn))附近,剪切模量G、體積模量B及剪切模量與體積模量的比值G/B均隨壓強(qiáng)P呈冪律標(biāo)度[3?5]:P～(Δ?)a?1,B～(Δ?)a?2,G～(Δ?)a?3/2,Z?Zc～(Δ?)1/2(Z為配位數(shù),Zc為J點(diǎn)配位數(shù)),G/B～P1/(2a?2),其中,a=2對(duì)應(yīng)線彈性接觸,a=5/2對(duì)應(yīng)赫茲接觸.因此,在J點(diǎn)附近剪切和體積模量比值G/B隨壓強(qiáng)P的標(biāo)度是驗(yàn)證jamming臨界性的重要標(biāo)記[6].最近,Wang等[7]發(fā)現(xiàn)當(dāng)P＞Pj(Pj為類(lèi)jamming轉(zhuǎn)變點(diǎn))時(shí),雙分散的赫茲接觸顆粒體系中G/B～P1/3,體系處于TL玻璃態(tài)(橫?？v模都存在的狀態(tài));當(dāng)P＜Pj時(shí),G/B～(P/T)?1/2,體系處于只能承載有效縱向聲子的L玻璃態(tài);當(dāng)P=0時(shí),體系處于無(wú)剛性的未堵塞態(tài).但是,在2D線彈性接觸的顆粒體系中是否也存在處于TL玻璃態(tài)和未堵塞態(tài)之間的L玻璃態(tài)仍然是一個(gè)未解決的問(wèn)題.另外,由于L玻璃態(tài)不能?chē)?yán)格定義橫波聲子,L玻璃態(tài)表現(xiàn)出與TL玻璃態(tài)非常不同的力學(xué)特性,目前關(guān)于L玻璃態(tài)力學(xué)特性的研究已經(jīng)成為一個(gè)新的熱點(diǎn).

顆粒物質(zhì)是典型的無(wú)序系統(tǒng),內(nèi)部本征的不均勻性使得顆粒系統(tǒng)的力學(xué)和輸運(yùn)特性與晶體的情況存在明顯不同[8].比如,顆粒材料中的振動(dòng)自由度不是總能用普通的平面波表征,而是存在一個(gè)在基本成分(顆粒)的微觀尺寸和宏觀樣品尺寸之間的特征介觀長(zhǎng)度尺度l?,低于l?的系統(tǒng)不再支持波長(zhǎng)λ小于l?的聲學(xué)平面波,連續(xù)體彈性理論失效[9].但是,λ?l?的長(zhǎng)波長(zhǎng)聲子仍然能在無(wú)序固體中存在,在這個(gè)長(zhǎng)度尺度上體系能被看作連續(xù)彈性介質(zhì).從這個(gè)角度講,顆粒材料等無(wú)序固體可以看作一系列由特征尺寸l?量級(jí)的區(qū)域結(jié)構(gòu)單元組成[10].近年來(lái),關(guān)于顆粒材料中介觀長(zhǎng)度尺度隨體積分?jǐn)?shù)的變化的研究已經(jīng)取得了很大的進(jìn)展.例如:文獻(xiàn)[11—13]通過(guò)研究振動(dòng)態(tài)密度曲線發(fā)現(xiàn)顆粒體系中存在區(qū)分德拜標(biāo)度和態(tài)密度平臺(tái)的過(guò)渡頻率ω?[5],與ω?對(duì)應(yīng)的縱波特征波長(zhǎng)vl為縱波聲速)及橫波特征波長(zhǎng)vt為橫波聲速)均隨到J點(diǎn)的距離(???c)發(fā)散:Vitelli[14]數(shù)值研究了通過(guò)排斥彈簧互相作用的無(wú)摩擦軟球jammed體系剪切聲波的長(zhǎng)波衰減,發(fā)現(xiàn)特征長(zhǎng)度??～(???c)?1/2.Mizuno等[15]數(shù)值研究發(fā)現(xiàn)3D無(wú)序體系偏離宏觀彈性性質(zhì)的特征尺度ξG～(?? ?c)?0.16.Wang等[7]數(shù)值確定顆粒材料等無(wú)序固體的特征長(zhǎng)度??～P?2/9.同時(shí)Ikeda和Berthier[6]的研究給出了不同的結(jié)論,即縱向和橫向振動(dòng)都在硬球區(qū)域中存在,且相關(guān)的長(zhǎng)度尺度和時(shí)間尺度均隨著密度減少遠(yuǎn)離J點(diǎn)而變微觀.由于對(duì)無(wú)序材料中的特征尺度的研究還存在爭(zhēng)議,值得進(jìn)一步研究.

聲學(xué)測(cè)量是一種探測(cè)顆粒系統(tǒng)內(nèi)部結(jié)構(gòu)和力學(xué)特性的有效手段[16?19].目前,人們已經(jīng)利用聲學(xué)手段對(duì)3維(3D)顆粒體系低壓下聲學(xué)特性進(jìn)行了大量研究,但是關(guān)于2D圓盤(pán)顆粒體系低壓強(qiáng)下是否也存在L玻璃態(tài)以及其低壓強(qiáng)下的聲學(xué)衰減、聲學(xué)非線性等行為的研究還較少.本文利用數(shù)值模擬手段對(duì)卸載過(guò)程中2D顆粒圓盤(pán)體系在低壓強(qiáng)下的聲學(xué)行為展開(kāi)了系統(tǒng)研究.首先,利用飛行時(shí)間法數(shù)值測(cè)量了卸載過(guò)程中不同壓強(qiáng)下二維單分散顆粒體系中的橫波聲速vt和縱波聲速vl,進(jìn)一步研究了剪切模量G與體積模量B的比值G/B隨壓強(qiáng)的變化規(guī)律,結(jié)果表明,在低壓強(qiáng)區(qū)域2D單分散圓盤(pán)顆粒體系處于L玻璃態(tài).其次,利用傅里葉變換研究了非線性系數(shù)、聲學(xué)衰減系數(shù)隨壓強(qiáng)P的變化以及衰減系數(shù)α隨頻率f的變化規(guī)律,并且進(jìn)一步研究了2D圓盤(pán)顆粒系統(tǒng)中散射相關(guān)的特征長(zhǎng)度??隨壓強(qiáng)的變化規(guī)律.本文的結(jié)果對(duì)于進(jìn)一步理解2D顆粒體系中是否存在L玻璃態(tài)以及L玻璃態(tài)的聲學(xué)衰減及非線性行為的規(guī)律有所幫助.

2 數(shù)值模擬

研究的顆粒體系由2048個(gè)質(zhì)量相同的無(wú)摩擦單分散圓盤(pán)顆粒組成,它們被隨機(jī)放置在邊長(zhǎng)為L(zhǎng)的正方形盒子中,采用周期性邊界,且不考慮重力的影響.顆粒與顆粒間的相互作用為單邊線性彈簧,即當(dāng)顆粒i和j間距rij=ri?rj小于它們的半徑之和時(shí)存在相互作用,其中ri,rj分別為顆粒i和j的位置矢量.接觸力的法向分量由給出,其中,kn是法向接觸剛度,ζ為阻尼比,nij=rij/|rij|是顆粒i中心指向顆粒j中心的單位矢量,mi和mj是顆粒i和顆粒j的質(zhì)量.在本文的模擬中,顆粒法向剛度系數(shù)為1.0×106N/m,切向剛度系數(shù)為0.

非熱(T=0)顆粒體系的具體制備過(guò)程如下:首先,在邊長(zhǎng)為L(zhǎng)的正方形盒子中隨機(jī)生成2048個(gè)半徑很小的單分散光滑圓盤(pán)顆粒.接著,一步增大顆粒半徑使體系的體積分?jǐn)?shù)達(dá)到某一固定值(P/T=∞).然后,采用共軛梯度使得體積能量趨于能量極小的亞穩(wěn)態(tài),可認(rèn)為體系為非熱體系(T=0).最后,通過(guò)逐漸減小顆粒半徑的方法對(duì)體系進(jìn)行卸載,生成一系列壓強(qiáng)的穩(wěn)定顆粒構(gòu)形. 本文分別取σ,m,ε/σ2和（其中,σ=Ri+Rj為顆粒的平均直徑,m為顆粒的質(zhì)量,為體系特征能量,Ri與Rj分別為體系中顆粒i與顆粒j的半徑為長(zhǎng)度、質(zhì)量、壓強(qiáng)和頻率的單位.

采用飛行時(shí)間(TOF)法測(cè)量顆粒體系的聲速,測(cè)量過(guò)程中分別使用單個(gè)顆粒作為激勵(lì)源和接收器.具體的測(cè)量過(guò)程如下:用單個(gè)正弦脈沖信號(hào)驅(qū)動(dòng)激勵(lì)源顆粒,使它從初速度為零開(kāi)始在水平(或垂直)方向上做一個(gè)完整的正弦振動(dòng)(v=0.0001ωsin(ωt),時(shí)間步長(zhǎng)t0=Tperiod/1000,Tperiod為周期,反比于圓頻率ω),同時(shí)記錄激勵(lì)和接收顆粒的v-t曲線.測(cè)量發(fā)射端信號(hào)、接收端信號(hào)第一個(gè)波峰的時(shí)間差tTOF=t2?t1(其中,t1和t2分別是激勵(lì)脈沖和接受脈沖的第一個(gè)波峰對(duì)應(yīng)的時(shí)間).進(jìn)而可以得到聲速vTOF=l/tTOF(l為激勵(lì)源和接收器顆粒間的距離).

采用傅里葉變換法來(lái)研究顆粒物質(zhì)中的聲學(xué)衰減和非線性現(xiàn)象.實(shí)驗(yàn)過(guò)程中,用連續(xù)的正弦信號(hào)激勵(lì)顆粒,并同時(shí)記錄激勵(lì)顆粒和接受顆粒的振動(dòng)速度隨時(shí)間的變化,如圖1所示.通過(guò)對(duì)接收和發(fā)射信號(hào)的傅里葉變換分析,從中提取出接收信號(hào)和發(fā)射信號(hào)的基頻、二倍頻幅值.可進(jìn)一步得到聲衰減系數(shù)其中,l為接收顆粒與激勵(lì)顆粒之間的距離,I0為激勵(lì)顆粒的振動(dòng)強(qiáng)度即初始強(qiáng)度,I(l)為接收顆粒的振動(dòng)強(qiáng)度即通過(guò)介質(zhì)顆粒衰減后的強(qiáng)度),及接收信號(hào)二倍頻幅值與接收信號(hào)基頻幅值平方的比值,進(jìn)而研究顆粒物質(zhì)的非線性特性.

圖1 激勵(lì)波形速度度幅A1(黑線)和接收波形的速度振幅A2(藍(lán)線)隨時(shí)間t的變化Fig.1.Time t dependence of velocity amplitude of exciting waveform A1(black line)and velocity amplitude of receiving waveform A2(blue line).

3 結(jié)果與討論

3.1 橫波、縱波聲速及材料剪切模量與體積模量比值G/B隨壓強(qiáng)的變化

圖2(a)和圖2(b)分別顯示了用飛行時(shí)間法測(cè)得的2D圓盤(pán)顆粒體系的橫波聲速vt和縱波聲速vl隨壓強(qiáng)的變化曲線.由圖2可知,vt(P)和vl(P)曲線均呈現(xiàn)分段行為:當(dāng)P＜10?4時(shí),橫縱波聲速均隨壓強(qiáng)的增加冪律減小,vt～P?0.446,vl～P?0.232;當(dāng)P＞10?4時(shí),橫波、縱波聲速均隨壓強(qiáng)增加冪律增大,vt～P0.202,vl～P0.338,與線彈性接觸理論[5,20]預(yù)言vt～P1/4一致.應(yīng)該指出的是,關(guān)于3D球形顆粒體系vt-P和vl-P曲線的分段行為在Somfai等[21]的數(shù)值模擬、Jia等[22]及張攀等[23]的實(shí)驗(yàn)中也有報(bào)道,Lherminier等[24]在研究聲波通過(guò)受壓的二維顆粒介質(zhì)時(shí)也發(fā)現(xiàn)了類(lèi)似現(xiàn)象.值得一提的是,圖2給出的vt(P)和vl(P)曲線的分段行為與理論的預(yù)言并不一致,即在本文的體系中vt(P)和vl(P)曲線呈現(xiàn)先冪律減小再冪律增大的規(guī)律.這可能對(duì)應(yīng)著卸載過(guò)程中在振動(dòng)源顆粒和接收顆粒之間形成了一個(gè)強(qiáng)力鏈,導(dǎo)致測(cè)得的橫波和縱波聲速都按冪律增加.vt(P)和vl(P)曲線的拐點(diǎn)實(shí)際上對(duì)應(yīng)著在振動(dòng)源顆粒和接收顆粒之間強(qiáng)力鏈的形成點(diǎn).關(guān)于整體壓強(qiáng)降低過(guò)程中,在體系局部形成強(qiáng)力鏈的過(guò)程的物理機(jī)制值得進(jìn)一步研究.

圖2 聲波在光滑顆粒物質(zhì)中的傳播速度與壓強(qiáng)的依賴關(guān)系,實(shí)線為擬合曲線 (a)橫波速度與壓強(qiáng)的冪律關(guān)系;(b)縱波速度與壓強(qiáng)的冪律關(guān)系;(c)剪切模量與體積模量的比值G/B與壓強(qiáng)P的變化關(guān)系Fig.2.Pressure dependence of velocity of acoustic wave propagating through frictionless granular matter,the solid line are the fitted line:(a)The velocity of transversal wave versus pressure exhibits picewise power law scaling behaviors;(b)the velocity of longitudinal wave versus pressure shows similar picewise power law as that of transversal wave;(c)the ratio of shear modulus and bulk modulus G/B versus pressure P.

材料剪切模量與體積模量的比值G/B是反映材料塑性的重要物理量,例如:對(duì)于液體,G/B=0;對(duì)于赫茲接觸固體,G/B～P1/3.為了探究卸載過(guò)程中2D顆粒體系塑性轉(zhuǎn)變,根據(jù)進(jìn)一步得到2D顆粒體系的圖2(c)顯示了2D顆粒體系剪切模量與體積模量比值G/B隨壓強(qiáng)P的變化規(guī)律.由圖2(c)可知,在本文研究的低壓強(qiáng)范圍內(nèi),2D顆粒體系的G/B隨P的變化也出現(xiàn)了類(lèi)似于在3D顆粒體系中[7,23]發(fā)現(xiàn)的隨壓強(qiáng)P增加G/B減小的現(xiàn)象,G/B～P?0.502,而與線彈性接觸理論[5,20]預(yù)言的在P＞Pj情況下G/B～P1/2不同.有趣的是,Wang等[7]的數(shù)值研究發(fā)現(xiàn)在卸載過(guò)程中硬球膠體玻璃分別經(jīng)歷高壓強(qiáng)的TL玻璃態(tài)到低壓強(qiáng)的L玻璃態(tài)及從L玻璃態(tài)到未堵塞態(tài)的兩次轉(zhuǎn)變,并且發(fā)現(xiàn)當(dāng)系統(tǒng)處在硬球玻璃的等效區(qū)域L玻璃區(qū)域,即Pg＜P＜Pj(Pg是玻璃化轉(zhuǎn)變點(diǎn))時(shí),赫茲接觸顆粒體系G/B～P?1/2,與張攀等[23]的實(shí)驗(yàn)結(jié)果及本文在2D顆粒體系的低壓強(qiáng)范圍內(nèi)中發(fā)現(xiàn)的G/B～P?0.502接近,暗示在低壓強(qiáng)范圍內(nèi)2D顆粒體系也處于L玻璃態(tài).值得一提的是,2D顆粒體系和3D顆粒體系在L玻璃態(tài)的中G/B隨壓強(qiáng)的標(biāo)度指數(shù)與顆粒作用勢(shì)無(wú)關(guān),具體原因有待進(jìn)一步研究.

3.2 衰減系數(shù)隨頻率的變化

通常,非均勻介質(zhì)的聲學(xué)衰減分為吸收衰減αa(ω)=Aω和由于聲阻抗不連續(xù)發(fā)生在顆粒邊界和缺位的彈性散射αs(ω).其中αs(ω)隨頻率的冪律關(guān)系由散射機(jī)制決定[25],例如:對(duì)于瑞利散射,即入射波長(zhǎng)頻率λ遠(yuǎn)大于散射體的平均尺寸D(λ ?D),αs(f)=S1D3ω4;對(duì)于隨機(jī)散射(λ≤D),αs(f)=S2Dω2;對(duì)于幾何散射(λ?D),αs(f)=S3/D.為了研究2D顆粒體系的聲學(xué)衰減機(jī)制,數(shù)值測(cè)量了壓強(qiáng)P=0.2,0.1,0.02,0.01,0.002,0.001,0.0002,0.0001時(shí)2D顆粒體系的橫波和縱波衰減系數(shù)隨頻率的變化,如圖3(a)和圖3(b)所示.為了避開(kāi)系統(tǒng)共振的影響,本文對(duì)遠(yuǎn)離共振峰的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了擬合,如圖3中的實(shí)線所示.由圖3(a)和圖3(b)可知,2D顆粒物質(zhì)的橫波、縱波衰減系數(shù)均隨頻率呈分段規(guī)律:當(dāng)f＜0.05時(shí),橫波、縱波的衰減系數(shù)均不隨頻率的增加而變化;當(dāng)0.05＜f＜0.35時(shí),橫縱波衰減系數(shù)與頻率成正比αT～f,αL～f,意味著在0.05—0.35頻率范圍(在P=0.02壓強(qiáng)下,對(duì)應(yīng)橫波波長(zhǎng)2.11—15.464,縱波波長(zhǎng)9.726—64.895)內(nèi),2D顆粒物質(zhì)中的衰減的主要機(jī)制為吸收衰減;當(dāng)f＞0.35時(shí),αT～f2,αL～f1.5,此時(shí)衰減機(jī)制中除了吸收衰減還有瑞利散射衰減[25].值得一提的是,在f＞0.35范圍內(nèi),橫波和縱波衰減隨頻率的變化不同(αT～f2,αL～f1.5),可能與橫波和縱波波長(zhǎng)不同導(dǎo)致散射衰減占總衰減的比例不同有關(guān)[26].

圖3 壓強(qiáng)P=0.2,0.1,0.02,0.01,0.002,0.001,0.0002,0.0001時(shí),衰減系數(shù)α與頻率f的變化 (a)橫波的衰減系數(shù)αT隨頻率f的變化;(b)縱波的衰減系數(shù)αL隨頻率f的變化;圖中用實(shí)線顯示不同壓強(qiáng)范圍內(nèi)衰減系數(shù)隨頻率的冪律關(guān)系Fig.3.The attenuation coefficient α versus frequency f at P=0.2,0.1,0.02,0.01,0.002,0.001,0.0002,0.0001:(a)The attenuation coefficient of transversal acoustic αTversus frequency f;(b)the attenuation coefficient of longitudinal acoustic αLversus frequency f.The solid line in the figure shows the picewise power law frequency dependence of attenuation coefficient.

3.3 衰減系數(shù)隨壓強(qiáng)的變化

衰減系數(shù)的標(biāo)度分析能夠提供對(duì)jammed固體結(jié)構(gòu)的洞察力.為了研究2D顆粒體系中散射特征長(zhǎng)度隨壓強(qiáng)的演化,本文測(cè)量了卸載過(guò)程2D顆粒物質(zhì)橫波、縱波長(zhǎng)波長(zhǎng)衰減系數(shù)隨壓強(qiáng)的變化,如圖4(a)和圖4(b)所示.由圖4(a)和圖4(b)可知,2D顆粒物質(zhì)的橫波、縱波衰減系數(shù)均隨壓強(qiáng)也呈分段現(xiàn)象.當(dāng)P＜10?4時(shí),橫縱波的衰減系數(shù)αT,αL幾乎不隨壓強(qiáng)變化;而當(dāng)P＞10?4時(shí),橫縱波的衰減系數(shù)αT,αL均隨壓強(qiáng)的增大而呈冪律減小,αT～P?0.099,αL～P?0.105.2D顆粒系統(tǒng)中聲衰減系數(shù)αT,αL隨壓強(qiáng)的變化與系統(tǒng)內(nèi)部結(jié)構(gòu)隨壓強(qiáng)調(diào)整有關(guān).即隨著壓強(qiáng)的增大,顆粒的平均接觸數(shù)增加,顆粒體系內(nèi)部結(jié)構(gòu)更不容易發(fā)生調(diào)整,使得橫波和縱波的聲學(xué)耗散相應(yīng)減小.這一點(diǎn)與Hong[27]對(duì)一維水平干顆粒能量衰減的研究及Liu等[28,29]、張攀等[23]對(duì)顆粒3D體系衰減耗散的研究定性一致.注意到彈性散射對(duì)衰減系數(shù)的貢獻(xiàn)與散射平均自由程ls(ω)成反比,α(ω)～1/ls,在顆粒介質(zhì)中,觀察到聲音衰減由橫向波支配.根據(jù)瑞利定律,ls(ω)～(vt/ω)4(1/D)3(其中,D是表征對(duì)散射負(fù)責(zé)的無(wú)序結(jié)構(gòu)的長(zhǎng)度尺度),可以得到α(ω)～(ω/vt)4D3.P＜10?4時(shí),vt∝P?0.446,P＞10?4時(shí),vt～P0.202,αT～P?0.099,將測(cè)得的2D顆粒體系的橫波聲速,代入α(ω)～(ω/vt)4D3式,可進(jìn)一步得到與散射相關(guān)特征長(zhǎng)度隨壓強(qiáng)的標(biāo)度:P＜10?4時(shí),??～D～P?0.595;P＞10?4時(shí),??～D～P0.236.在低壓強(qiáng)下的結(jié)果與文獻(xiàn)[5,14]中長(zhǎng)波長(zhǎng)的剪切波在J點(diǎn)附近存在??～P?1/2并在J點(diǎn)處發(fā)散的結(jié)果一致.綜合考慮圖2中vt(P)和vl(P)曲線的分段行為(P=10?4為拐點(diǎn)),我們認(rèn)為這里發(fā)現(xiàn)的與散射相關(guān)特征長(zhǎng)度隨壓強(qiáng)的分段標(biāo)度可能與vt(P)和vl(P)曲線的分段行為具有相同的結(jié)構(gòu)起源,即在卸載(整體壓強(qiáng)下降)過(guò)程中在振動(dòng)源顆粒和接收顆粒之間形成了一個(gè)強(qiáng)力鏈有關(guān).對(duì)比圖4(a)和圖4(b)可知,橫波衰減隨壓強(qiáng)變化的冪指數(shù)與縱波衰減隨壓強(qiáng)的變化的冪指數(shù)基本一致,說(shuō)明縱波與橫波的衰減行為并無(wú)太大不同.

圖4 衰減系數(shù)隨壓強(qiáng)的依賴關(guān)系 (a)在1×10?4—0.2的壓強(qiáng)范圍內(nèi),橫波衰減系數(shù)隨壓強(qiáng)冪律減小αT～P?0.099,實(shí)線為冪律擬合曲線;(b)在1×10?4—0.2的壓強(qiáng)范圍內(nèi),縱波衰減系數(shù)隨壓強(qiáng)的冪律依賴αL～P?0.105,實(shí)線為冪律擬合曲線Fig.4.The pressure dependence of attenuation coeffi-cient,the line is a power law fit:(a)The power law decrease of attenuation coefficient of transversal acoustic αTwith the pressure P increasing from 1×10?4to 0.2,αT～P?0.099;(b)the power law decrease of attenuation coefficient of longitudinal acoustic αLwith pressure P increasing from 1×10?4to 0.2,αL～P?0.105.

3.4 非線性系數(shù)隨壓強(qiáng)的變化

顆粒材料是典型的非均勻系統(tǒng),體系內(nèi)部的弱接觸使得體系表現(xiàn)出強(qiáng)烈的非線性.由于系統(tǒng)非線性的影響,當(dāng)用正弦信號(hào)激勵(lì)顆粒介質(zhì)時(shí),顆粒介質(zhì)中聲波除了一階聲波外,還會(huì)出現(xiàn)二階甚至更高階的諧波.非線性系數(shù)（其中,A1是接收基頻振幅A2是接收二倍頻的振幅,v代表聲速,ω代表頻率)是一個(gè)表征非線性大小的重要參數(shù),非線性系數(shù)越大,二次諧波成分越大,波形畸變?cè)絿?yán)重.為了研究壓強(qiáng)變化對(duì)系統(tǒng)非線性的影響,本文數(shù)值測(cè)量了2D顆粒體系橫波、縱波非線性系數(shù)隨壓強(qiáng)的變化,如圖5(a)和圖5(b)所示.由圖5(a)和圖5(b)可知,2D顆粒體系橫波、縱波二倍頻振幅與基頻振幅平方的比值均隨壓強(qiáng)呈分段關(guān)系:當(dāng)P＜10?5時(shí),橫縱波的不隨壓強(qiáng)變化;當(dāng)P＞10?5時(shí),橫波、縱波二倍頻振幅與基頻振幅平方的比值均隨壓強(qiáng)增大冪律減小:利用本文測(cè)得的可進(jìn)一步得到橫波、縱波非線性系數(shù)隨壓強(qiáng)的冪律標(biāo)度:當(dāng)P＜10?4時(shí),橫波非線性系數(shù)βT∝P?0.230,縱波則與壓強(qiáng)無(wú)明顯的依賴關(guān)系.當(dāng)P＞10?4時(shí),兩者均隨壓強(qiáng)增大呈冪律減小βT～P?0.703,βL～P?0.684.顆粒體系非線性系數(shù)隨壓強(qiáng)增加而減小可能與體系結(jié)構(gòu)隨壓強(qiáng)調(diào)整有關(guān),即隨著壓強(qiáng)的增加,顆粒間的接觸緊密,隨著接觸網(wǎng)絡(luò)的變化,非線性系數(shù)減小.值得注意的是,2D顆粒體系表現(xiàn)出與3D赫茲接觸顆粒體系[30]類(lèi)似的非線性隨壓強(qiáng)的冪律衰減.但是,當(dāng)P＞10?4時(shí),2D體系的非線性系數(shù)隨壓強(qiáng)變化的冪指數(shù)的絕對(duì)值略大于3D體系β～P?0.51的冪指數(shù).

需要強(qiáng)調(diào)的是,以單一顆粒激勵(lì)作為聲源可能會(huì)測(cè)量到力鏈效應(yīng).本文中測(cè)量到的vt(P)和vl(P)曲線以及與散射相關(guān)特征長(zhǎng)度隨壓強(qiáng)變化曲線的分段標(biāo)度很可能與卸載(整體壓強(qiáng)下降)過(guò)程中在振動(dòng)源顆粒和接收顆粒之間形成了一個(gè)強(qiáng)力鏈有關(guān).另外,由單粒子源輻射引起的擴(kuò)散衰減的削減可能與在源顆粒和接收顆粒之間存在的多個(gè)衰減通路的統(tǒng)計(jì)平均有關(guān).

圖5 (a)橫波非線性系數(shù)隨壓強(qiáng)的變化,在低壓強(qiáng)下βT～P?0.230,在高壓強(qiáng)下βT～P?0.703,實(shí)線為擬合曲線;(b)縱波非線性系數(shù)隨壓強(qiáng)的變化,在高壓強(qiáng)下βL～P?0.684;插圖為橫縱波的二倍頻振幅與基頻振幅平方的比值隨壓強(qiáng)P的變化關(guān)系,實(shí)線為擬合曲線Fig.5.(a)The pressure dependence of the nonlinear coeffcient of transverse wave,βT～P?0.230,at lower pressure,βT～P?0.703,at higher pressure,and the solid line is a power law fit;(b)the pressure dependence of the nonlinear coefficient of longitudinal wave,βL～P?0.684,at higher pressure.The inset is the ratio of the second harmonic amplitude of the transverse and longitudinal wave to the square of the fundamental frequency amplitudewith pressure P,and the line is a power law fit.

4 結(jié) 論

數(shù)值研究了在二維光滑圓盤(pán)顆粒體系中傳播過(guò)程的橫波、縱波的聲速、剪切模量與體積模量的比值G/B、聲衰減及非線性效應(yīng)隨壓強(qiáng)的變化規(guī)律,得到如下結(jié)論.

1)2D顆粒體系中,橫波和縱波聲速隨壓強(qiáng)的變化均出現(xiàn)了分段現(xiàn)象,當(dāng)P＜10?4時(shí),橫縱波聲速均隨壓強(qiáng)的增大而減小,vt～P?0.446,vl～P?0.232;而當(dāng)P＞10?4時(shí),橫波、縱波聲速均隨壓強(qiáng)呈現(xiàn)冪律關(guān)系為vt～P0.202,vl～P0.338.由此可以得到2D圓盤(pán)顆粒體的G/B隨壓強(qiáng)P的變化關(guān)系G/B～P?0.502對(duì)應(yīng)著L玻璃態(tài).

2)2D顆粒體系中橫波和縱波衰減隨頻率的變化規(guī)律不同:對(duì)于橫縱波,在低頻率下,有α～f,暗示2D顆粒體系的橫縱向衰減主要是起因于阻尼導(dǎo)致的吸收衰減及少部分散射衰減.而在高頻率下,橫波經(jīng)歷了αT～f到αT～f2的轉(zhuǎn)變,縱波經(jīng)歷了αL～f到αL～f1.5的轉(zhuǎn)變.表明隨頻率的增加,橫縱波衰減中的散射衰減成分增加.

3)2D顆粒體系的聲衰減系數(shù)隨壓強(qiáng)也呈分段規(guī)律:當(dāng)P＜10?4時(shí),橫波、縱波聲衰減系數(shù)αT和αL幾乎都不隨壓強(qiáng)變化;當(dāng)P＞10?4時(shí),橫波、縱波聲衰減系數(shù)與壓強(qiáng)呈現(xiàn)冪律規(guī)律,αT～P?0.099,αL～P?0.105,進(jìn)而得到2D顆粒體系散射相關(guān)的特征長(zhǎng)度??隨壓強(qiáng)P的冪律關(guān)系:當(dāng)P＜10?4時(shí)??～P?0.595;當(dāng)P＞10?4時(shí)??～P0.236.

4)2D顆粒體系的橫波的非線性響應(yīng)要略大于縱波,且橫波、縱波非線性系數(shù)隨壓強(qiáng)的變化也呈分段現(xiàn)象:當(dāng)P＜10?4時(shí),橫波非線性系數(shù)βT～P?0.230,βL不隨壓強(qiáng)變化;而當(dāng)P＞10?4時(shí),橫波、縱波非線性系數(shù)均隨壓強(qiáng)P的增加呈冪律減小,βT～P?0.703,βL～P?0.684.在高壓強(qiáng)下非線性系數(shù)的減小暗示高壓強(qiáng)下顆粒間的接觸緊密,沿強(qiáng)力鏈傳播的聲波產(chǎn)生了更少的非線性響應(yīng).

[1]Liu A J,Nagel S R 1998Nature396 21

[2]Henkes S,Chakraborty B 2005Phys.Rev.Lett.95 198002

[3]O’Hern C S,Silbert L E,Nagel S R 2003Phys.Rev.E68 011306

[4]O’Hern C,Langer S,Liu A,Nagel S 2002Phys.Rev.Lett.88 075507

[5]Xu N 2011Front.Phys.6 109

[6]Ikeda A,Berthier L 2015Phys.Rev.E92 012309

[7]Wang X,Zheng W,Wang L,Xu N 2015Phys.Rev.Lett.114 035502

[8]Coulais C,Behringer R P,Dauchot O 2014Soft Matter10 1519

[9]Karimi K,Maloney C E 2015Phys.Rev.E92 022208

[10]Sussman D M,Goodrich C P,Liu A J,Nagel S R 2015Soft Matter11 2745

[11]Wyart M,Nagel S R,Witten T A 2005Europhys.Lett.72 486

[12]Wyart M,Silbert L,Nagel S 2005Phys.Rev.E72 051306

[13]Silbert L,Liu A,Nagel S 2005Phys.Rev.Lett.95 098301

[14]Vitelli V 2010Soft Matter6 3007

[15]Mizuno H,Silbert L E,Sperl M 2016Phys.Rev.Lett.116 068302

[16]Merkel A,Tournat V,Gusev V 2014Phys.Rev.E90 023206

[17]Zhang Q,Li Y,Hou M,Jiang Y,Liu M 2012Phys.Rev.E85 031306

[18]Zhang Q,Li Y C,Liu R,Jiang Y M,Hou M Y 2012Acta Phys.Sin.61 234501(in Chinese)[張祺,李寅閶,劉銳,蔣亦民,厚美瑛2012物理學(xué)報(bào)61 234501]

[19]Zheng H P 2014Chin.Phys.B23 054503

[20]Vitelli V,Xu N,Wyart M,Liu A J,Nagel S R 2010Phys.Rev.E81 021301

[21]Somfai E,Roux J N,Snoeijer J,van Hecke M,van Saarloos W 2005Phys.Rev.E72 021301

[22]Jia X,Caroli C,Velicky B 1999Phys.Rev.Lett.82 1863

[23]Zhang P,Zhao X D,Zhang G H,Zhang Q,Sun Q C,Hou Z J,Dong J J 2016Acta Phys.Sin.65 024501(in Chinese)[張攀,趙雪丹,張國(guó)華,張祺,孫其誠(chéng),侯志堅(jiān),董軍軍2016物理學(xué)報(bào)65 024501]

[24]Lherminier S,Planet R,Simon G,Vanel L,Ramos O 2014Phys.Rev.Lett.113 098001

[25]West B J,Shlesinger M F 1984J.Stat.Phys.36 779

[26]Langlois V,Jia X 2015Phys.Rev.E91 022205

[27]Hong J 2005Phys.Rev.Lett.94 108001

[28]Wang P J,Li Y D,Xia J H,Liu C S 2008Phys.Rev.E77 060301

[29]Wang P J,Xia J H,Li Y D,Liu C S 2007Phys.Rev.E76 041305

[30]Brunet T,Jia X,Johnson P A 2008Geophys.Res.Lett.35 L19308

Numerical study on acoustic behavior of two-dimensional granular system?

Liu Xiao-Yu1)Zhang Guo-Hua1)?Sun Qi-Cheng2)?Zhao Xue-Dan1)Liu Shang1)

1)(Department of Physics,University of Science and Technology Beijing,Beijing 100083,China)
2)(State Key Laboratory of Hydroscience and Engineering,Tsinghua University,Beijing 100084,China)

15 May 2017;revised manuscript

17 July 2017)

The transversal and longitudinal wave velocities,the acoustic attenuation coefficients,the nonlinear coefficients at different pressures and the acoustic attenuation coefficient as a function of frequency in a two-dimensional(2D)monodisperse disc system are numerically calculated.The results show that the transversal and longitudinal wave velocities both exhibit a piecewise power law with pressureP.WhenP＜10?4,the velocity decreases with the increase of pressure in the 2D disc granular system,and whenP＞10?4,the transversal wave velocityVtand longitudinal wave velocityVlshow the scaling power laws,i.e.,vt～P0.202andvl～P0.338,respectively.The ratio of the shear modulus to the bulk modulusG/Bshows a power law scaling with the pressure,G/B～P?0.502,implying that the system lies in an L glass state at low pressure,similar to that of a three-dimensional(3D)spherical granular system.The attenuation coefficients(αT,αL)of the horizontal excitation and vertical excitation also show the picecewise behaviors with the change of frequencyf.Whenf＜0.05,neither of the two attenuation coefficients changes with frequencyf.Whenf＞0.05,α∝fαT,αL∝f.And whenf＞0.35,αT∝f2andαL∝f1.5.In addition,the nonlinear coefficient and the attenuation coefficient of the 2D disc granular system under the vertical and horizontal excitation both also show a piecewise law behavior with pressure,similar to that of the acoustic velocity.WhenP＜10?4,only the transversal nonlinear coefficient changes according toβT∝P?0.230,while the other coefficient has no change.WhenP＞10?4,the attenuation coefficients and nonlinear coefficients decrease according to their power law with the increase of pressure,i.e.,βT∝P?0.703,βL∝P?0.684,αT∝P?0.099,αL∝P?0.105.The characteristic length??,which characterizes the disordered structure responsible for the scattering,also decreases according to power law with the increase of pressure,whenP＜10?4,?? ∝P?0.595;whenP＞10?4,?? ∝P0.236.

granular matter,acoustic velocity,nonlinear,acoustic attenuation

PACS:45.70.–n,46.40.Cd,43.40.+sDOI:10.7498/aps.66.234501

*Project supported by the National Natural Science Foundation of China(Grant Nos.11272048,11572178,91634202).

?Corresponding author.E-mail:zhguohua@sas.ustb.edu.cn

?Corresponding author.E-mail:qcsun@tsinghua.edu.cn

(2017年5月15日收到;2017年7月17日收到修改稿)

數(shù)值測(cè)量了卸載過(guò)程中二維單分散圓盤(pán)顆粒系統(tǒng)的橫波、縱波聲速、聲衰減系數(shù)、非線性系數(shù)隨壓強(qiáng)的變化以及聲衰減系數(shù)隨頻率的變化.結(jié)果表明,二維(2D)圓盤(pán)顆粒體系的橫波、縱波聲速均隨壓強(qiáng)呈分段冪律標(biāo)度:當(dāng)壓強(qiáng)P＜10?4時(shí),橫波、縱波聲速隨壓強(qiáng)的增大而減小;當(dāng)P＞10?4時(shí),有vt～P0.202,vl～P0.338.進(jìn)一步得到其剪切模量和體積模量的比值G/B也隨壓強(qiáng)呈冪律標(biāo)度,G/B～P?0.502,暗示在低壓強(qiáng)下,與三維(3D)球形顆粒體系類(lèi)似,2D圓盤(pán)顆粒體系也處于L玻璃態(tài).水平激勵(lì)和垂直激勵(lì)下2D圓盤(pán)顆粒系統(tǒng)的衰減系數(shù)隨頻率變化也呈現(xiàn)分段行為:當(dāng)頻率f＜0.05時(shí),衰減系數(shù)不隨f變化;當(dāng)f＞0.05時(shí),橫波縱波的衰減系數(shù)α～f;當(dāng)f＞0.35時(shí),橫波衰減系數(shù)αT～f2,縱波衰減系數(shù)αL～f1.5.此外,豎直水平激勵(lì)下的2D圓盤(pán)顆粒系統(tǒng)的非線性系數(shù)和衰減系數(shù)隨壓強(qiáng)也呈現(xiàn)與聲速類(lèi)似的分段規(guī)律:當(dāng)P＜10?4時(shí),橫波非線性系數(shù)βT～P?0.230,其余都不隨壓強(qiáng)變化.當(dāng)P＞10?4時(shí),兩者均隨壓強(qiáng)增大呈冪律減小:βT～P?0.703,βL～P?0.684,αT～P?0.099,αL～P?0.105.進(jìn)而得到2D圓盤(pán)顆粒系統(tǒng)中散射相關(guān)的特征長(zhǎng)度??隨壓強(qiáng)呈冪律標(biāo)度,當(dāng)P＜10?4時(shí),??～P?0.595;當(dāng)P＞10?4時(shí),??～P0.236.

10.7498/aps.66.234501

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):11272048,11572178,91634202)資助的課題.

?通信作者.E-mail:zhguohua@sas.ustb.edu.cn

?通信作者.E-mail:qcsun@tsinghua.edu.cn